Empirical Green's tensor retrieved from ambient noise cross-correlations at The Geysers geothermal field, Northern California

Empirical Green's tensor retrieved from ambient noise cross-correlations at The Geysers... Summary We retrieve empirical Green's functions in the frequency range (∼0.2–0.9 Hz) for interstation distances ranging from ∼1 to ∼30 km (∼0.22 to ∼6.5 times the wavelength) at The Geysers geothermal field, Northern California, from coherency of ambient seismic noise being recorded by a variety of sensors (broad-band, short-period surface and borehole sensors, and one accelerometer). The applied methodology preserves the intercomponent relative amplitudes of the nine-component Green's tensor that allows us to directly compare noise-derived Green's functions (NGFs) with normalized displacement waveforms of complete single-force synthetic Green's functions (SGFs) computed with various 1-D and 3-D velocity models using the frequency-wavenumber integration method and a 3-D finite-difference wave propagation method, respectively. These comparisons provide an effective means of evaluating the suitability of different velocity models to different regions of The Geysers, and assessing the quality of the sensors and the NGFs. In the T-Tangential, R-Radial, Z-Vertical reference frame, the TT, RR, RZ, ZR and ZZ components (first component: force direction, second component: response direction) of NGFs show clear surface waves and even body-wave phases for many station pairs. They are also broadly consistent in phase and intercomponent relative amplitudes with SGFs for the known local seismic velocity structure that was derived primarily from body-wave traveltime tomography, even at interstation distances less than one wavelength. We also find anomalous large amplitudes in TR, TZ, RT and ZT components of NGFs at small interstation distances (≲4 km) that can be attributed to ∼10°–30° sensor misalignments at many stations inferred from analysis of longer period teleseismic waveforms. After correcting for sensor misalignments, significant residual amplitudes in these components for some longer interstation distance (≳8 km) paths are better reproduced by the 3-D velocity model than by the 1-D models incorporating known values and fast axis directions of crack-induced VS anisotropy in the geothermal field. We also analyse the decay of Fourier spectral amplitudes of the TT component of NGFs at 0.72 Hz with distance in terms of geometrical spreading and attenuation. While there is considerable scatter in the NGF amplitudes, we find the average decay to be consistent with the decay expected from SGF amplitudes and with the decay of tangential component local-earthquake ground-motion amplitudes with distance at the same frequency. Seismic anisotropy, Seismic interferometry, Seismic noise, Surface waves and free oscillations, Wave propagation, Crustal structure INTRODUCTION The Geysers in the Mayacamas Mountains of Northern California are the largest complex of geothermal power plants in the world today. Spread over an area of ∼115 km2, it currently produces ∼835 MW of electricity (http://www.energy.ca.gov/tour/geysers/, last accessed August 2016). It is a vapour-dominated geothermal reservoir (Allis & Shook 1999) in which approximately 20 million gallons of reclaimed wastewater are injected on a daily basis (http://www.geysers.com/water.aspx, last accessed August 2016) from neighbouring communities (Lake County, since 1998; City of Santa Rosa, since 2004). The water is injected to tap the heat in the volcanic reservoir rocks that are at ∼240 to ∼300° C (Truesdell et al.1991; Lowenstern & Janik 2003; Majer & Peterson 2007). While the injection of wastewater has successfully replenished the resource and the steam production that was in severe decline during the late 1980s, it has also been linked to elevated rates of microseismicity (MW ≤ 3) in the reservoir when compared to regional rates (Ludwin et al.1982). Many studies have demonstrated strong temporal correlation of microseismicity with both steam production and water injection at The Geysers, on both localized and field-wide scales (Oppenheimer 1986; Stark 1991; Majer & Peterson 2007; Johnson et al.2016; Trugman et al.2016). Detailed analysis of microseismicity has also revealed significant temporal variations in earthquake hypocenter distributions, in b-value and in the number of large magnitude (MW > 4) earthquakes (Trugman et al.2016). Investigations of focal mechanisms and moment tensors of these earthquakes using body-wave polarities, amplitude ratios, and low-frequency displacement waveforms have shed some light on various aspects of their source mechanisms (Oppenheimer 1986; Ross et al.1999; Guilhem et al.2014; Boyd et al.2015). However, imperfect knowledge of subsurface seismic velocity structure and assumption of 1-D velocity models are believed to introduce considerable uncertainties in synthetic ray paths (Oppenheimer 1986) and in synthetic Green's functions (SGFs) required for source inversion. 3-D P and S-wave traveltime tomography studies have revealed considerable heterogeneities in the subsurface (Julian et al.1996; Ross et al.1999; Gritto et al.2013a) that must be taken into account in the inversion of polarities and seismic waveforms for recovery of accurate source mechanisms. In addition to the temporal changes in microseismicity, temporal changes in bulk medium properties of the reservoir such as fluid pressure and VP/VS ratio have also been observed from levelling and GPS studies, time-lapse body-wave travel time tomography, etc. (Foulger et al.1997; Mossop & Segall 1999; Gritto et al.2013b). The response of a thermal reservoir to water injection or production, generally indicated by (1) the spatiotemporal behaviour and mechanisms of seismicity, and/or (2) changes in seismic properties of the subsurface, is of great interest to further our understanding of induced seismicity and geothermal systems. This knowledge can help to better characterize and manage the seismic hazard by optimizing fluid injection into the subsurface for energy or storage purposes. In this study, we establish a framework for ambient seismic noise cross-correlation at The Geysers, which will aid in improving our understanding of the seismic velocity structure. The cross-correlations of ambient seismic noise obtained at pairs of seismic sensors have been shown to converge over a period of time to the empirical Green's functions (GFs) or the seismic medium's response at one receiver to a unit force applied at the other receiver (Shapiro & Campillo 2004; Shapiro et al.2005; Bensen et al.2007). In numerous studies, the body and surface waves in the noise-derived Green's functions (NGFs) have been used for tomography at various scales (e.g. Shapiro et al.2005; Lin et al.2008, 2014; Nakata et al.2015). Additionally, the coda of NGFs can be used for spatiotemporal monitoring of small subsurface seismic velocity changes (e.g. Sens-Schönfelder & Wegler 2006; Brenguier et al.2014; Taira et al.2015) that have previously been observed at The Geysers with body-wave tomography. However, the non-uniformity in ambient noise source distribution results in systematic differences between recovered NGFs and the exact GFs in the real Earth (Tsai 2009). This can lead to errors in the estimates of subsurface seismic properties inferred from NGFs and limits their utility and reliability. In this study, we retrieve NGFs, represented by interstation coherency (Prieto et al.2011), in the frequency range of ∼0.2–0.9 Hz for interstation distances from ∼1 to ∼30 km at the Geysers. We use seismic data that are continuously recorded by a wide variety of sensors in and around the geothermal reservoir. The minimum VS is ∼2 km s−1 and the interstation distances range from ∼0.22 λ (near-field) to ∼6.5 λ (far-field) assuming a phase velocity ∼2.3 km s−1 at a frequency of ∼0.5 Hz, where λ is the frequency-dependent surface-wave wavelength and distances ≳1λ–2λ are assumed to be far-field. For pairs of three-component sensors (T: Tangential, R: Radial, Z: Vertical), we retrieve nine components of NGFs (TT, TR, etc. where the first component is direction of the applied force and the 2nd component is the displacement direction). We directly compare NGF waveforms to normalized single-force displacement SGFs, computed using various 1-D velocity models utilized for various monitoring applications at The Geysers. We evaluate the similarity of NGF and SGF waveforms in terms of waveform fits, phase and relative intercomponent amplitudes. NGFs with little or no coherent energy in the ballistic wave arrival time window or NGFs that bear little or no resemblance to SGFs expected for a wide variety of reasonable and expected velocity models are interpreted to be contaminated with errors possibly due to poor sensor coupling, or non-uniformity of ambient noise source distribution. After discarding these erroneous NGFs, we use the remaining NGFs to evaluate the applicability of various 1-D velocity models to different sub-regions of The Geysers as NGFs contain information about real Earth 3-D wave propagation (Ma et al.2008). We examine if NGFs for station pairs in a particular sub-region are systematically better fit by SGFs of a particular velocity model compared to SGFs of other velocity models. We interpret the phase differences between NGFs and SGFs in terms of the deficiencies in the velocity models. While some phase errors in the NGFs are expected due to the non-uniformity of ambient noise source distribution, our interstation paths are well distributed in azimuth, which allows us to draw inferences from spatially coherent phase differences pervasive over multiple azimuths. We show that waveforms, phases and intercomponent relative amplitudes of the primary non-zero components of retrieved NGFs (i.e. TT, RR, RZ, ZR, ZZ) are similar to those of SGFs, even at distances < λ. We compare our inferences on the features of subsurface seismic velocities at The Geysers from the NGF-SGF comparisons with the inferences drawn in many previous body-wave traveltime studies. We also examine if the NGFs compare better to SGFs computed using a 3-D velocity model of the geothermal field than with SGFs from the 1-D velocity models. Our study additionally includes: (1) an analysis of long period teleseismic waveforms to verify sensor orientations, (2) application of the Optimal Rotation Algorithm (Roux 2009) on the NGFs to detect any dominant noise source illumination direction, and (3) comparison with SGFs computed using the 1-D models incorporating known values and fast axis directions of crack-induced VS anisotropy in the geothermal field. We evaluate the contribution of the above-mentioned factors to the significant non-zero amplitudes in TR, TZ, RT and ZT components of NGFs for many stations pairs. We also compare the decay of Fourier spectral amplitudes of the TT component of NGFs at 0.72 Hz as a function of distance with the decay of synthetic spectral amplitudes for velocity models incorporating strong, weak and known values of anelastic attenuation, and with the decay of tangential component local-earthquake ground-motion amplitudes. BACKGROUND AND MOTIVATION Noise cross-correlation tomography studies typically involve phase velocity measurements on the fundamental-mode Rayleigh waves obtained from cross-correlation of ambient noise recorded on vertical components of broad-band stations at far-field interstation distances (e.g. Yao et al.2006; Lin et al.2008, 2014). The theoretical proof for retrieval of ballistic surface waves from cross-correlation of diffused waves at two receivers was first provided by Snieder (2004) under a stationary phase approximation valid at far-field distances (Yokoi & Margaryan 2008; Zhan & Ni 2010). Subsequently, using source-receiver reciprocity, Wapenaar & Fokkema (2006) derived an expression for the velocity GF response at a receiver to a single force applied at another receiver in an inhomogeneous anisotropic loss-less medium from cross-correlation of responses of sources distributed over an enclosing surface. The variations of their basic formulation involve variations in the nature of the boundary and the outside medium, types of sources, types of responses cross-correlated or retrieved, distance of the enclosing source surface from the receivers and whether the cross-correlation operates on the responses of individual sources or the net ground motion field at the receivers. In practice, the application of most temporal and spectral normalization methods removes the absolute amplitude information in component-pair NGFs leading to their ‘empirical’ nature (Bensen et al.2007). However, significant relative amplitude information in addition to the phase information can still be retrieved from NGFs, although many aspects are still under considerable debate as discussed in the Amplitude decay section. Cross-correlation of spectrally whitened noise is equivalent to computing interstation coherency (Bensen et al.2007; Prieto et al.2011). For components i, j of velocity v recorded at stations A, B at positions xA, xB respectively, coherency γij(xA, xB, ω) as a function of frequency ω is defined as   \begin{eqnarray} - {G_{ij}}\left( {{x_A},{x_B},\omega } \right)& \approx & {\gamma _{ij}}\ \left( {{x_A},{x_B},\omega } \right) \nonumber\\ &=& \left\langle {\frac{{v_i^*\left( {{x_A},\omega } \right)\ {v_j}\left( {{x_B},\omega } \right)}}{{\left\{ {\left| {{v_i}\left( {{x_A},\omega } \right)} \right|} \right\}\left\{ {\left| {{v_j}\left( {{x_B},\omega } \right)} \right|} \right\}}}} \right\rangle , \end{eqnarray} (1) where Gij(xA, xB, ω) is the displacement NGF with single force, displacement directions i, j at stations A, B that are acting as source and receiver, respectively; <> implies stacking results for data recorded over multiple time windows known as ensemble averaging, {} implies spectral amplitude smoothing (three-point ∼0.02 Hz moving window average in our study). The directions i, j can be radial (R; positive outwards), transverse (T; positive clockwise 90° from radial) or vertical (Z; up positive). In eq. (1), in the following and throughout this text, we define single-force GF Gij(xA, xB, ω) as the displacement response to an input step force and directly relate it to interstation coherency obtained from the cross-correlation of spectrally whitened ambient noise. This is different than the traditional definition of a GF in which it is defined as the displacement response to an input impulsive force (Aki & Richards 2002). For the traditional definition (i.e. impulsive force as an input), the GF can be related to the time derivative of coherency as in many other studies, that is, Lin et al. (2008), Ma et al. (2008), etc. For a particular component pair, γ has been shown to preserve geometrical spreading and anelastic attenuation information as a function of interstation distance and structure (Prieto et al.2011; Lawrence et al.2013). Additionally, if the same normalization factors are used for all components of a station, relative intercomponent amplitudes for a particular station pair can be preserved enabling extraction of attributes such as the Rayleigh wave ellipticity from the amplitude ratio of R and Z components of NGFs (Lin et al.2014). The amplitude spectrum of the noise at the source station can also be used to normalize the noise spectra at both stations, which yields the impulse response function (Prieto et al.2011; Denolle et al. 2013), and is known as deconvolution interferometry (Nakata et al.2011). The impulse response function has been used to extract site response of basins and 3-D structures at long periods, ∼4–10 s (Prieto & Beroza 2008; Denolle et al.2014). The retrieved phase information is the same for different spectral normalization methods (Prieto et al.2011). While interstation distances ≳λ are considered to be ‘far-field’ for surface waves (Lin et al.2013), phase velocity measurements from NGFs are usually restricted to interstation distances ≳3λ in order to avoid bias at shorter distances caused by inhomogeneous noise source distributions (Lin et al.2008, 2014). Various numerical studies have shown that phase velocity measurements derived from NGFs at interstation distances ≲3λ–5λ can have significant error, even for seismic noise sources that are azimuthally widely distributed (Kimman & Trampert 2010; Zhan & Ni 2010). However, reliable phase velocities can be extracted from NGFs at distances down to ∼1λ using sophisticated methods that are able to reduce biases caused by inhomogeneously distributed noise sources. These methods include measurement of phase velocities by phase-front tracking on dense arrays (Lin et al.2013) and azimuthal averaging of multiple coherency measurements as in spatial autocorrelation (Tsai & Moschetti 2010). The studies on relative amplitude extraction from NGFs, mentioned earlier, have also utilized NGFs at primarily far-field distances assuming Rayleigh waves on RZ, RR, ZR and ZZ components and Love waves on TT components, which are the non-zero components of the GF tensor for 1-D isotropic media. The TR, TZ, RT and ZT components of the GF tensor can be non-zero in the presence significant 3-D structure or anisotropy. In NGFs, these components can also be non-zero due to inhomogeneous noise source distribution [e.g. Parkfield, California (Durand et al.2011)] and have not been analysed in detail in crustal structure studies [e.g. some paths in Southern California (Denolle et al.2013)]. Sensor misalignments or incorrect knowledge of sensor gains and orientations may lead to incorrectly oriented NGFs and anomalous amplitudes in the TR, TZ, RT and ZT components as well. To the best of our knowledge, the retrieval of relative amplitude information in NGFs at distances < λ has not yet been explored. METHODOLOGY Data Fig. 1 shows a map of The Geysers, the seismic stations employed in this study and the approximate outline of the steam field of the reservoir area (Gritto et al.2013a). Most of our data are recorded by the Lawrence Berkeley National Laboratory (LBNL) network at The Geysers consisting of ∼30 short-period (4.5 Hz) three-component geophones (OYO-GS11D), one Nanometrics Titan accelerometer (station DRH) and three short-period (8.0 Hz) three-component borehole sensors (OYO-GS11D8) at depths of ∼35–150 m (stations DEB, JKB and SRB). We include data recorded by short-period (1.0 Hz) vertical component sensors (L4) and by one broad-band STS2 sensor (station GDXB) at stations operated by the United States Geological Survey (USGS) in and around the geothermal field. Data recorded by the two closest broad-band stations (an STS1 at HOPS and an STS2 at MNRC) of the Berkeley Digital Seismic Network (BDSN) are also included to provide a reference for quality control of data recorded at all other stations in our study. Figure 1. View largeDownload slide Map shows the outline of the reservoir area (red polygon), faults (solid black lines) and seismic stations plotted on the gradient of topography at The Geysers. Black triangles: LBNL 4.5 Hz geophones, green squares: LBNL 8.0 Hz borehole geophones, red square: LBNL accelerometer, blue triangles: USGS 1.0 Hz vertical component sensors, and blue square: USGS broad-band sensor. Inset shows the location of The Geysers (red square) and BDSN stations (red triangles) with respect to the San Francisco Bay area in Northern California, USA. Figure 1. View largeDownload slide Map shows the outline of the reservoir area (red polygon), faults (solid black lines) and seismic stations plotted on the gradient of topography at The Geysers. Black triangles: LBNL 4.5 Hz geophones, green squares: LBNL 8.0 Hz borehole geophones, red square: LBNL accelerometer, blue triangles: USGS 1.0 Hz vertical component sensors, and blue square: USGS broad-band sensor. Inset shows the location of The Geysers (red square) and BDSN stations (red triangles) with respect to the San Francisco Bay area in Northern California, USA. Cross-correlation analysis Our data pre-processing and cross-correlation workflows are derived from Bensen et al. (2007), Seats et al. (2012) and Lin et al. (2014). Most importantly, we use the same temporal and spectral normalization factors for all components of a three -component station in order to preserve the relative amplitudes between the components (e.g. Lin et al.2014). Other details of the methodology adopted are provided in  Appendix A1. The details of the subsequent quality control procedure are described in  Appendix A2. After quality control, we average the causal and anticausal sides of the final stacked cross-correlations or coherencies, extracting the symmetric component of NGFs in the T-R-Z reference frame. Fig. 2 shows ZZ and TT components of NGFs in the interstation distance range 0–25 km in which most of our station pairs are located. The waveforms are stacked in 0.5 km interstation distance bins. Rayleigh and Love waves can easily be identified in components ZZ and TT, respectively, with apparent velocities of ∼2.5 and ∼2.8 km s−1, respectively. ZZ components at distances of ∼17–20 km also show faint but coherent body-wave or higher-mode Rayleigh wave energy that is faster than the fundamental-mode Rayleigh waves but slower than direct P-waves, traveling at ∼3.6 km s−1 (composition of waveforms further discussed in the Results section). We also estimate an empirical expression for distance-dependent approximate start time of coda waves following the passage of the surface waves by fitting a fourth degree polynomial to times at which envelopes of NGFs drop below a level of 0.2 times the peak absolute amplitude. This empirical coda-wave start time is used for selecting waveform segments that contain body-wave and surface-wave phases in order to evaluate waveform fits between NGFs and SGFs. Figure 2. View largeDownload slide NGFs filtered between ∼0.2–0.9 Hz and stacked in 0.5 km interstation distance bins. Dashed red curves are reference traveltime curves for various constant apparent velocities (shifted to the left by 0.7 s for clarity). Dashed blue curve shows the empirically estimated approximate start time of coda waves after the passage of surface waves. The NGF waveforms are proportional to displacement response for an input step force. Figure 2. View largeDownload slide NGFs filtered between ∼0.2–0.9 Hz and stacked in 0.5 km interstation distance bins. Dashed red curves are reference traveltime curves for various constant apparent velocities (shifted to the left by 0.7 s for clarity). Dashed blue curve shows the empirically estimated approximate start time of coda waves after the passage of surface waves. The NGF waveforms are proportional to displacement response for an input step force. Velocity models Here we briefly describe the velocity structure of the reservoir area. Gritto et al. (2013a) conducted a tomographic campaign to jointly invert for the 3-D P- and S-wave velocity structures and the hypocentre locations of microseismicity at The Geysers using body-wave arrival times of ∼32,000 earthquakes at the LBNL stations. The starting 3-D velocity model was based on 3-D interfaces in the reservoir (interface of top steam entry and interface of felsite) coupled with P- and S-wave velocity estimates that were derived from a vertical seismic profiling (VSP) experiment in the centre of the reservoir (O'Connell & Johnson 1991). The top of the steam interface is considered the top of the reservoir and separates mélange (a complex assemblage of Franciscan greywacke, greenstone and serpentinite) from the deeper metamorphosed reservoir rocks (Franciscan greywacke and metagraywacke) at depths of ∼1–3 km. The top of the felsite interface demarks the transition from metamorphosed reservoir rock to the underlying granitic pluton at depths greater than ∼3–5 km. On average, VP, VS values increase, respectively, from ∼3.6 and ∼2 km s−1 at the surface, to ∼5.9 and ∼3.3 km s−1 at depth of ∼6 km, with typical values of ∼4.8 and ∼2.8 km s−1 for reservoir rocks, and ∼5.5 and ∼3.0 km s−1 for the underlying felsite (O'Connell & Johnson 1991). The reservoir area is bounded by NNW–SSE trending faults, namely the Mercuryville fault to the southwest and the Collayomi fault to the northeast. There is widespread lateral heterogeneity in and around the reservoir, with up to ∼20 per cent variations in VP at shallow (less than ∼2.5 km) depths (Julian et al.1996; Ross et al.1999). On average, the velocities in the northwestern section of The Geysers have been found to be lower by ∼10 per cent than in the central and southeastern section (Eberhart-Phillips 1986; Julian et al.1996; Gritto et al.2013a). The VP/VS ratio, believed to be relatively insensitive to lithology but quite sensitive to the saturation of rocks (Gritto & Jarpe 2014) and to the compressibility of pore fluids, shows a large negative anomaly (up to ∼–9 per cent) at the centre and southeast of the centre of the reservoir (Julian et al.1996; Foulger et al.1997; Ross et al.1999). This anomaly exists at depths up to ∼2.5 km and is anticorrelated with a positive VS anomaly (Gritto et al.2013a). We select four 1-D velocity models (Fig. 3) at The Geysers to compute single force SGFs for comparison with the NGFs of all station pairs not including stations HOPS and MNRC. The models are, (1) REF, (2) BACK1, (3) VSP0, and (4) AVG1. Model REF is a slightly modified version of the 1-D average of the 3-D model of Julian et al. (1996) for the NW Geysers and was used by Guilhem et al. (2014) for modelling 0.5–2.5 Hz waveforms of MW ∼ 3.3–3.9 earthquakes at distances up to a few km. Model BACK1 is a 1-D depth average of the 3-D starting velocity model used in Gritto et al. (2013a), while model AVG1 is a 1-D depth average of their final 3-D velocity model within the reservoir area. VSP0 is the velocity profile shown in fig. 2 in Gritto et al. (2013a) and is derived from a VSP profile at the felsite peak in the SE Geysers (O'Connell & Johnson 1991). The depth-dependent density and anelastic attenuation quality factors (QP and QS) are borrowed from Guilhem et al. (2014). In our study, we ignore the surface topography at The Geysers and assume our reference datum to be ∼800 m above the mean sea level, which is the average elevation of stations that we are using. For all station pairs that include station HOPS or MNRC, we use the California Central Coast Ranges velocity model GIL7 (Stidham et al.1999) and its modified version GIL7.1 (Fig. 3) to compute SGFs (modifications explained in section NGFs with BDSN stations HOPS and MNRC). Upper crustal QP and QS values in GIL7 and its modified versions are reduced to values comparable to those of REF. Figure 3. View largeDownload slide 1-D velocity models used in our study. Figure 3. View largeDownload slide 1-D velocity models used in our study. Fig. 4 shows a few example depth slices of G3D1, a modified version of the 3-D velocity model of Gritto et al. (2013a) used to compute SGFs accounting for the 3-D structure in the reservoir in our study. G3D1 consists of a 3-D velocity structure surrounded by a 1-D background model (BACK1). The horizontal extent of the 3-D structure is approximately restricted by the outline of the reservoir steam field at the surface (red polygon in Fig. 1) (except towards the northwest where the LBNL network extends beyond the outline; see stations HER, AL2 and DRH in Fig. 1), and decreases with increasing depth. The modifications to the model are explained as follows. For nodes with VP/VS ratio $$ < \sqrt 2 $$, VP and VS values are increased and decreased respectively in small increments, until they satisfy the minimum VP/VS ratio criteria of $$\sqrt 2 $$ of our 3-D wave propagation software. The velocities are tapered at the edges of the model so that they smoothly transition to model BACK1. The model is extended laterally assuming the values of BACK1, and extended in depth assuming a half-space with VP ∼ 5.86 km and VS ∼ 3.35 km at depths ≳ 5.3 km, which are comparable to the deepest velocities in BACK1 and velocities in model GIL7 at the same depths. Figure 4. View largeDownload slide Depth slices of model G3D1 used to compute SGFs. Dashed line is the outline of the geothermal steam field and black triangles are locations of stations used in this study (same as in Fig. 5). Outside this domain, we extrapolate the velocity model to match the 1-D model BACK1. Figure 4. View largeDownload slide Depth slices of model G3D1 used to compute SGFs. Dashed line is the outline of the geothermal steam field and black triangles are locations of stations used in this study (same as in Fig. 5). Outside this domain, we extrapolate the velocity model to match the 1-D model BACK1. Synthetic Green's functions For all 1-D models, single force SGFs were computed using the frequency–wavenumber integration (FKI) method based on Haskell (1964) and Wang & Herrmann (1980), as provided in Herrmann (2013a). The Herrmann (2013a) FKI code computes complete three-component seismograms consisting of all the terms of the elastic wave equation solution (including the near- and intermediate-field terms) and all body-wave and surface-wave phases for isotropic 1-D layered velocity models with anelastic attenuation. We use 10 m as source depth for all stations at the surface (the numerical integration method used requires a non-zero depth difference between the source and the receiver) and depths reported in Northern California Earthquake Data Center (NCEDC) as source or receiver depths for borehole stations. For computation of all SGFs (both 1-D and 3-D) in this study, we use impulse or Dirac source-time functions for numerical stability. The ground motions thus obtained are displacement waveforms for an impulse input or effectively velocity waveforms for a step input. The velocity waveforms are subsequently integrated to displacement SGFs (for an input step force). We use model G3D1 with the 3-D seismic wave propagation code SW4 to compute SGFs accounting for the 3-D structure in the reservoir. SW4 is a 3-D finite difference seismic wave propagation code that solves the elastic wave equation in displacement formulation with fourth order accuracy in space and time (Sjögreen & Petersson 2012; Petersson & Sjögreen 2014a,b). Absorbing supergrid boundary conditions are used at the far-field boundaries and anelastic attenuation is also allowed (Petersson & Sjögreen 2014b). Our grid spacing is 80 m, which allows ∼27 grid points per shortest wavelength for the minimum VS ∼ 2.0 km s−1 in our model and the maximum frequency of our interest ∼0.9 Hz. For the maximum VP ∼ 6.4 km s−1 in the model, the time step is set to ∼1.18 × 10−2 s for numerical stability. Our model dimensions are ∼36.6 km (NS) x ∼41.2 km (EW) and ∼14.0 km (Z), ∼42 million grid points, with supergrid absorbing layer thickness of 75 grid points. Fig. 5 shows a snapshot of an example wavefield computed using SW4 in response to a single force applied at station DRH. The displacement time histories for this wavefield can be compared to NGFs with DRH as the ‘source’ station. The effect of faster than background velocities in the central section of The Geysers is reflected in the bulge in the otherwise circular ∼0.9 Hz wave front along with reduced amplitudes. Note the absence of any significant reflections from the boundaries of the domain showing the effectiveness of the supergrid absorbing boundary conditions. Figure 5. View largeDownload slide The left subplot shows snapshot of the vertical component velocity wavefield at the surface at elapsed time 9.86 s in response to a vertical force (1e + 10 N) applied at location of station DRH (red star). The amplitudes of the wavefield are normalized by the instantaneous peak absolute amplitude (indicated at the bottom right corner) and plotted on a red (-1) to blue ( + 1) colour scale. The domain shown here corresponds to the computational domain used in SW4 and the dotted line marks the supergrid absorbing layer boundary. Explanation of other symbols on the map is the same as in Fig. 4. The right subplot shows the time history of the same wavefield at station GDXB (red square on the map). The red tick is the instant corresponding to the snapshot shown on the map. A Gaussian source time function with fundamental frequency ∼0.9 Hz (freq parameter in SW4 ∼ 5.65) was used for this particular wavefield. Figure 5. View largeDownload slide The left subplot shows snapshot of the vertical component velocity wavefield at the surface at elapsed time 9.86 s in response to a vertical force (1e + 10 N) applied at location of station DRH (red star). The amplitudes of the wavefield are normalized by the instantaneous peak absolute amplitude (indicated at the bottom right corner) and plotted on a red (-1) to blue ( + 1) colour scale. The domain shown here corresponds to the computational domain used in SW4 and the dotted line marks the supergrid absorbing layer boundary. Explanation of other symbols on the map is the same as in Fig. 4. The right subplot shows the time history of the same wavefield at station GDXB (red square on the map). The red tick is the instant corresponding to the snapshot shown on the map. A Gaussian source time function with fundamental frequency ∼0.9 Hz (freq parameter in SW4 ∼ 5.65) was used for this particular wavefield. Sensor orientations Rotation of the NGF tensor to the T-R-Z reference frame with respect to the source-receiver azimuth requires precise knowledge of sensor orientations. Among the stations employed in our study, the exact orientation of the borehole sensors is unknown (http://www.ncedc.org/egs/, last accessed July 2016). Moreover, preliminary comparisons of NGFs and SGFs (described in the following section) indicated possible sensor misalignments, that is, non-zero angles between the horizontal components of the three-component sensors and the geographic NS–EW reference frame. We compare three-component long period waveforms (∼0.02–0.1 Hz) of regional and teleseismic earthquakes at all stations to the corresponding waveforms at the broad-band USGS station GDXB to verify sensor orientations and determine the degree of misalignment, if any (e.g. Grigoli et al.2012). The details of the methodology and results are presented in the Supporting Information (section Sensor orientations, Figs S1 and S2, and Tables S1 and S2). Our results provide the first orientation estimates for the borehole sensors and show surprisingly large alignment angles (>10°) at many of the surface stations as well (stations ACR, AL6, CLV, DRK, DVB, FUM, HVC, LCK, MNS, NEG and TCH) with standard deviations < 5° over multiple measurements. Procedure for comparison between NGFs and SGFs We find that proper filtering of SGFs is a critical step for comparison with NGFs at small interstation distances. As a result of the band-limiting taper applied during spectral whitening and cross-correlation (see Appendix A1), the NGFs are acausal. For positive correlation lag time, NGFs at short interstation distances will have non-zero acausal amplitudes leaking to t < 0 following zero-phase bandpass filtering that contaminate the symmetric time-reversed NGFs at negative correlation time and vice versa for t > 0 (e.g. Cho et al.2007). Therefore, we zero pad the raw SGFs at t < 0, bandpass-filter them between 0.2 and 1.6 Hz with a zero-phase eight-pole Butterworth filter (same as the taper applied during spectral normalization), time-reverse the SGFs and take the symmetric component. This is explained in detail in Supporting Information Fig. S3. All NGFs and SGFs are subsequently bandpass-filtered between 0.2 Hz and 0.9 Hz with a causal two-pole Butterworth filter. For comparison of the NGF and the SGF tensors, we choose waveform segments between 1 s prior to the theoretical P-wave arrival time or t = 0 (whichever is later) and 5 s after the empirical coda-wave start time. The tensors are normalized by the overall L2 norm of all available components, which preserves relative intercomponent amplitudes and makes the scaling factor invariant under the rotation of the station components. However, this also creates a slight bias against the SGFs of 1-D velocity models, because TR, TZ, RT and ZT components of SGFs are identically zero, whereas corresponding components in NGFs are always non-zero. This is manifested by slightly greater amplitudes for the five primary, non-zero components- TT, RR, RZ, ZR and ZZ of the SGF tensors compared to the same for the NGF tensors for many station pairs, especially at greater distances at which the amplitudes in the other NGF components are dominated by noise in the NGF time-series (see section Contribution of 3-D structure and anisotropy). The overall similarity between the two tensors is evaluated by two estimates of goodness-of-fit, namely by the variance reduction (VR; e.g. Guilhem et al.2014), and by the zero-lag normalized correlation coefficient (CC). P- and S-wave arrival times are computed using velocity model REF for all station pairs excluding stations MNRC and HOPS, and models GIL7 or GIL7.1 for all station pairs including stations MNRC or HOPS, respectively, and are shown on all the following figures that contain waveforms, for reference. We also examine the time delay of each component of SGFs with respect to the corresponding NGF component (e.g. Ma et al.2008) if the normalized CC (not zero-lag) for that component exceeds 0.7. Positive time delay implies that the NGF arrives earlier than the SGF and the actual velocities are faster than the model velocities (at these frequencies, waveforms are primarily sensitive to VS). Synthetic tests Prior to comparing NGF and SGF tensors at The Geysers, we perform tests on ‘synthetic noise’ to evaluate the extent to which phase and relative amplitude information can be retrieved from coherency of seismic noise. Our synthetic tests are modified after synthetic tests done by Herrmann (2013b). We design a 100 km × 100 km domain with stations placed at various positions at the centre allowing us to obtain synthetic noise Green's functions (SNGFs), that is, NGFs estimated from coherency of synthetic noise, for interstation distances from 1 to 35 km (Supporting Information Fig. S4). For constructing synthetic noise records at the stations, we sum filtered (0.2–0.9 Hz) three-component velocity waveforms that are generated by randomly oriented force vectors (amplitude range −1 to +1) at random locations (but at least 50 m away from stations) on the surface with 20 sources acting together every three seconds. Different from the synthetic tests of Herrmann (2013b), we use sources overlapping in time as in Kimman & Trampert (2010) and the FKI method described earlier to compute complete single force responses on the REF model instead of surface wave responses computed by modal summation. The results of the synthetic tests are described in detail in the Supporting Information section Synthetic tests and Supporting Information Fig. S5. Based on the synthetic tests, we find that interstation coherency of velocity noise recorded at two stations compares well with normalized displacement SGFs (for an input step force) for the methodology and scenarios (distances and frequencies, layered attenuating medium) in our study (VR ≳ 80 per cent for distances ∼λ and higher VR at shorter and longer distances) under the assumption of an idealized homogenous noise distribution on the Earth's surface. We also observe better recovery of the Rayleigh wave cross-term components, RZ and ZR, compared to the RR and ZZ components, similar to results of observational (van Wijk et al.2011) and numerical studies (Haney et al.2012). Our failure to retrieve the exact GFs (i.e. VR = 100 per cent) can be attributed to the absence of depth sources, the absence of deformation rate or dipole sources, cross-terms of different phases, incomplete cancellation of overlapping sources and ubiquitous distribution of noise sources instead of their clustering near the stationary phase region (Kimman & Trampert 2010). RESULTS In the context of noise cross-correlations, ‘distances’ in the current and following sections implies interstation distances. In the comparisons of NGFs and SGFs, each SGF tensor corresponds to the velocity model that returns the highest VR with NGFs of that specific station pair, unless stated otherwise. First we discuss NGFs involving stations within the steam field or adjacent to it (i.e. all stations excluding HOPS and MNRC). Effects of sensor misorientations on NGFs In our preliminary analysis, we observed that the NGFs of some stations pairs at very small distances (<λ) had large amplitudes on the TR, TZ, RT and ZT components (hereinafter referred to as T-[R,Z] components) that are unlikely the result of 3-D velocity structure or anisotropy. These anomalous amplitudes were restricted to NGFs involving some specific stations that also returned fairly large alignment angles in the analysis of teleseismic waveforms. Fig. 6(a) shows some examples with surface stations FUM and DVB and borehole station JKB. We consider, sensor misalignments as a likely cause of these non-zero amplitudes. We use alignment angles obtained from earthquake waveforms to rotate the NGF tensors to correct orientations. For the NGFs in Fig. 6(a), rotation to the correct orientation reduces the large amplitudes on the T-[R,Z] components and significantly improves the VR between the SGFs and the NGFs (Fig. 6b). On an average, correcting the NGFs for all alignment angles > 10° increases VR between the SGFs and the NGFs, and reduces the mean root mean square amplitude of the T-[R, Z] components for 76.4 per cent and 79.8 per cent of the 407 interstation pairs, respectively. Figure 6. View largeDownload slide (a) Nine-component NGF (black) and SGF (red) tensors for specific station pairs with large anomalous amplitudes on the T-[R,Z] components. Each panel (3 × 3 group of subplots) is specific to a station pair (names of stations shown at the top of the panels). Each subplot in a panel is specific to a component pair (component pair names are at top right corner of each subplot; distance and azimuth are at the top left corner of the TZ subplots). The first station/component corresponds to the force location/direction and the second station/component corresponds to the receiver location/direction. For each station pair, the velocity model used to compute the best fitting SGF tensor is indicated at the bottom left corner of TZ subplot. For the overall tensor comparisons of a station pair, two goodness-of-fit estimates, VR and normalized zero-lag CC are mentioned at the top left corner of TR subplots. For all non-zero component pairs, the normalized CC along with the time lag between the NGF and SGF waveforms is indicated at the bottom right corners of the subplots, if the CC exceeds 0.7. For reference, the blue and green marks are the theoretical P- and S-wave arrival times, respectively, while the grey mark corresponds to 5 s after the empirical coda start time. For a particular station pair, following normalization of amplitudes of the NGF and SGF tensors, only relative amplitudes are meaningful. (b) Same as (a) but after the NGF tensors are corrected for sensor misalignments. The rotation angles for different stations are indicated on the arrows between (a) and (b). NGFs in all the following figures have been similarly corrected for sensor misorientations using alignment angles obtained from the analysis of earthquake waveforms. Figure 6. View largeDownload slide (a) Nine-component NGF (black) and SGF (red) tensors for specific station pairs with large anomalous amplitudes on the T-[R,Z] components. Each panel (3 × 3 group of subplots) is specific to a station pair (names of stations shown at the top of the panels). Each subplot in a panel is specific to a component pair (component pair names are at top right corner of each subplot; distance and azimuth are at the top left corner of the TZ subplots). The first station/component corresponds to the force location/direction and the second station/component corresponds to the receiver location/direction. For each station pair, the velocity model used to compute the best fitting SGF tensor is indicated at the bottom left corner of TZ subplot. For the overall tensor comparisons of a station pair, two goodness-of-fit estimates, VR and normalized zero-lag CC are mentioned at the top left corner of TR subplots. For all non-zero component pairs, the normalized CC along with the time lag between the NGF and SGF waveforms is indicated at the bottom right corners of the subplots, if the CC exceeds 0.7. For reference, the blue and green marks are the theoretical P- and S-wave arrival times, respectively, while the grey mark corresponds to 5 s after the empirical coda start time. For a particular station pair, following normalization of amplitudes of the NGF and SGF tensors, only relative amplitudes are meaningful. (b) Same as (a) but after the NGF tensors are corrected for sensor misalignments. The rotation angles for different stations are indicated on the arrows between (a) and (b). NGFs in all the following figures have been similarly corrected for sensor misorientations using alignment angles obtained from the analysis of earthquake waveforms. We also consider DEB and SRB that were replacement borehole stations for surface stations DVB and SSR at depths ∼150 m and ∼140 m, and at horizontal distances of ∼30 m and ∼70 m, respectively from the surface sites. Both the surface and borehole stations were simultaneously operational for ∼25 days. The NGFs and the SGFs for the two surface-borehole sensor pairs are shown in Supporting Information Fig. S6. At these small distances, GFs are primarily proportional to filtered source-time functions with the largest near-equal amplitudes on the TT, RR and ZZ components (e.g. Wielandt & Forbriger 1999). Orientation corrections to the NGFs significantly reduce the amplitudes in the off-diagonal components of these tensors (Supporting Information Fig. S6b). Although there can be other factors responsible for the off-diagonal terms, utilizing the NGF tensors may be a useful alternative approach for assessing orientations of borehole sensors. The successful validation of angles measured independently from long period teleseismic earthquake waves on higher frequency NGF tensors provides confidence in our results. Alignment angles at three stations (FUM, NEG and TCH) have also been independently obtained at the sites (Ramsey Haught, private communication, 2016), which are reasonably close to our estimates within the uncertainties (Supporting Information Table S2). In the following sections, all NGFs have been corrected using alignment angles obtained from the analysis of earthquake waveforms. NGF versus SGF at various distances Fig. 7 shows some examples of NGFs and SGFs at distances ≲2 km or ∼0.44 λ. At these small distances, we observe simple waveforms, strong NGFs, low relative amplitudes of noise trailing the ballistic phases with respect to the peak amplitudes, good fits at VR ≳ 40 per cent, and consistency in relative amplitudes of different components of the tensors. Small time delays are sometimes observed between NGFs and SGFs (e.g. AL3-AL5 shown Fig. 7); however, these delays are not consistent across all components, and based on the synthetic tests it is unclear whether they are real or due to bias introduced by the noise distribution or by filtering or processing artefacts. For NGFs at distances >2 km and ≲8 km (or >0.44λ and ≲1.7λ) shown in Fig. 8, we observe a decrease in the quality of retrieved NGFs as evidenced by the increase in relative amplitude of trailing noise. The waveforms also grow slightly more complex and the amplitude misfit between NGFs and SGFs increases, with VR decreasing to ∼20 per cent. We start to see clear inadequacies in the velocity models; for example, in HBW-HER in Fig. 8, the TT component NGF is slower than the SGF indicating that the 3-D model, which is not resolved in this region and therefore is essentially the same as the starting average 1-D model (BACK1) in this area, is too fast for this path. While these distances are too small for reliable phase measurements, large reliable time shifts measured on multiple components of the NGFs with good signal-to-noise ratio can be utilized as ‘broad-band phase delays’ in waveform tomography assuming that multiple scattering in the Earth's crust alleviates the ill-effects of inhomogeneous noise source distribution (Fichtner 2014; Lee et al.2014). Figure 7. View largeDownload slide Same as Fig. 6(b) but for station pairs at interstation distances ≤ 2 km. Figure 7. View largeDownload slide Same as Fig. 6(b) but for station pairs at interstation distances ≤ 2 km. Figure 8. View largeDownload slide Same as Fig. 6(b) but for station pairs at interstation distances > 2 km and ≤ 8 km. Figure 8. View largeDownload slide Same as Fig. 6(b) but for station pairs at interstation distances > 2 km and ≤ 8 km. Fig. 9 shows examples for distances >8 km. These are characterized by lower signal-to-noise ratio, larger amplitudes of trailing noise and large amplitudes on the TR, TZ, RT and ZT components. The simple waveforms on the TT components consist of Love waves that depend on VS structure only and show good fits, whereas the waveforms on the RR, RZ, ZR and ZZ components consisting of P–SV and Rayleigh waves depend on both, VP and VS structure and exhibit rapidly increasing misfits with distance. However, the RR components of the NGFs in Fig. 9 compare well to the SGFs. P–SV body-wave energy arriving before the theoretical S-wave arrival time can be clearly identified above the noise floor in the ZR components of AL6-NEG, AL2-STY, AL3-DES and AL2-DVB. Similar to results of Lin et al. (2008), but for a different distance and frequency range, we obtain higher amplitude coherency in the TT component than in the RR and the ZZ components at distances >2λ that compare well to the SGFs. Therefore, Love waves on the TT component NGFs should be used to constrain the shear wave velocity models in addition to Rayleigh waves whenever three-component sensors are available, as the former are generally retrieved with better signal-to-noise ratio. Similar to our synthetic tests and observations in other studies (van Wijk et al.2011), we generally find good fits between ZR and RZ component NGFs and SGFs at all distances. Figure 9. View largeDownload slide Same as Fig. 6(b) but for station pairs at interstation distances > 8 km. Figure 9. View largeDownload slide Same as Fig. 6(b) but for station pairs at interstation distances > 8 km. Compiling overall results, the average peak-to-peak amplitude ratio of the primary components (TT, RR, RZ, ZR, ZZ) of NGFs to the corresponding components of SGFs for all station pairs is ∼0.7. The standard deviation is ∼0.18 in log10 units (i.e. difference of a factor of ∼1.5, bigger or smaller) with respect to the mean amplitude ratio. To select only robust time delays between individual components of NGFs and corresponding components of the best-fitting SGFs for interpretation, we require that VR for a component pair is ≥ 35 per cent following correction for the time delay estimated from normalized cross-correlation. This criterion ensures similarity in both waveform shape and amplitude. Out of 238 station pairs that have at least 2 such components, only ∼19 per cent of station pairs have large contrasting time delays on different components (maximum time delay ≥0.1 s and minimum time delay ≤ –0.1 s) indicating possible errors in the recovered NGFs. For many station pairs such as HBW-HER in Fig. 8, large robust time delays on different components of the tensors are of the same sign, indicating an overall consistency in observations. Contribution of body waves and surface waves To study the contributions of different body-wave and surface-wave phases to the observed NGF waveforms, we compare pure surface-wave SGFs computed using the modal summation method as provided in Herrmann (2013a) to complete FKI SGFs. We use the REF velocity model for this analysis; the synthetic phase velocities and depth-dependent eigenfunctions are computed using codes provided in Herrmann (2013a). Phases that are present in the complete FKI synthetics but absent in the surface-wave synthetics are interpreted to be body-wave phases. Synthetic phase velocity dispersion curves based on the REF model for fundamental- and first higher-mode Rayleigh waves are plotted in Fig. 10(a). Fig. 10(b) shows an enlarged view of the fast arrivals at ∼3.6 km s−1 in the ZZ component NGFs at distances 17–20 km in Fig. 2. Fig. 10(c) shows SGFs at similar distances containing the fundamental-mode (‘Mode 0’) and first higher-mode Rayleigh waves (‘Mode 1’) and their combination (‘Mode 0+1’). For a source at ∼10 m depth, the peak amplitudes of the higher-mode are ∼14 times smaller than peak amplitudes of the fundamental-mode and therefore, the higher mode has a very small contribution to the total waveforms containing both modes. The amplitudes of the second higher-mode Rayleigh wave are even smaller and therefore, it is not discussed here. The recovery of higher-mode surface waves on land is more difficult than the recovery of fundamental-mode surface waves primarily due to the near-surface concentration of ambient noise sources (Kimman & Trampert 2010). However, some studies have reported success in this regard (e.g. Nishida et al.2008; Lin et al.2013). The body-wave energy observed in the FKI synthetics is considerably weaker than energy observed in the ZZ component NGFs between the theoretical P-wave and S-wave arrival times. Compared to the SGFs, the NGFs appear to be depleted in high frequencies (≳0.5 Hz) as well. More detailed analyses including the analysis of particle motions and SGFs computed with other velocity models are required to positively identify this phase in the ZZ component. However, since it has greater amplitudes and is faster (∼3.6 km s−1) than the first higher-mode Rayleigh wave expected for the REF model (phase velocity ∼3.0–3.3 km s−1), it is likely to be body-wave energy. Figure 10. View largeDownload slide (a) Synthetic phase velocity dispersion curves for the fundamental-mode and first higher-order mode Rayleigh waves for the REF velocity model. (b) Record section of the ZZ component NGFs at interstation distances of 17–20 km. For reference, the blue and green solid lines denote the theoretical P- and S-wave arrival times, respectively. The two dashed red lines mark the arrival times for apparent medium velocities of 3.6 and 2.5 km s−1. (c) REF model ZZ component SGFs at interstation distances of 17–20 km for different Rayleigh waves modes (‘Mode 0’ and ‘Mode 1’) computed using the modal summation technique, compared to complete synthetics computed using FKI. The type of modes/method is indicated at the top right corner of each subplot. ‘Mode 0’ and ‘Mode 1’ indicate fundamental and first higher-mode Rayleigh waves, respectively. ‘Mode 0+1’ synthetics include both ‘Mode 0’ and ‘Mode 1’. The range of maximum absolute amplitudes of traces in each subplot is indicated in the bottom left corner. For example, the maximum absolute amplitude of ZZ component SGF computed using FKI varies from 4.7 × 10–4 to 5.4 × 10–4 cm over distances 17–20 km. Figure 10. View largeDownload slide (a) Synthetic phase velocity dispersion curves for the fundamental-mode and first higher-order mode Rayleigh waves for the REF velocity model. (b) Record section of the ZZ component NGFs at interstation distances of 17–20 km. For reference, the blue and green solid lines denote the theoretical P- and S-wave arrival times, respectively. The two dashed red lines mark the arrival times for apparent medium velocities of 3.6 and 2.5 km s−1. (c) REF model ZZ component SGFs at interstation distances of 17–20 km for different Rayleigh waves modes (‘Mode 0’ and ‘Mode 1’) computed using the modal summation technique, compared to complete synthetics computed using FKI. The type of modes/method is indicated at the top right corner of each subplot. ‘Mode 0’ and ‘Mode 1’ indicate fundamental and first higher-mode Rayleigh waves, respectively. ‘Mode 0+1’ synthetics include both ‘Mode 0’ and ‘Mode 1’. The range of maximum absolute amplitudes of traces in each subplot is indicated in the bottom left corner. For example, the maximum absolute amplitude of ZZ component SGF computed using FKI varies from 4.7 × 10–4 to 5.4 × 10–4 cm over distances 17–20 km. Comparing surface-wave (‘Mode 0/1/0+1’) SGFs and complete FKI SGFs at distances of 14–17 km for the ZR and RR components (Supporting Information Fig. S7), we positively identify strong P-wave energy in the RR and ZR component NGFs for many station pairs (e.g. AL3-DES in Fig. 9) in light of the negligible contribution of the first higher-mode Rayleigh wave compared to the fundamental mode. NGFs with USGS 1.0 Hz sensors Fig. 11 shows GF comparisons for station pairs with one station being a USGS vertical component short period (∼1.0 Hz) station, for which only the TZ, RZ and ZZ components are available (with the vertical component station being the receiver). The NGF-SGF comparisons for these station pairs are important as four of these stations (GAXB, GBG, GGPB and GSG) lie well outside the reservoir area (>3 km away from the reservoir outline and the LBNL network) and none of our 1-D velocity models are calibrated to fit these paths. Moreover, a single background 1-D model is assumed in model G3D1 for almost all regions outside the reservoir area and an artificial, unrealistic boundary exists between the 1-D and the 3-D models along the reservoir outline, except to the northwest. Model REF, which is a 1-D average of a 3-D model of NW Geysers, fits most paths to stations GPM and GSG that lie to the NW and NE of The Geysers, respectively. However, large time delays, −0.4 to −0.75 s, are clearly observed at station GSG (see waveforms for ACR-GSG, JKR-GSG, STY-GSG, FUM-GSG in Fig. 11) consistent with body-wave tomography studies that have determined lower velocities to the northeast of The Geysers attributed to rocks of the Clear Lake volcanics and interbedded sedimentary deposits of the Great Valley sequence (Hearn et al.1981; Eberhart-Phillips 1986; Julian et al.1996). The delays are also consistent with a 0.25 s P-wave station correction at GSG with respect to GBG (Eberhart-Phillips & Oppenheimer 1984). The retrieval of robust NGFs from these sensors demonstrates their potential in supplementing temporary deployments of higher quality broad-band sensors for ambient noise Rayleigh wave tomography studies (e.g. Porritt et al.2011). This is especially true for regions like Northern California, where they are a significant part of pre-existing seismic infrastructure in the Northern California Seismic Network (http://www.ncedc.org/ncsn/, last accessed July 2016). Figure 11. View largeDownload slide Same as Fig. 6(b), but for all distances and stations pairs with one USGS vertical component 1.0 Hz sensor. Interstation distance and azimuth, and the overall goodness-of-fit estimates are at the top of TZ subplots. Figure 11. View largeDownload slide Same as Fig. 6(b), but for all distances and stations pairs with one USGS vertical component 1.0 Hz sensor. Interstation distance and azimuth, and the overall goodness-of-fit estimates are at the top of TZ subplots. We were also able to retrieve robust NGFs with the accelerometer at station DRH (Supporting Information Fig. S8) that compare very well to SGFs at distances up to 20 km. The RR component NGFs for station pairs CLV-DRH, HVC-DRH and ACR-DRH also show strong P-waves (see Supporting Information Fig. S7a for comparison). Low gain accelerometers have been successfully used in some ambient noise tomography studies (e.g. Cho et al.2007). CONTRIBUTION OF 3-D STRUCTURE AND ANISOTROPY Even after the orientation corrections, the NGFs for some longer distance paths (≳8 km or ∼1.7 λ) show significant non-zero amplitudes on the T-[R,Z] components that appear robust with respect to the trailing noise (Fig. 12). These waveforms are reproduced well by the SGFs computed using the 3-D velocity model indicating the effects of significant 3-D structure at distances ≳2 λ. We visually identified these interstation paths, which appear to be approximately parallel to the southwestern boundary of the reservoir area along the extent of the high-velocity region in the southeast section of the reservoir. While most of these paths involve stations SSR and SRB, the teleseismic wave analysis described in the Supporting Information indicates that the sensor installed at SSR was correctly oriented. It is possible that off-great-circle wave propagation or multipathing effects along the large velocity contrast cause these large amplitudes on the T-[R,Z] components. However, the 3-D velocity model is poorly resolved outside the reservoir boundary (Gritto et al.2013a). There are other interstation paths for which we see robust non-zero amplitudes on these components but they are not reproduced by the 3-D velocity model. Figure 12. View largeDownload slide Same as Fig. 6(b), but for some station pairs with significant non-zero amplitudes on the T-[R,Z] components of NGFs even after corrections for sensor misalignments. These amplitudes are reproduced well by the 3-D velocity model. Figure 12. View largeDownload slide Same as Fig. 6(b), but for some station pairs with significant non-zero amplitudes on the T-[R,Z] components of NGFs even after corrections for sensor misalignments. These amplitudes are reproduced well by the 3-D velocity model. Assuming far-field noise source incidence only, if the station pairs are not aligned with the dominant noise source incidence direction, non-zero amplitudes can be expected in T-[R,Z] components of NGFs (Roux 2009; Durand et al. 2011). It is possible to minimize the amplitudes in these components by rotating the reference frames at the two stations so that the new reference frames are approximately aligned with the azimuth of dominant noise source incidence direction [the method is referred to as Optimal Rotation Algorithm (ORA); Roux 2009]. We quantify the relative amplitudes of T-[R,Z] components of the NGFs at The Geysers, in terms of the misfit parameter φ, defined as the ratio of sum of squares of amplitudes of the four off-diagonal components to the sum of squares of all components of the NGF tensor (eq. 17 in Saade et al. 2015). We examined φ as a function of azimuth of the station pairs to investigate any systematic relationship. Our results, which are fully described in the Supporting Information section Application of Optimal Rotation Algorithm, indicate that the non-zero signals in T-[R,Z] components of our NGFs at large distances are just leading or trailing noise in the NGF time-series rather than coherent signals resulting from systematic azimuthal deviation of interstation paths from a particular noise source direction. It is possible that seismic anisotropy also contributes to the non-zero amplitudes in these components at larger distances by giving rise to quasi-Love or quasi-Rayleigh waves (Maupin & Park 2007). The shallow subsurface (top ∼3–5 km) at The Geysers can be considered as a horizontal transversely isotropic (HTI) medium with an average ∼4 per cent VS anisotropy (up to ∼11 per cent at some locations) estimated from shear-wave splitting measurements on local earthquakes records (Majer et al.1988; Elkibbi & Rial 2005; Elkibbi et al.2005). The VS anisotropy is defined as $$\frac{{{V_{S,{\rm{fast}}}} - {V_{S,{\rm{slow}}}}}}{{{V_{S,{\rm{fast}}}}}}$$ × 100 as in Elkibbi et al. (2005). The anisotropy at The Geysers is believed to be caused by stress-induced alignment of fractures and cracks. While the fast axes are generally parallel or sub-parallel to the N-to-NE direction of the regional maximum compressive stress, a smaller set of NW–SE oriented fast axes are also observed, primarily in the southeast Geysers and likely related to local fault shearing effects (Elkibbi & Rial 2005; Elkibbi et al.2005). The 3-D model of Gritto et al. (2013a,b) was inverted utilizing a well distributed set of earthquakes and stations, and therefore can be considered as an isotropic average of the velocity structure at The Geysers, with a maximum of 19 per cent ± 7 per cent variation in VS at depths between 0.5 km and 4.0 km. Therefore, seismic wave propagation at The Geysers at large distances is possibly controlled, to first order, by the average 3-D isotropic structure in the reservoir. In order to test the amplitudes caused by anisotropy and expected in the T-[R,Z] components of the NGFs, we also compute SGFs for station pairs at The Geysers incorporating an HTI medium in the 1-D velocity models. The 3-D finite-difference seismic wave propagation code SW4 provides the capability to model seismic waveforms for layered anisotropic but purely elastic media. Therefore, we restrict this analysis to stations pairs at distances <10 km, so that we can ignore the effects of anelastic attenuation (further discussed in section Amplitude decay). We also ignore the considerable scatter in the degree of anisotropy and in the fast axis directions across The Geysers that were observed by Elkibbi & Rial (2005) and Elkibbi et al. (2005). We use the formulation of Hudson (1981, 1982) as provided in eqs (1) to (4) in Crampin (1984) to estimate the five anisotropic elastic constants for an HTI medium—A, C, F, L and N in Love notation (Saade et al.2015). The differences between (A, C) and (L, N) determine the VP and VS anisotropy, respectively. A review on the $$\eta \ = \frac{F}{{A - 2L}}\ $$ parameter can be found in Kawakatsu et al. (2015). The effective anisotropic elastic tensor in a cracked medium can be modelled as the sum of the elastic tensor for the uncracked solid and first and second order perturbation terms that depend on crack density (ε), crack aspect ratio (d), Lame's constants for the uncracked solid, and Lame's constants for the weak crack inclusions (λ', μ'). This formulation assumes a weak distribution of parallel disconnected penny-shaped cracks with dimensions and spacing much smaller than the seismic wavelengths under consideration and is valid for small values of ε < 0.1 (Crampin 1984; Peacock & Hudson 1990; Cheng 1993). The Geysers is characterized by a vapour-dominated geothermal field with a considerable volume of water injected into the reservoir. For this study, we assume that the cracks are filled with water ( VP = 1.5 km s− 1, λ' = 2.25e9 Nm,  μ' = 0). It is important to note that the presence of steam in water can drastically reduce VP of the two-phase mixture, especially if the water and steam phases are in thermodynamic equilibrium. However, this effect is expected to decrease with increasing pressures and temperatures (Liu & Kieffer 2009). Majer et al. (1988) investigated one site at The Geysers and found little or no evidence for VP anisotropy. Therefore, we assume thin cracks, d ∼ 0, which gives A = C. An assumed value of ε ∼ 0.036 leads to an approximate VS anisotropy ∼4 per cent in the layers of the top ∼3.1 km of the 1-D velocity models. We obtain values of η ∼ 0.82–0.89, compared to η = 1.1 as adopted by Saade et al. (2015) based on scaling relationships (Montagner & Anderson 1989). We test two directions of crack normals, N60° W and N125° W, transverse to the average directions of the fast axes determined by Elkibbi & Rial (2005). For the computation in SW4, the crack normal is assumed to be in direction 1 and the station paths are rotated accordingly. Fig. 13(a) shows NGF and SGF comparisons for 2 station pairs at distance < 10 km for which the significant non-zero amplitudes on the T-[R,Z] components of the NGFs are reproduced well by the SGFs computed using the 3-D velocity model. Fig. 13(b) shows the same NGFs as Fig. 13(a) but with SGFs computed with an elastic (no anelastic attenuation) and anisotropic version of the 1-D velocity model VSP0, incorporating ∼4 per cent VS anisotropy in the top ∼3.1 km layers with crack normal direction N125° W for station pairs FUM-SSR and DRK-SSR. The VR and the amplitudes on the T-[R,Z] components for the 1-D anisotropic models are clearly smaller than those for the 3-D model. The 1-D anisotropic model SGFs yield higher VR compared to the isotropic 1-D or 3-D models for only 13 out of ∼270 nine-component station pairs (∼5 per cent) at distances < 10 km; for 6 out of 13 station pairs, the improvement in VR was negligible (<1 per cent). In the case of the SGFs dominated by surface waves in an anisotropic medium, ORA can also be used to minimize the amplitudes in the T-[R,Z] components (Durand et al.2011; Saade et al.2015). We also perform a synthetic test, in which we compute SGFs assuming the anisotropic REF model with one reference station at the centre of a circle and other stations placed at radial-distance intervals of 0.5 km up to a total distance of 10 km and 2.0° azimuthal spacing. For these GFs, the maximum value of φ ranges from ∼4e-4 to ∼1e-2 for distances between 0.5 to 10 km. These values are smaller than φ for our NGFs at The Geysers (Supporting Information Fig. S9a) at distances <10 km. This observation, combined with the failure of ORA to significantly reduce φ in our NGFs (Supporting Information Fig. S9b), indicates that the large residual amplitudes in the T-[R,Z] components of the NGFs for most station pairs (<10 km) at The Geysers correspond to just leading or trailing noise in the NGF time-series. These non-zero signals are unlikely to be real coherent signals resulting from systematic effects of non-uniform illumination from far-field noise sources or currently known subsurface anisotropy in the geothermal field. For at least some station pairs that have robust signals in the T-[R,Z] components, the 3-D model provides better fits than the anisotropic 1-D models. We leave analysis of the effects of anisotropy at distances > 10 km for future studies. Figure 13. View largeDownload slide (a) Same as Fig. 12 but for two different station pairs at interstation distances <10 km. (b) Same as (a) but for SGFs computed using an elastic and anisotropic version of the 1-D velocity model VSP0 with ∼4 per cent VS anisotropy to a depth of ∼3.1 km and crack normal in the direction N125° W. ‘–A’ has been added to the model names to indicate anisotropy. SGFs computed using the 3-D velocity model fit the significant non-zero amplitudes on the T-[R,Z] components of NGFs better than SGFs computed using the 1-D anisotropic velocity model. Figure 13. View largeDownload slide (a) Same as Fig. 12 but for two different station pairs at interstation distances <10 km. (b) Same as (a) but for SGFs computed using an elastic and anisotropic version of the 1-D velocity model VSP0 with ∼4 per cent VS anisotropy to a depth of ∼3.1 km and crack normal in the direction N125° W. ‘–A’ has been added to the model names to indicate anisotropy. SGFs computed using the 3-D velocity model fit the significant non-zero amplitudes on the T-[R,Z] components of NGFs better than SGFs computed using the 1-D anisotropic velocity model. EVALUATION OF VELOCITY MODELS Fig. 14 shows a summary of ∼385 interstation noise coherency paths with the best-fitting VR between the NGFs and the SGFs ≥ –20 per cent, grouped by the best-fitting velocity models, that is, the models that provide the SGFs with the highest VR for the NGFs along these paths. The VR threshold of −20 per cent is a quality check on both, the quality of the NGFs (poor quality NGFs will return poor VR upon comparison with SGFs), and the ability of velocity models to fit the NGFs (paths along which the velocity models are found to be deficient [e.g. paths to station GSG] are filtered out). Out of the four 1-D models tested, REF provides the highest VR for the most number of paths, primarily across NW and central Geysers (Fig. 14a). This can be attributed to its origin from a 3-D model of the NW Geysers region. Similarly, VSP0, which corresponds to a velocity profile in the SE Geysers, is the best-fitting velocity model for most NGFs across the SW- and the SE Geysers (Fig. 14c). The 3-D model G3D1 provides the best-fitting SGFs for a slightly higher number of NGF paths compared to REF (Fig. 14e). To find regions where the 3-D model fits the NGF waveforms significantly better than the 1-D velocity models in a relative sense, we examine paths for which the VR of the 3-D model SGFs exceeds those of the 1-D models by an arbitrarily chosen threshold of 7 per cent or more (Fig. 14f). Most of these paths are across the transition zone between the slower NW section and the faster SE section of The Geysers. The 3-D model also incorporates a prominent low VP/VS ratio anomaly in this area with estimates (down to ∼1.42) considerably lower than those in the 1-D models (Figs 3 and 4). Therefore, it is not surprising that the 3-D model fits long distance paths across this complex region better than the 1-D models. Models AVG1 and BACK1 provide minor improvements over other velocity models for a smaller number of paths in the NW Geysers (Figs 14b and d). As expected, the 1-D models fit most of the paths to the stations outside the reservoir area (stations GAXB, GBG, GGPB, GSG) better than G3D1 (59 out of 62). Figure 14. View largeDownload slide First five subplots (a–e) show interstation paths grouped by best-fitting velocity models. The velocity model and the number of paths are mentioned in the top and bottom left corner, respectively. Subplot (f) ‘3D BEST’ shows paths for which VR of NGFs with the 3-D model SGFs exceed the VR with the 1-D model SGFs by 7 per cent or more (34 out of 124). For all the paths shown in this figure, the best-fitting VR is ≥–20 per cent. Figure 14. View largeDownload slide First five subplots (a–e) show interstation paths grouped by best-fitting velocity models. The velocity model and the number of paths are mentioned in the top and bottom left corner, respectively. Subplot (f) ‘3D BEST’ shows paths for which VR of NGFs with the 3-D model SGFs exceed the VR with the 1-D model SGFs by 7 per cent or more (34 out of 124). For all the paths shown in this figure, the best-fitting VR is ≥–20 per cent. Next we examine phase differences represented by robust time delays (defined in section NGF versus SGF at various distances) between NGFs and SGFs computed with different velocity models for some station pairs. We first consider NGFs and SGFs for intermediate distance (between 2 km and 12 km) interstation paths in ∼NW and Central Geysers (Figs 14a and b), for which REF and BACK1 are the best-fitting velocity models. Barring some large outliers (<–1.0 s or >1 s; ≲2 per cent of measurements removed), the robust time delays for these pairs of NGFs and SGFs are ∼–0.02 ± 0.17 s (mean ± 1 standard deviation; for ∼240 measurements). For NGFs along the same paths, SGFs computed using the faster velocity model VSP0 lead to time delays of ∼–0.12 ± 0.18 s, indicating slightly earlier arrival of SGFs with respect to NGFs. Similarly, time delays between NGFs between SGFs computed with the best-fitting model VSP0 for intermediate distance paths in ∼SE Geysers (Fig. 14c) are ∼+0.12 ± 0.11 s (for ∼75 measurements). The actual velocities could be slightly greater than VSP0 velocities. For NGFs along the same paths, SGFs computed using the slower velocity model REF lead to time delays of ∼+0.26 ± 0.16 s, indicating more delayed arrival of SGFs with respect to NGFs. Our interstation paths in Figs 14(a)–(c) are well distributed in azimuth which is expected to help reduce the effects of phase errors in NGFs due to non-uniform noise source distribution. The meaningful grouping of interstation paths in various sub-regions of The Geysers (Figs 14a–c) indicates an overall consistency of the results. NGFs WITH BDSN STATIONS HOPS AND MNRC We first analyse the simple waveforms of Love waves on the TT components of the NGFs from BDSN stations HOPS and MNRC to stations at the periphery of the reservoir area to gain insight into the shear wave velocity structure in the region to the northwest and to the east of The Geysers (Fig. 15a). The REF model SGFs are too fast compared to the NGFs (Figs 15b and c), which is consistent with body-wave traveltime studies that suggest that the reservoir area has faster than regional seismic velocities (Majer & McEvilly 1979; Eberhart-Phillips 1986). The SGFs of the regional velocity model, GIL7 (Fig. 3) can reasonably fit the Love waves on the paths from the eastern boundary of the reservoir to MNRC (Fig. 15b). GIL7 SGFs are also consistent with the NGFs of paths from some of these stations to station GSG (Fig. 15e), for which the REF model was determined to be too fast (Fig. 11). But we find that GIL7 is too slow for paths from the northwestern boundary of the reservoir area to HOPS (not shown here) with broad-band phase delays from ∼1.3 s to ∼2.1 s. Therefore, we iteratively perturb the VS values in the top ∼5 km of GIL7 to improve the waveform fits of the TT components of the NGFs for these paths. Primarily deceasing the VS values in the top ∼4 km layers of GIL7 by ∼0.1–0.3 km s−1 (model GIL7.1 shown in Fig. 3) produces satisfactory waveform fits (Fig. 15c). In addition to Love waves, clear SH waves can also be identified in the NGFs of some paths to station HOPS (from HBW, RGP, MCL, BRP) that compare well to the SGFs. Figure 15. View largeDownload slide (a) Map showing interstation paths (red and green lines) from stations on the periphery of the reservoir area (white polygon) to BDSN stations HOPS and MNRC (red triangles), respectively. Meaning of other symbols is same as in Fig. 1; (b) TT component NGFs (black traces) and SGFs (red traces) computed using REF (left subplot) and GIL7 (right subplot) velocity models for paths from stations at the eastern boundary of The Geysers to station MNRC arranged in the order of increasing interstation distance. Station names and azimuths are shown near the traces; (c) same as (b) but for paths from stations at the northwestern boundary of The Geysers to HOPS. The velocity models are REF (left subplot) and GIL7.1 (right subplot); (d) Observed Love wave phase velocity dispersion curves measured on TT component NGFs for paths in Figs 15(a)–15(c) plotted against phase velocities predicted for models REF, GIL7.1 and GIL7; (e) same as Fig. 11 but for station pairs including GSG and stations at the northeastern boundary of the reservoir and using velocity model GIL7; (f) same as Fig. 9 but for a station pair including MNRC. Figure 15. View largeDownload slide (a) Map showing interstation paths (red and green lines) from stations on the periphery of the reservoir area (white polygon) to BDSN stations HOPS and MNRC (red triangles), respectively. Meaning of other symbols is same as in Fig. 1; (b) TT component NGFs (black traces) and SGFs (red traces) computed using REF (left subplot) and GIL7 (right subplot) velocity models for paths from stations at the eastern boundary of The Geysers to station MNRC arranged in the order of increasing interstation distance. Station names and azimuths are shown near the traces; (c) same as (b) but for paths from stations at the northwestern boundary of The Geysers to HOPS. The velocity models are REF (left subplot) and GIL7.1 (right subplot); (d) Observed Love wave phase velocity dispersion curves measured on TT component NGFs for paths in Figs 15(a)–15(c) plotted against phase velocities predicted for models REF, GIL7.1 and GIL7; (e) same as Fig. 11 but for station pairs including GSG and stations at the northeastern boundary of the reservoir and using velocity model GIL7; (f) same as Fig. 9 but for a station pair including MNRC. We also compare the Love wave phase velocities extracted from the TT component in the NGFs for these paths (in Figs 15b and c) with synthetic phase velocities for the 1-D velocity models (Fig. 15d). We first apply multiple filter analysis (Dziewonski et al.1969; Herrmann 1973, 2013b) on the NGFs to measure group velocities at periods ∼1.1 to ∼5.0 s (or ∼0.2 to ∼0.9 Hz) and then determine phase velocities using the reference dispersion curves to resolve the 2πN ambiguity, in which N is an integer (Bensen et al.2007; Lin et al.2008). In this frequency passband, observed phase velocities for paths to station HOPS are clearly greater than those for station MNRC by ∼0.1 to ∼0.3 km s−1, and they agree well with phase velocities predicted by GIL7.1 and GIL7 models, respectively. All of these phase velocities are slower than the REF model velocities by ∼0.3 to ∼0.7 km s−1. While it is possible that other velocity models (different from GIL7 and GIL7.1) might provide equal or better fits to these NGFs, it appears that to the first order, shallow crustal shear wave velocities in the region to the northwest of The Geysers are higher than those to the east of The Geysers, and they are both lower than velocities within the reservoir area. Comparing the RR, ZZ, ZR and RZ components suggests that the NGFs at most stations at The Geysers with BDSN stations HOPS and MNRC at distances >30 km are characterized by complicated waveforms, poor fits with synthetic waveforms of GIL7 and its modified versions, and higher relative amplitudes of trailing noise. However, for some stations pairs (e.g. FUM-MNRC), good quality NGFs are retrieved that show remarkable similarity to the SGFs computed with GIL7 with slightly modified values (VS in top 3 layers increased by 100, 50 and 20 m s−1, respectively) in the shallow crust (Fig. 15f). Note the large amplitudes in the TR and RT component of the NGFs in Fig. 15(f), in which FUM has already been corrected for sensor misalignment (Fig. 6). AMPLITUDE DECAY We study the decay of ground motion amplitudes in the reservoir area by analysing the NGF amplitudes as a function of frequency f and interstation distance r between 1 and 25 km. Observations and interpretation of TT component amplitudes First, we focus on the TT component as they consist of relatively simple SH and Love waves, and usually have the largest relative amplitudes among all components of NGFs (see Figs 6–9, 12 and 13, and Supporting Information Figs S8 and S9c). We use the NGF waveform segments in the data window prior to the causal filtering and the tensor amplitude normalization steps described in subsection Procedure for comparison between NGFs and SGFsin the Methodology section. We extract spectral amplitudes from the smoothed Fourier amplitude spectra of the waveform segments at multiple frequencies, f = 0.25 Hz, 0.42 Hz and 0.72 Hz. REF model Love wave group velocities UL(f) at these frequencies are ∼2.44, 2.29 and 2.22 km s−1, respectively. The Fourier Spectral Amplitudes (FSAs) are extracted by interpolating the spectra at 7 points around the frequency of interest with half-width ∼0.1 Hz and averaging the amplitudes with a 7-point Hanning window. We compare the NGF FSAs with the FSAs of the equivalent SGFs of the 3-D and the 1-D REF models that are computed the exact same way. A second group of FSAs of raw synthetic Green's functions (RSGFs) is extracted from the raw unfiltered velocity GFs of the same duration prior to the integration, zero-phase bandpass filtering, time-reversal and symmetric component steps in subsection Procedure for comparison between NGFs and SGFs (Supporting Information Fig. S3). This allows us to evaluate the degree of contamination of the SGF amplitudes with respect to the ‘true’ amplitudes and compare the decay characteristics better to earthquake ground-motion amplitudes. Without the contamination, the distance scaling of RSGF spectral amplitudes should be approximately similar to that of NGF amplitudes in a narrow frequency passband. To distinguish the effects of distance-dependent geometrical spreading, G(r), and total anelastic attenuation, represented by the quality factor Q(f), we also compare the FSAs of SGFs computed for velocity models with weak and strong anelastic attenuation. SW4 allows computation of synthetic seismograms for purely elastic 3-D velocity models. Synthetic seismograms for the weak anelastic attenuation 1-D models were computed by setting QP and QS to high values (1000 and 500, respectively) in all layers. For the strong anelastic attenuation 1-D models, QP and QSare set to ∼15 for depths < 5.8 km and ∼25 for greater depths. In this paper, we primarily focus on FSAs at 0.72 Hz, as the effects of anelastic attenuation are expected to be more prominent at higher frequencies. Fig. 16(a) shows the TT component of the RSGF FSAs at 0.25 Hz as a function of r. The scatter in the amplitudes for the 3-D velocity model is caused by smooth 3-D variations in the seismic velocities. At 0.25 Hz and for r < 25 km, the ground motion decay is very similar for both, low Q and high Q models, implying that the effect of anelastic attenuation can be neglected at these low frequencies and small distances. Directly approximating G(r) by the amplitude decay, we obtain G(r)  ∼  r− 0.82 to G(r)  ∼  r− 1.14, implying dominance of body waves. This decay is stronger than the surface wave geometrical spreading(r− 0.5) expected at distances ≳1λ–2λ (1λ ∼ 11.4 km for Love wave phase velocity cL∼2.85 km s−1 at ∼0.25 Hz) for surface sources. We also plot the RSGF FSAs of the 3-D elastic model and 1-D weak attenuation AVG1 model at 0.72 Hz. AVG1 is an average 1-D model of the 3-D velocity model within the reservoir area. Using various approximations, Menon et al. (2014) have shown that azimuthally averaging the coherence estimates of different station pairs for a 3-D inhomogeneous but elastic velocity structure can introduce apparent anelastic attenuation in the decay of average coherence amplitudes if one were simply fitting an exponential decay model to the amplitudes. The similarity in decay characteristics of the FSAs obtained from the average 1-D and 3-D weak attenuation models at both 0.25 and 0.72 Hz suggests that it might be possible to detect amplitude decay caused by anelastic attenuation following proper averaging of a 3-D velocity structure across all paths (including the ones outside the reservoir), at least for the wavelengths and heterogeneity scales analysed in this study. Figure 16. View largeDownload slide (a) TT component RSGF FSAs for different velocity models and frequencies on log10 scale as a function of distance. The decay curves are scaled so that they have similar amplitudes at ∼1–2 km. 1-D and 3-D indicate FSAs of synthetics computed with 1-D and 3-D velocity models, respectively. ‘No Q’ and ‘Q∼15’ indicate synthetics computed with weak and strong attenuation models, respectively. 1-D RSGF decay curve at 0.72 Hz used the AVG1 model unlike other 1-D model decay curves that use the REF model. (b) TT component NGF and SGF FSAs at 0.25 Hz as a function of interstation distance. The NGF FSAs (grey +) are uniformly scaled up by a constant such that the maximum value is 0. The frequency is indicated at the top left corner of the plots. Black diamonds and error bars represent the mean of NGF amplitudes and their standard deviation ± 1 σ in bins for all bins with more than five data points. The bin widths are 1 and 2 km for distances ≤10 km and >10 km, respectively. ‘SGF 3D’ and ‘SGF 1-D’ are FSAs of equivalent SGFs of the 3-D and 1-D REF velocity models, respectively, plotted for the same station pairs. The FSAs of SGFs are uniformly scaled by a constant to minimize the L1 norm of their difference with the binned NGF FSAs. (c) Similar to (b) but for FSAs at 0.72 Hz. We also plot synthetic FSAs of weak and strong attenuation models. The SGF decay curves are scaled such that their amplitude at ∼1.5 km is similar to the average amplitude of NGF FSAs in the 1–2 km range. Scaling the binned NGF decay curve instead, to directly fit the REF model SGF decay curve by minimizing the L1 norm (black dashed line), doesn't change the absolute amplitudes significantly. The RSGF decay curve (red dashed line) is scaled to fit the REF model SGF decay curve at distances ≳ 6 km. (d) The data (green +) are tangential component FSAs of earthquake records at 0.72 Hz and the red diamonds and error bars represent the mean amplitudes and their standard deviation ± 1 σ in the bins. For comparison, we also plot binned TT component NGF FSAs at 0.72 Hz (black diamonds and error bars; same as [c]) that are scaled to minimize L1 norm of binned FSAs at distances ≳ 3 km (∼1 wavelength). ‘N’ at the top right corner indicates the number of data points. Figure 16. View largeDownload slide (a) TT component RSGF FSAs for different velocity models and frequencies on log10 scale as a function of distance. The decay curves are scaled so that they have similar amplitudes at ∼1–2 km. 1-D and 3-D indicate FSAs of synthetics computed with 1-D and 3-D velocity models, respectively. ‘No Q’ and ‘Q∼15’ indicate synthetics computed with weak and strong attenuation models, respectively. 1-D RSGF decay curve at 0.72 Hz used the AVG1 model unlike other 1-D model decay curves that use the REF model. (b) TT component NGF and SGF FSAs at 0.25 Hz as a function of interstation distance. The NGF FSAs (grey +) are uniformly scaled up by a constant such that the maximum value is 0. The frequency is indicated at the top left corner of the plots. Black diamonds and error bars represent the mean of NGF amplitudes and their standard deviation ± 1 σ in bins for all bins with more than five data points. The bin widths are 1 and 2 km for distances ≤10 km and >10 km, respectively. ‘SGF 3D’ and ‘SGF 1-D’ are FSAs of equivalent SGFs of the 3-D and 1-D REF velocity models, respectively, plotted for the same station pairs. The FSAs of SGFs are uniformly scaled by a constant to minimize the L1 norm of their difference with the binned NGF FSAs. (c) Similar to (b) but for FSAs at 0.72 Hz. We also plot synthetic FSAs of weak and strong attenuation models. The SGF decay curves are scaled such that their amplitude at ∼1.5 km is similar to the average amplitude of NGF FSAs in the 1–2 km range. Scaling the binned NGF decay curve instead, to directly fit the REF model SGF decay curve by minimizing the L1 norm (black dashed line), doesn't change the absolute amplitudes significantly. The RSGF decay curve (red dashed line) is scaled to fit the REF model SGF decay curve at distances ≳ 6 km. (d) The data (green +) are tangential component FSAs of earthquake records at 0.72 Hz and the red diamonds and error bars represent the mean amplitudes and their standard deviation ± 1 σ in the bins. For comparison, we also plot binned TT component NGF FSAs at 0.72 Hz (black diamonds and error bars; same as [c]) that are scaled to minimize L1 norm of binned FSAs at distances ≳ 3 km (∼1 wavelength). ‘N’ at the top right corner indicates the number of data points. Fig 16(b) compares NGF and SGF FSAs at 0.25 Hz. Since the absolute amplitudes of the NGFs are unknown, we scale the FSA decay curves by increasing or decreasing all the data in the log scale by some constant for comparison. The shape of the decay of the SGF amplitudes is different from the decay seen for the RSGF amplitudes, because of the amplitude contamination as explained previously. The degree of contamination will be lower for higher frequencies and for gently decaying band-limiting tapers applied to cross-correlation spectral amplitudes such as a cosine taper. While there is considerable scatter in the NGF amplitudes, the mean amplitudes binned by distance show systematic decay with increasing distance. The variability in the amplitudes, indicated by their standard deviation in the bins (∼0.3 in log10 units or a factor of ∼2), is likely caused by significant differences in time periods over which daily NGFs of different station pairs were stacked (∼25 d to ∼3 yr) or by effects of inhomogeneous noise source distribution or even by variability in the coupling of the sensors with the ground at different stations. For example, in Fig. 16(b), amplitudes of NGFs for pairs with station JKR are systematically lower than mean amplitudes at the same distances. The amplitude decay characteristics of the NGFs provide useful information specifically for surface sources and they are expected to be free from earthquake radiation and source effects (but not from effects of inhomogeneous noise source distributions) that are likely present in earthquake ground motions at low frequencies and small epicentral distances. While the NGF amplitudes generally exhibit similar decay behaviour as the SGF amplitudes, they appear to be systematically higher at r ≳ 10 km. Fig. 16(c) compares SGF FSAs at 0.72 Hz for a variety of 1-D and 3-D models with NGF FSAs. In this figure, the decay curves are scaled such that the mean of the NGF FSAs at 1–2 km distance is similar to the SGF FSAs of different models at ∼1.5 km. The decay of mean NGF amplitudes binned by distance is inside the domain spanned by the SGF amplitudes for 1-D velocity models with weak and strong anelastic attenuation. The presence of NGF amplitudes well outside this domain would have led to the obvious conclusion that amplitude decay characteristics obtained from coherency of the ambient noise at The Geysers are physically unrealistic and incorrect. If we assume that the relative interstation NGF amplitudes are correct and the REF model is a realistic 1-D representation of the velocity structure at The Geysers, the difference between the average NGF FSAs binned by distance and the SGF FSAs for the weak attenuation REF model at f = 0.72 Hz can be attributed to anelastic attenuation. Fitting this difference in amplitudes by a factor of e− αr, where the attenuation coefficient α is defined as $$\alpha \ = \frac{{\pi f}}{{Q( f )U( f )}}\ $$, we obtain α ∼ 0.03 km–1 and Q ∼ 33 at 0.72 Hz (in Herrmann 2013a, the attenuation coefficient's symbol is γ instead of α). The corresponding values of α and Q at 0.42 Hz are ∼0.02 km–1 and ∼27, respectively. By employing synthetic amplitudes from a realistic low-attenuation model as reference, we avoid any assumption regarding the functional form of geometrical spreading which can strongly depend on crustal structure (Burger et al.1987; Bowman & Kennett 1991). Notwithstanding the scatter in the NGF FSAs, we find satisfactory agreement between the average NGF and the REF model SGF FSA decay curves at 0.72 Hz in Fig. 16(c). Scaling the average NGF decay curve to directly fit the REF-model SGF decay curve by minimizing the L1 norm between the amplitudes does not lead to any significant change in amplitudes. We reach similar conclusions for FSAs at 0.42 Hz (Supporting Information Fig. S10a). Measured values of α and Q are also comparable to Love wave α values predicted by the REF-model (∼0.022 and ∼0.04 km–1 at 0.42 and 0.72 Hz, respectively) and constant S-wave Q values in the REF-model (∼25 to ∼40 in the top ∼4.5 km), respectively. The similarity between the RSGF and the SGF FSAs at distances r≳ 3 km shows that the contamination of the NGFs from the acausal amplitudes of the anti-causal component can be ignored at higher frequencies and large distances. For FSAs at 0.42 Hz, they are similar at distances r ≳ 8 km (Supporting Information Fig. S10a). Comparison with earthquake ground motion amplitudes We also examine the path attenuation of horizontal component ground motions of earthquakes at The Geysers. The details of the processing of earthquake records are described in  Appendix A3. Supporting Information Fig. S11 shows examples of waveforms and Fourier amplitude spectra of an earthquake. Fig. 16(d) shows tangential component earthquake FSAs at 0.72 Hz normalized by seismic moments of individual earthquakes calculated using MD–M0 relationships valid for central California (Bakun 1984). At these low frequencies, MW < 3 earthquakes can be treated as point sources and effects of source spectrum shape can be ignored. The earthquake ground motion amplitudes are not theoretically comparable to NGF amplitudes as they are a function of force-couple GFs (unlike single force GFs recovered by ambient noise cross-correlation), earthquake radiation patterns, earthquake depths and station azimuths. At distances ≳ 1λ − 2λ, the seismic phases that dominate waveforms of empirical surface-focus GFs and shallow earthquakes should be similar. Therefore, if the empirical GFs have been correctly retrieved from ambient noise cross-correlations in terms of the relative amplitudes, they should exhibit similar amplitude decay behaviour as shallow-earthquake ground motions at far-field distances (e.g. Zhang & Yang 2013). Scattering at high frequencies and averaging across multiple paths and azimuths alleviates the effects of earthquake source radiation pattern. Most of the earthquakes used in our study are very shallow (94 out of 121 earthquakes are < 3 km deep; all earthquakes < 4 km deep; Supporting Information Table S3) with respect to the shortest wavelength examined in this section (λ ∼ 3.2 km at ∼0.72 Hz). In Fig. 16(d), significant variability (standard deviation ∼0.3 log10 units) can be observed in earthquake FSAs similar to other studies (e.g. Atkinson 2004). It is likely that many factors including effects of 3-D structure, site amplification, and location and magnitude uncertainties (∼0.2 units; Peggy Hellweg, private communication, 2017) contribute to the scatter observed in the amplitudes. However, we are find that the mean earthquake FSAs binned by distance show similar decay as the NGF FSAs at 0.72 for distances ≳ 3 km. This provides strong evidence that our NGFs are recovering realistic path attenuation of ground motions. At closer distances, earthquake and NGF amplitudes are not comparable, because NGF amplitudes are contaminated and earthquake ground motion amplitudes are strongly dependent on depth and the source radiation pattern. Mean FSAs at slightly lower frequency 0.42 Hz seem to agree with each other at distances ≳ 4 km (Supporting Information Fig. S10b). Following the large difference between the shape of the TT component SGF and RSGF attenuation curves at 0.25 Hz (Figs 16a and b), we do not attempt the comparison between the tangential component earthquake and TT component NGF attenuation curves at 0.25 Hz. RR component amplitudes We also compare the RR component SGF and RSGF FSAs with NGF FSAs at 0.25 and 0.42 Hz (see Supporting Information section RR Amplitudes, Supporting Information Figs S12 and S13). The attenuation curves decay some oscillatory features that are likely an artefact of the Fourier Transform of the near-field term in the solution of elastic wave equation for a radial force and radial velocity response (Aki & Richards 2002). While there is considerable scatter in the NGF amplitudes, the average RR component NGF decay curves at 0.25 and 0.42 Hz indicate at least partial recovery of the near-field term. We recommend that future noise cross-correlation studies, especially ones involving broad-band sensors at short interstation distances, should explore any evidence for near-field terms. Bias in NGF amplitudes Accuracy of relative amplitude information retrieved from ambient noise cross-correlations has been widely discussed and debated in both theoretical and numerical studies (e.g. Cupillard & Capdeville 2010; Tsai 2011; Lawrence et al.2013). Commonly, a wave equation solution of the form $${e^{ - \alpha r}}{J_0}( {\frac{{\omega r}}{c}} )$$ is assumed in which J0 is the zero-order Bessel function of the first kind with c being the frequency-dependent phase velocity. This solution is valid for surface waves that are usually retrieved from noise cross-correlation, and consists of an implicit geometrical spreading factor in the Bessel function (∝r− 0.5 in the far-field) and attenuation in the form of e− αr. For intrinsic attenuation in a homogenous medium, Tsai (2011) showed that retrieval of this solution from coherency measurements is possible only under the conditions of a uniform noise source distribution. The ambient noise source distribution at The Geysers is evidently not uniform as none of our NGFs are symmetric. In case that the noise source distribution is not uniform but is ubiquitous (present both inside and outside an array), then azimuthal averaging of coherency over same-distance station pairs may provide correct relative amplitudes (e.g. Lawrence et al.2013; Zhang & Yang 2013; Weemstra et al.2015). Whatever the source may be, the scatter in the observed NGF amplitudes requires some averaging for meaningful interpretation, which is also true for earthquake ground motions on the same spatial scale (Fig. 16d). Otherwise, to avoid loss of local resolution from averaging, one can resort to more sophisticated techniques such as the C3 method for analysing NGF amplitudes at receivers along a line (Zhang & Yang 2013). In case of far-field noise source distribution, which is expected for continental regions in primary and secondary microseismic passbands (Stehly et al.2006), correct intrinsic attenuation cannot be retrieved even if the illumination is uniform from all directions (Tsai 2011; Weemstra et al.2015). In our study, the amplitude decay shown by the TT component RSGFs for the weak attenuation REF model at 0.72 Hz can be assumed to be equivalent to geometrical spreading, which gives us r− 0.53 for 1λ ≲ r ≲ 13 km (λ ∼ 3.5 km for cL ∼2.5 km s−1 at 0.72 Hz) and r− 1.17 for r ≳ 16 km, which are steeper than r− 0.5 (similar to SGF in Fig. 16c). Nevertheless, we compare our results for the TT component NGFs at 0.72 Hz with various conclusions drawn analytically by Tsai (2011) for different noise source distributions, assuming plane-wave incidence (r ≫ c/ω) and weak intrinsic attenuation (Supporting Information Fig. S14). We do not attempt to correct for deviations from the analytical formulations arising from orientations of components and polarizations of incident waves or effects of near-field distances (Aki 1957; Haney et al.2012). Since we observe a decay that is stronger than the expected G(r), we rule out the possibility of a single point noise source distribution that would have led to little or no decay in amplitudes with distance (Cupillard & Capdeville 2010; table 1 in Tsai 2011). In case of uniform or one-sided far-field noise source distribution, a stronger than G(r) decay is expected, which is observed in our NGF amplitudes (Cupillard & Capdeville 2010; eqs 25 and 29 in Tsai 2011). The apparent attenuation factors for the two distributions are (I0(αr))− 1 and $${( {I_0^2( {\alpha r} ) - L_0^2( {\alpha r} )} )^{ - 0.5}}$$, respectively, where I0 is the modified Bessel function of first kind and order 0 and L0 is the modified Struve function of order 0. The apparent attenuation decays are significantly weaker than the observed NGF amplitude decay. Given the scatter in the NGF amplitudes, it is difficult to constrain α, and it is possible that the decay observed in our NGF amplitudes might just be an artefact of a non-uniform far-field noise source distribution. However, 0.72 Hz is higher than the frequency band of secondary microseisms generated near the coast (Stehly et al.2006) and our analysis of T-[R,Z] components of NGFs didn't reveal any preferred noise-source incidence direction. Another possible scenario includes that all background noise is generated from the anthropological activities in the geothermal field. While the long semi-axis of the geothermal field at r' ∼12 km is not significantly smaller than the attenuation distance 1/α ∼ 1/0.04 ∼ 25 km at 0.72 Hz, using eq. 44 in Tsai (2011) for truncated near-field source distribution (apparent attenuation ∼ $${( {1 - \frac{{{r^2}}}{{4r{'^2}}}} )^{0.5}}$$), we obtain a decay that is still weaker than the observed decay for r < 20 km. This noise source distribution also predicts little or no recovery of coherence for stations outside the geothermal field whereas we obtain robust NGFs for many station pairs located at the edge of the reservoir area (e.g. SRB-DRH in Fig. 12). We are unable to discern any obvious stronger-than-actual amplitude decay at distances ≲ 2λ expected from eq. (1) in which ensemble averaging is done after spectral normalization as opposed to ensemble averaging the cross-spectrum and the amplitude spectra separately prior to the normalization (Tsai 2011; Weemstra et al.2014). It has been suggested that the approach followed in eq. (1) might be helpful if the background noise in highly non-stationary (Weemstra et al.2014). Given the difficulty in analytically estimating coherency for realistic noise source distributions, Lawrence et al. (2013) employed numerical tests and showed that azimuthal averaging of same-distance coherency measurements, similar to averaging of FSAs in our study, can yield correct attenuation estimates under a wide range of noise source distributions. We note that while averaging should be performed within spatial dimensions with slowly and smoothly varying background medium properties (Weemstra et al.2015), we average over the entire reservoir area, which doesn't seem to produce any obvious artefact in the average decay of amplitudes. Among other factors, incoherent noise locally observed at stations contributes to the autocorrelations in the denominator in coherency expression (eq. 1) and may play an important role as the contribution of coherent wavefield weakens with distance (e.g. Tsai 2011; Lawrence et al.2013). In recent studies, Weemstra et al. (2015) have shown that scattering attenuation in a non-dissipative medium, illuminated uniformly from far-field sources, can be correctly recovered from cross-spectrums. At The Geysers geothermal field, given the high temperatures (Lowenstern & Janik 2003) and considerable heterogeneities and fractures present in the subsurface (Lockner et al.1982; Thompson 1989; Gunderson 1991; O'Connell & Johnson 1991; Sammis et al.1992; Elkibbi et al.2005; Jeanne et al.2014), both intrinsic (Romanowicz 1995) and scattering attenuation are expected to be high. The QS values adopted in the REF model, ∼25 to ∼40 in the top ∼4.5 km and ≳60 at greater depths, were estimated using the NetMoment method (Hutchings 2001; Viegas & Hutchings 2011). They are considerably lower than QS values expected from standard VS-QS relationships (Brocher 2008), for example, QS ∼150–300 for VS ∼ 2.0–3.1 km s−1 in the top ∼4.5 km in the REF model, implying stronger attenuation, and possible contribution of scattering attenuation that might be better recovered under more realistic noise-source distributions. Scattering outside the study region can also act to homogenize the noise source illumination, which would aid in the better recovery of amplitudes from NGFs. In methodological aspects, some studies prefer using the same spectral normalization factors for all station pairs to obtain more appropriate relative amplitudes (e.g. Denolle et al.2013; Bowden et al.2015). However, studies applying spectral whitening at multiple stages have been successful in retrieving reliable estimates of anelastic attenuation as well (Handel et al.2016). The scatter in the NGF amplitudes observed at all frequencies necessitates a rigorous analysis of uncertainties and trade-offs in the ground motion attenuation parameters extracted from NGF FSAs, which is beyond the scope of this study. An investigation of ambient noise source distribution at The Geysers along with more NGF amplitude data are also required to properly resolve the individual contributions of possible biases, G(r) and anelastic attenuation to the overall observed decay of NGF amplitudes with distance. Notwithstanding the incompleteness of our analysis, the similarity between the average attenuation curves of the NGF FSAs, earthquake ground motion FSAs and the SGF FSAs predicted from an appropriate velocity model is very compelling and indicates that NGF amplitudes at The Geysers are believable at the distances and frequencies examined here. CONCLUSIONS AND FUTURE WORK When three-component sensors are deployed for long durations, ambient noise cross-correlation techniques provide a robust empirical, nine-component, NGF tensor that can be used in a variety of different applications. The main conclusions of our study are summarized below: We were able to retrieve NGFs in the frequency range (∼0.2–0.9 Hz) for a range of interstation distances from ∼1–30 km (∼0.22 λ–6.5 λ) at The Geysers using a variety of sensors in and around the reservoir area. For many station pairs, the NGFs are found to be similar to the single force displacement SGFs computed using pre-existing and revised 1-D velocity models in terms of waveforms, phase and the relative amplitudes of all components of the tensors, even at distances < λ. We identify both body-wave and surface-wave phases in the NGF waveforms. The direct comparison of NGFs with SGFs helps to evaluate the quality of the retrieved NGFs and the suitability of different 1-D velocity models to various sub-regions of The Geysers. SGFs computed with the faster model VSP0 and the slower model REF preferentially provide better waveform fits to the NGFs for interstation paths across SE and NW Geysers, respectively. We are also able to confirm the results of previous body-wave traveltime tomography studies such as the low velocities in the region around station GSG and higher velocities of the reservoir area relative to the surrounding region. Large anomalous amplitudes on the off-diagonal T-[R,Z] components of NGFs helped us detect sensor misalignments for many stations of The Geysers. We confirm the alignment angles estimated from the analysis of long period teleseismic waveforms by a significant reduction in these anomalous amplitudes upon rotation of the NGF tensors to correct orientations. The comparison between NGFs and SGFs computed using 1-D anisotropic velocity models and a 3-D isotropic velocity model of the reservoir suggests that robust amplitudes on the T-[R,Z] components of NGFs for some longer distance paths likely result from wave propagation effects caused by strong heterogeneity in 3-D structure. The 3-D model that was derived from body-wave traveltimes using an infinite frequency geometrical ray tomography approach (Gritto et al.2013a) can fit low frequency NGF (0.2–0.9 Hz) waveforms remarkably well. We were unable to detect any dominant ambient noise source illumination direction applying ORA to our NGFs. While there is considerable scatter in the TT component NGF FSAs, we find their average decay with distance to be similar to the decay expected from SGF amplitudes and with the decay of tangential component local-earthquake ground-motion amplitudes at the same frequencies (∼0.25–0.72 Hz). The flattening of RR component FSAs at distances ∼9–16 km suggests possible recovery of the near-field term. The similarity of the NGF and the SGF waveforms computed with appropriate velocity models in this study indicates that the full nine-component NGF tensor should be used in waveform tomography studies (e.g. Lee et al.2014) whenever multicomponent stations are available in dense networks. As demonstrated, higher amplitude TT components provide strong reliable constraints on VS structure. ZR and RZ component NGFs are possibly less susceptible to effects of directional noise source incidence compared to the commonly used ZZ component. Broad-band phase delays from NGF waveforms as in the case of the NGFs to station GSG can be potentially useful at short interstation distances. Calibration of starting velocity models using NGF waveforms rather than earthquake waveforms can be desirable in regions like The Geysers, where earthquakes might have non-double couple source mechanisms (Guilhem et al.2014). At The Geysers, possibilities of future work include further refinement of the velocity models, monitoring of temporal changes in the coda of NGFs and characterization of background noise source distribution. Given the excellent distribution of earthquakes and stations within the geothermal field at The Geysers, the 3-D velocity model derived from body-wave tomography is well resolved in the reservoir area; therefore, it is generally successful fitting low-frequency NGF waveforms. However, the model is not well resolved at depths ≳ 4.2 km owing to the shallow focal depths of the earthquakes ( ≲ 4.5 km) or at very shallow depths ( ≲ 1.0 km). Long period Rayleigh waves (∼5–8 s; cR ∼ 2.6–2.9 km s−1 for the GIL7 model) derived from noise cross-correlation at interstation distances (∼40–70 km; r ≳ 3λ) can be used to investigate the velocity structure of the deeper reservoir rocks and the underlying felsite (depth sensitivity ∼λ/2 ∼ 6–12 km). We have demonstrated good results with some of the short period (∼1.0 Hz) sensors of the NCSN that are located outside The Geysers. Retrieval of longer period measurements will also benefit from temporary broad-band sensors deployed at The Geysers (Specht et al.2014) that are expected to yield more stable and higher quality NGFs than the ∼4.5 Hz geophones used in this study. We didn't find any significant energy above 1.0 Hz in our NGFs and the presence of multiple spectral peaks (see Appendix A2) between 1.7–2.5 Hz precluded any meaningful interpretation at these frequencies. Therefore, retrieval of high frequency (≳ 1 Hz) empirical GFs from cross-correlation of coda waves of the numerous local earthquakes in the reservoir area can be attempted (e.g. Campillo & Paul 2003) to investigate shallow velocity structure (cR ∼ 2.2 km s−1 at ∼1.0 Hz for REF model; depth sensitivity ∼λ/2 ∼ 1.1 km). A refined version of the 3-D velocity model can subsequently be applied to higher frequency (up to ∼2.5 Hz) earthquake waveforms and should increase confidence in source inversion studies by facilitating the inclusion of more stations at longer distances (e.g. Guilhem et al.2014). For any future study attempting improvement of velocity models at The Geysers, the results obtained in this study provide a framework for appropriate initial models for inversion. Possible temporal changes in the coda of noise cross-correlations can be investigated to look for possible changes in subsurface VS related to micro/macroseismicity or injection/production activities. Variations of the order of few per cent in VS similar to ∼4–6 per cent variations in VP/VS ratio inferred from earthquake body-wave travel times (Foulger et al.1997; Gritto & Jarpe 2014) should be readily detectable in the coda of NGFs. However, it must be noted that many stations don't return stable or robust NGFs during the entire time period. Characterization of ambient noise source distribution and incidence direction in the frequency range ∼0.2–0.9 Hz at The Geysers and its temporal variation using beamforming analysis, etc. (e.g. Menon et al.2014) is an important task. As discussed previously, ambient noise source distribution has important implications for the decay of NGF amplitudes and amplitudes in the T-[R,Z] components of NGFs. It would be interesting to see if effects on NGFs expected from the inferred noise source distribution (Fichtner 2014) can be reconciled with our observations. Acknowledgements A.N. is indebted to R.B. Herrmann (St. Louis University) and F.-C. Lin (University of Utah) for their help and guidance on many noise cross-correlation and surface wave concepts, and for reviewing the first draft of this paper. We thank an anonymous reviewer, D.C. Bowden (Caltech), M.A. Denolle (Harvard University) and Editor L. Métivier for carefully reviewing this paper, and for their comments and suggestions, which helped to significantly improve this study. We are thankful to V.H. Lai (Caltech), P. Hellweg [Berkeley Seismological Lab (BSL)], H. Rademacher (BSL) and V.C. Tsai (Caltech) for helpful comments and/or discussions. We also thank Ramsey Haught for visiting the stations and checking the sensor orientations. A.N. thanks V. Maupin (University of Oslo) for pointing out the reference for the elastic tensor for crack-induced anisotropy. We are grateful for permission to use the Baribu cluster at BSL for running SW4. This study was supported by National Science Foundation award EAR-1053211 and partially supported by the France-Berkeley Fund 2014-0051. Data used in this study come from Berkeley Digital Seismic Network (BDSN; dx.doi.org/10.7932/BDSN), operated by the University of California Berkeley Seismological Laboratory, Northern California Seismic Network (NCSN), and Lawrence Berkeley National Laboratory (LBNL) Short Period Network at The Geysers, which are all archived at the Northern California Earthquake Data Center (NCEDC; dx.doi.org/10.7932/NCEDC; http://www.ncedc.org, last accessed July 2016). Codes for computing SGFs for 1-D velocity models by FKI and the modal summation method, codes for extracting Love wave group and phase velocities from NGFs and codes for forward modelling the same using 1-D velocity models are all available in the software Computer Programs in Seismology available at http://www.eas.slu.edu/eqc/eqccps.html, last accessed June 2016. SW4 is hosted by the Computational Infrastructure for Geodynamics (CIG) which is supported by the National Science Foundation award NSF-0949446. SW4 can be downloaded from CIG website https://geodynamics.org/cig/software/sw4/, last accessed June 2016. REFERENCES Aki K., 1957. Space and time spectra of stationary stochastic waves, with special reference to microtremors, Bull. Earthq. Res. Inst. Univ. Tokyo , 25, 415– 457. Aki K., Richards P.G., 2002. Quantitative Seismology , 2nd edn, University Science Books, pp. 27– 72. Allis R.G., Shook G.M., 1999. An alternative mechanism for the formation of The Geysers vapor-dominated reservoir, in Proc. 24th Workshop on Geotherm. Reservoir Engg. , Stanford University, California. Atkinson G., 2004. Empirical attenuation of ground-motion spectral amplitudes in Southeastern Canada and the Northeastern United States, Bull. seism. Soc. Am. , 94( 3), 1079– 1095. https://doi.org/10.1785/0120030175 Google Scholar CrossRef Search ADS   Bakun W.H., 1984. Seismic moments, local magnitudes and coda-duration magnitudes for earthquakes in Central California, Bull. seism. Soc. Am. , 74( 2), 439– 458. Bensen G.D., Ritzwoller M.H., Barmin M.P., Levshin A.L., Lin F.-C., Moschetti M.P., Shapiro N.M., Yang Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, Geophys. J. Int. , 169( 3), 1239– 1260. https://doi.org/10.1111/j.1365-246X.2007.03374.x Google Scholar CrossRef Search ADS   Bowden D.C., Tsai V.C., Lin F.-C., 2015. Site amplification, attenuation, and scattering from noise correlation amplitudes across a dense array in Long Beach, CA, Geophys. Res. Lett. , 42( 5), 1360– 1367. https://doi.org/10.1002/2014GL062662 Google Scholar CrossRef Search ADS   Bowman J.R., Kennett B.L.N., 1991. Propagation of Lg waves in the North Australian Craton: influence of crustal velocity gradients, Bull. seism. Soc. Am. , 81( 2), 592– 610. Boyd O.S., Dreger D.S., Lai V.H., Gritto R., 2015. A systematic analysis of seismic moment tensor at The Geysers geothermal field, California, California, Bul. seism. Soc. Am. , 105( 6), 2969– 2986. https://doi.org/10.1785/0120140285 Google Scholar CrossRef Search ADS   Brenguier F., Campillo M., Takeda T., Aoki Y., Shapiro N.M., Briand X., Emoto K., Miyake H., 2014. Mapping pressurized volcanic fluids from induced crustal seismic velocity drops, Science , 345( 6192), 80– 82. https://doi.org/10.1126/science.1254073 Google Scholar CrossRef Search ADS PubMed  Brocher T.M., 2008. Key elements of regional seismic velocity models for long period ground motion simulations, J. Seismol. , 12( 2), 217– 221. https://doi.org/10.1007/s10950-007-9061-3 Google Scholar CrossRef Search ADS   Burger R.W., Somerville P.G., Barker J.S., Herrmann R.B., Helmberger D.V., 1987. The effect of crustal structure on strong ground motion attenuation relations in Eastern North America, Bull. seism. Soc. Am. , 77( 2), 420– 439. Campillo M., Paul A., 2003. Long-range correlations in the diffuse seismic coda, Science , 299( 5606), 547– 549. https://doi.org/10.1126/science.1078551 Google Scholar CrossRef Search ADS PubMed  Cheng C.H., 1993. Crack models for a transversely isotropic medium, J. geophys. Res. , 98( B1), 675– 684. https://doi.org/10.1029/92JB02118 Google Scholar CrossRef Search ADS   Cho K.H., Herrmann R.B., Ammon C.J., Lee K., 2007. Imaging the upper crust of the Korean peninsula by surface-wave tomography, Bull. seism. Soc. Am. , 97( 1B), 198– 207. https://doi.org/10.1785/0120060096 Google Scholar CrossRef Search ADS   Crampin S., 1984. Effective anisotropic elastic constants for wave propagation through cracked solids, Geophys. J. Int. , 76( 1), 135– 145. https://doi.org/10.1111/j.1365-246X.1984.tb05029.x Google Scholar CrossRef Search ADS   Cupillard P., Capdeville Y., 2010. On the amplitude of surface waves obtained by noise correlation and the capability to recover the attenuation: a numerical approach, Geophys. J. Int. , 181, 1687– 1700. Denolle M.A., Dunham E.M., Prieto G.A., Beroza G.C., 2013. Ground motion prediction of realistic earthquake sources using the ambient seismic field, J. geophys. Res. , 118( 5), 2102– 2118. https://doi.org/10.1029/2012JB009603 Google Scholar CrossRef Search ADS   Denolle M.A., Miyake H., Nakagawa S., Hirata N., Beroza G.C., 2014. Long-period seismic amplification in the Kanto Basin from the ambient seismic field, Geophys. Res. Lett. , 41( 7), 2319– 2325. https://doi.org/10.1002/2014GL059425 Google Scholar CrossRef Search ADS   Durand S., Montagner J.P., Roux P., Brenguier F., Nadeau R.M., Ricard Y., 2011. Passive monitoring of anisotropy change associated with the Parkfield 2004 earthquake, Geophys. Res. Lett. , 38( 13), n/a-n/a. https://doi.org/10.1029/2011GL047875 Google Scholar CrossRef Search ADS   Dziewonski A.M., Bloch S., Landisman M., 1969. A technique for the analysis of transient seismic signals, Bull. seism. Soc. Am. , 59, 427– 444. Eberhart-Phillips D., 1986. Three-dimensional velocity structure in northern California Coast Ranges from inversion of local earthquake arrival times, Bull. seism. Soc. Am. , 76( 4), 1025– 1052. Eberhart-Phillips D., Oppenheimer D.H., 1984. Induced seismicity in The Geysers geothermal area, California, J. geophys. Res. , 82( B2), 1191– 1207. https://doi.org/10.1029/JB089iB02p01191 Google Scholar CrossRef Search ADS   Elkibbi M., Rial J.A., 2005. The Geysers geothermal field: results from shear-wave splitting analysis in a fractured reservoir, Geophys. J. Int. , 162( 3), 1024– 1035. https://doi.org/10.1111/j.1365-246X.2005.02698.x Google Scholar CrossRef Search ADS   Elkibbi M., Yang M., Rial J.A., 2005. Crack-induced anisotropy models in The Geysers geothermal field, Geophys. J. Int. , 162( 3), 1036– 1048. https://doi.org/10.1111/j.1365-246X.2005.02697.x Google Scholar CrossRef Search ADS   Fichtner A., 2014. Source and processing effects on noise correlations, Geophys. J. Int. , 197( 3), 1527– 1531. https://doi.org/10.1093/gji/ggu093 Google Scholar CrossRef Search ADS   Foulger G.R., Grant C.C., Ross A., Julian B.R., 1997. Industrially induced changes in earth structure at The Geysers geothermal area, California, Geophys. Res. Lett. , 24( 2), 135– 137. https://doi.org/10.1029/96GL03152 Google Scholar CrossRef Search ADS   Grigoli F., Cesca S., Dahm T., Krieger L., 2012. A complex linear least-squares method to derive relative and absolute orientations of seismic sensors, Geophys. J. Int. , 188( 3), 1243– 1254. https://doi.org/10.1111/j.1365-246X.2011.05316.x Google Scholar CrossRef Search ADS   Gritto R., Yoo S.-H., Jarpe S.P., 2013a. 3D seismic tomography at The Geysers geothermal field, CA, USA, in Proc. Thirty-Eighth Workshop on Geothermal Reservoir Engineering , Stanford University, California, 11–13 February 2013. Gritto R., Yoo S.-H., Jarpe S., 2013b. Seismic Imaging of reservoir structure at The Geysers geothermal reservoir, in AGU Fall Meeting 2013 S33D-2460 , San Francisco, USA, 9–13 December 2013. Gritto R., Jarpe S.P., 2014. Temporal variations of VP/VS-ratio at The Geysers geothermal field, USA, Geothermics , 52, 112– 119. https://doi.org/10.1016/j.geothermics.2014.01.012 Google Scholar CrossRef Search ADS   Guilhem A., Hutchings L., Dreger D.S., Johnson L.R., 2014. Moment tensor inversions of M ∼ 3 earthquakes in the Geysers geothermal fields, California, J. geophys. Res. , 119( 3), 2121– 2137. https://doi.org/10.1002/2013JB010271 Google Scholar CrossRef Search ADS   Gunderson R.P., 1991. Porosity of reservoir greywacke at The Geysers, in Monograph on The Geysers Geothermal Field , Geothermal Resources Council Spl. Rept. No. 17, pp. 89– 93. Handel A., Ohrnberger M., Krüger F., 2016. Extracting near-surface QL between 1–4 Hz from higher-order noise correlations in the Euroseistest area, Greece, Geophys. J. Int. , 207( 2), 655– 666. https://doi.org/10.1093/gji/ggw295 Google Scholar CrossRef Search ADS   Haney M., Mikesell T.D., van Wijk K., Nakahara H., 2012. Extension of the spatial autocorrelation (SPAC) method to mixed-component correlations of surface waves, Geophys. J. Int. , 191( 1), 189– 206. https://doi.org/10.1111/j.1365-246X.2012.05597.x Google Scholar CrossRef Search ADS   Haskell N.A., 1964. Radiation pattern of surface waves from point sources in a multi-layered medium, Bull. seism. Soc. Am. , 54( 1), 377– 393. Hearn B.C.Jr, Donnelly-Nolan J.M, Goff F.E., 1981. The Clear Lake volcanics: Tectonic setting and magma sources, U.S. Geol. Surv. Prof. Pap.  1141, 25– 45. Herrmann R.B., 1973. Some aspects of band-pass filtering of surface waves, Bull. seism. Soc. Am. , 63( 2), 663– 671. Herrmann R.B., 2013a. Computer programs in seismology: An evolving tool for instruction and research, Seismol. Res. Lett. , 84( 6), 1081– 1088. https://doi.org/10.1785/0220110096 Google Scholar CrossRef Search ADS   Herrmann R.B., 2013b. Update to do_mft for the determination of phase velocities from empirical Green's functions from noise cross-correlation, http://www.eas.slu.edu/eqc/eqc_cps/TUTORIAL/EMPIRICAL_GREEN/index.html, last accessed June 2016. Hudson J.A., 1981. Wave speeds and attenuation of elastic waves in material containing cracks, Geophys. J. Int. , 64( 1), 133– 150. https://doi.org/10.1111/j.1365-246X.1981.tb02662.x Google Scholar CrossRef Search ADS   Hudson J.A., 1980. Overall properties of a cracked solid, Math. Proc. Camb. Phil. Soc. , 88( 02), 371– 384. https://doi.org/10.1017/S0305004100057674 Google Scholar CrossRef Search ADS   Hutchings L., 2001. Program NetMoment; simultaneous calculation of moment, source corner frequency, and site specific t* from network recordings, Lawrence Livermore National Laboratory Report  UCRL-ID-135693-REV-1. Jeanne P., Rutqvist J., Hartline C., Garcia J., Dobson P.F., Walter M., 2014. Reservoir structure and properties from geomechanical modeling and microseismicity analyses associated with an enhanced geothermal system at The Geysers, California, Geothermics , 51, 460– 469. https://doi.org/10.1016/j.geothermics.2014.02.003 Google Scholar CrossRef Search ADS   Johnson C.W., Totten E.J., Bürgmann R., 2016. Depth migration of seasonally induced seismicity at The Geysers geothermal field, Geophys. Res. Lett. , 43( 12), 6196– 6204 https://doi.org/10.1002/2016GL069546 Google Scholar CrossRef Search ADS   Julian B.R., Ross A., Foulger G.R., Evans J.R., 1996. Three-dimensional seismic image of a geothermal reservoir: The Geysers, California, Geophys. Res. Lett. , 23( 6), 685– 688. https://doi.org/10.1029/95GL03321 Google Scholar CrossRef Search ADS   Kawakatsu H., Montagner J.-P., Song T.-R.A., 2015. On DLA's η, in The Interdisciplinary Earth: A Volume in Honor of Don L. Anderson , Geol. Soc. Am. Spl. Pap. 514 and Amer. Geophys. Union Spl. Publ. 71, pp. 33– 38, eds Foulger G.R., Lustrino M, King S.D, Geological Society of America. Google Scholar CrossRef Search ADS   Kimman W.P., Trampert J., 2010. Approximations in seismic interferometry and their effects on surface waves, Geophys. J. Int. , 182, 461– 476. Lawrence J.F., Denolle M.A., Seats K.J., Prieto G.A., 2013. A numeric evaluation of attenuation from ambient noise correlation functions, J. geophys. Res. , 118( 12), 6134– 6145. https://doi.org/10.1002/2012JB009513 Google Scholar CrossRef Search ADS   Lee E.-J., Chen P., Jordan T.H., Maechling P.B., Denolle M.A., Beroza G.C., 2014. Full-3-D tomography for crustal structure in Southern California based on the scattering-integral and the adjoint-wavefield methods, J. geophys. Res. , 119( 8), 6421– 6451. https://doi.org/10.1002/2014JB011346 Google Scholar CrossRef Search ADS   Lin F.-C., Moschetti M.P., Ritzwoller M.H., 2008. Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps, Geophys. J. Int. , 173( 1), 281– 298. https://doi.org/10.1111/j.1365-246X.2008.03720.x Google Scholar CrossRef Search ADS   Lin F.-C., Li D., Clayton R.W., Hollis D., 2013. High-resolution 3D shallow crustal structure in Long Beach, California: application of ambient noise tomography on a dense seismic array, Geophysics , 78( 4), Q45– Q56. https://doi.org/10.1190/geo2012-0453.1 Google Scholar CrossRef Search ADS   Lin F.-C., Tsai V.C., Schmandt B., 2014. 3-D crustal structure of the western United States: application of Rayleigh-wave ellipticity extracted from noise cross-correlations, Geophys. J. Int.  198( 2), 656– 670. https://doi.org/10.1093/gji/ggu160 Google Scholar CrossRef Search ADS   Liu X., Kieffer S.W., 2009. Thermodynamics and Mass Transport in Multicomponent, Multiphase H 2 O Systems of Planetary Interest, Annu. Rev. Earth Planet. Sci. , 37( 1), 449– 477. https://doi.org/10.1146/annurev.earth.031208.100109 Google Scholar CrossRef Search ADS   Lockner D.A., Summers R., Moore D., Byerlee J.D., 1982. Laboratory Measurements of Reservoir Rock from the Geysers Geothermal Field, California, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. , 19( 2), 65– 80. https://doi.org/10.1016/0148-9062(82)91632-1 Google Scholar CrossRef Search ADS   Lowenstern J.B., Janik C.J., 2003. The origins of reservoir liquids and vapors from The Geysers geothermal field, California (USA), in Volcanic, Geothermal and Ore-forming Fluids: Rulers and Witnesses of Processes within the Earth , pp. 181– 195, eds Simmons S.F., Graham I., Society of Economic Geologists Special Publication 10. Ludwin R.S., Cagnetti V., Bufe C.G., 1982. Comparison of seismicity in The Geysers geothermal area with the surrounding region, Bull. seism. Soc. Am. , 72( 3), 863– 871. Ma S., Prieto G.A., Beroza G.C., 2008. Testing community velocity models for southern California using the ambient seismic field, Bull. seism. Soc. Am. , 98( 6), 2694– 2714. https://doi.org/10.1785/0120080947 Google Scholar CrossRef Search ADS   Majer E.L., McEvilly T.V., 1979. Seismological investigations at The Geysers geothermal field, Geophysics , 44( 2), 246– 269. https://doi.org/10.1190/1.1440965 Google Scholar CrossRef Search ADS   Majer E.L., Peterson J.E., 2007. The impact of injection on seismicity at The Geysers, California geothermal field, Int. J. Rock Mech. Min. Sci. , 44( 8), 1079– 1090. https://doi.org/10.1016/j.ijrmms.2007.07.023 Google Scholar CrossRef Search ADS   Majer E.L., McEvilly T.V., Eastwood F.S., Myer L.R., 1988. Fracture detection using P-wave and S-wave vertical seismic profiling at The Geysers, Geophysics , 53., 76– 84. https://doi.org/10.1190/1.1442401 Google Scholar CrossRef Search ADS   Maupin V., Park J., 2007. Theory and observations—wave propagation in anisotropic media, in Treatise on Geophysics, Vol. 1: Seismology and the Structure of the Earth , . 289– 321, eds Romanowicz B., Dziewonski A., Elsevier. Menon R., Gerstoft P., Hodgkiss W.S., 2014. On the apparent attenuation in the spatial coherence estimated from seismic arrays, J. geophys. Res. , 119( 4), 3115– 3132. https://doi.org/10.1002/2013JB010835 Google Scholar CrossRef Search ADS   Montagner J.P., Anderson D.L., 1989. Petrological constraints on seismic anisotropy, Phys. Earth planet. Inter. , 54( 1–2), 82– 105. https://doi.org/10.1016/0031-9201(89)90189-1 Google Scholar CrossRef Search ADS   Mossop A., Segall P., 1999. Volume strain within The Geysers geothermal field, J. geophys. Res. , 104( B12), 29 113–29 131. https://doi.org/10.1029/1999JB900284 Google Scholar CrossRef Search ADS   Nakata N., Snieder R., Tsuji T., Larner K., Matsuoka T., 2011. Shear wave imaging from traffic noise using seismic interferometry by cross-coherence, Geophysics , 76( 6), SA97– SA106. https://doi.org/10.1190/geo2010-0188.1 Google Scholar CrossRef Search ADS   Nakata N., Chang J., Lawrence J.F., Boué P., 2015. Body wave extraction and tomography at Long Beach, California, with ambient-noise interferometry, J. geophys. Res. , 120( 2), 1159– 1173. https://doi.org/10.1002/2015JB011870 Google Scholar CrossRef Search ADS   Nishida K., Kawakatsu H., Obara K., 2008. Three-dimensional crustal S wave velocity structure in Japan using microseismic data recorded by Hi-net tiltmeters, J. geophys. Res. , 113( B10), B10302. https://doi.org/10.1029/2007JB005395 Google Scholar CrossRef Search ADS   O'Connell D.R.H., Johnson L.R., 1991. Progressive inversion for hypocenters and P wave and S wave velocity structure: Application to the Geysers, California, Geothermal Field, J. geophys. Res. , 96( B4), 6223– 6236. https://doi.org/10.1029/91JB00154 Google Scholar CrossRef Search ADS   Oppenheimer D.H., 1986. Extensional tectonics at The Geysers geothermal area, California, J. geophys. Res. , 91( B11), 11 463–11 476. https://doi.org/10.1029/JB091iB11p11463 Google Scholar CrossRef Search ADS   Peacock S., Hudson J.A., 1990. Seismic properties of rocks with distributions of small cracks, Geophys. J. Int. , 102( 2), 471– 484. https://doi.org/10.1111/j.1365-246X.1990.tb04479.x Google Scholar CrossRef Search ADS   Petersson N.A., Sjögreen B., 2014a. SW4-v1.1 Users Guide , https://geodynamics.org/cig/software/sw4/, last accessed, December 2015. Petersson N.A., Sjögreen B., 2014b. Super-grid modeling of the elastic wave equation in semi-bounded domains, Commun. Commut. Phys. , 16( 04), 913– 955. https://doi.org/10.4208/cicp.290113.220514a Google Scholar CrossRef Search ADS   Porritt R.W., Allen R.M., Boyarko D.C., Brudzinski M.R., 2011. Investigation of Cascadia segmentation with ambient noise tomography, Earth planet. Sci. Lett. , 309(1–2), 67– 76. https://doi.org/10.1016/j.epsl.2011.06.026 Google Scholar CrossRef Search ADS   Prieto G.A., Beroza G.C., 2008. Earthquake ground motion prediction using the ambient seismic field, Geophys. Res. Lett. , 35( 14), L14304. https://doi.org/10.1029/2008GL034428 Google Scholar CrossRef Search ADS   Prieto G.A., Denolle M.A., Lawrence J.F., Beroza G.C., 2011. On amplitude information carried by the ambient seismic field, C. R. Geosci. , 343( 8–9), 600– 614. https://doi.org/10.1016/j.crte.2011.03.006 Google Scholar CrossRef Search ADS   Romanowicz B., 1995. A global tomographic model of shear attenuation in the upper mantle, J. geophys. Res. , 100( B7), 12 375–12 394. https://doi.org/10.1029/95JB00957 Google Scholar CrossRef Search ADS   Ross A., Foulger G.R., Julian B.R., 1999. Source processes of industrially-induced earthquakes at The Geysers geothermal area, California, Geophysics  64( 4), 1877–1889. Roux P., 2009. Passive seismic imaging with directive ambient noise: application to surface waves and the San Andreas Fault in Parkfield, CA, Geophys. J. Int. , 179( 1), 367– 373. https://doi.org/10.1111/j.1365-246X.2009.04282.x Google Scholar CrossRef Search ADS   Saade M., Montagner J.P., Roux P., Cupillard P., Durand S., Brenguier F., 2015. Influence of seismic anisotropy on the cross correlation tensor: numerical investigations, Geophys. J. Int. , 201( 2), 595– 604. https://doi.org/10.1093/gji/ggu470 Google Scholar CrossRef Search ADS   Sammis C.G., An L., Ershaghi I., 1992. Determining the 3-D fracture structure in the Geysers geothermal reservoir, in Proc. Seventeenth Workshop on Geothermal Reservoir Engineering , Stanford University, California, 29–31 January 1992. Seats K.J., Lawrence J.F., Prieto G.A., 2012. Improved ambient noise correlation functions using Welch's method, Geophys. J. Int. , 188( 2), 513– 523. https://doi.org/10.1111/j.1365-246X.2011.05263.x Google Scholar CrossRef Search ADS   Sens-Schönfelder C., 2008. Synchronizing seismic networks with ambient noise, Geophys. J. Int. , 174( 3), 966– 970. https://doi.org/10.1111/j.1365-246X.2008.03842.x Google Scholar CrossRef Search ADS   Sens-Schönfelder C., Wegler U., 2006. Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia, Geophys. Res. Lett. , 33( 21), L21302. https://doi.org/10.1029/2006GL027797 Google Scholar CrossRef Search ADS   Shapiro N.M., Campillo M., 2004. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise, Geophys. Res. Lett. , 31( 7), L07614. https://doi.org/10.1029/2004GL019491 Google Scholar CrossRef Search ADS   Shapiro N.M., Campillo M., Stehly L., Ritzwoller M.H., 2005. High-resolution surface-wave tomography from ambient seismic noise, Science , 307( 5715), 1615– 1618. https://doi.org/10.1126/science.1108339 Google Scholar CrossRef Search ADS PubMed  Sjögreen B., Petersson N.A., 2012. A fourth order accurate finite difference scheme for the elastic wave equation in second order formulation, J. Sci. Comput. , 52( 1), 17– 48. https://doi.org/10.1007/s10915-011-9531-1 Google Scholar CrossRef Search ADS   Snieder R., 2004. Extracting the Green's function from the correlation of coda waves: a derivation based on stationary phase, Phys. Rev. E , 69( 4), 046610. https://doi.org/10.1103/PhysRevE.69.046610 Google Scholar CrossRef Search ADS   Specht S., Jousset P., Zang A., Gritto R., Bruhn D., 2014. Velocity structure of The Geysers geothermal area (California) from ambient noise cross-correlation, in AGU Fall Meeting 2014 S51A-4425 , San Francisco, USA, 15–19 December 2014. Stark M.A., 1991. Microearthquakes—A tool to track injected water in The Geysers reservoir, in Monograph on The Geysers Geothermal Field , Geothermal Resources Council Spl. Rept. No. 17, pp. 111– 117. Stehly L., Campillo M., Shapiro N.M., 2006. A study of the seismic noise from its long-range correlation properties, J. geophys. Res. , 111( B10), B10306. https://doi.org/10.1029/2005JB004237 Google Scholar CrossRef Search ADS   Stidham C., Antolik M., Dreger D.S., Larsen S., Romanowicz B., 1999. Three-dimensional structure influences on the strong motion wavefield of the 1989 Loma Prieta earthquake, Bull. seism. Soc. Am. , 89, 1184– 1202. Takagi R., Nishida K., Aoki Y., Maeda T., Masuda K., Takeo M., Obara K., Shiomi K., Sato M., Saito K., 2015. A single bit matters: coherent noise of seismic data loggers, Seismol. Res. Lett. , 86( 3), 901– 907. https://doi.org/10.1785/0220150030 Google Scholar CrossRef Search ADS   Taira T., Brenguier F., Kong Q., 2015. Ambient noise-based monitoring of seismic velocity changes associated with the 2014 M w 6.0 South Napa earthquake, Geophys. Res. Lett. , 42( 17), 6997– 7004. https://doi.org/10.1002/2015GL065308 Google Scholar CrossRef Search ADS   Thompson R.C., 1989. Structural stratigraphy and intrusive rocks at The Geysers geothermal field, Trans. Geotherm. Resour. Counc. , 13, 481– 485. Truesdell A.H., Haizlip J.R., Box W.TJr., Amore F.D’, 1991. A geochemical overview of The Geysers geothermal reservoir, in Monograph on The Geysers Geothermal Field , Geothermal Resources Council Spl. Rept. No. 17, 121– 132. Trugman D.T., Shearer P.M., Borsa A.A., Fialko Y., 2016. A comparison of long-term changes in seismicity at The Geysers, Salton Sea, and Coso geothermal fields, J. geophys. Res. , 121( 1), 225– 247. https://doi.org/10.1002/2015JB012510 Google Scholar CrossRef Search ADS   Tsai V.C., 2009. On establishing the accuracy of noise tomography travel-time measurements in a realistic medium, Geophys. J. Int. , 178( 3), 1555– 1564. https://doi.org/10.1111/j.1365-246X.2009.04239.x Google Scholar CrossRef Search ADS   Tsai V.C., 2011. Understanding the amplitudes of noise correlation measurements, J. geophys. Res. , 116( B9), B09311, doi:10.1029/2011JB008483. https://doi.org/10.1029/2011JB008483 Google Scholar CrossRef Search ADS   Tsai V.C., Moschetti M.P., 2010. An explicit relationship between time-domain noise correlation and spatial autocorrelation (SPAC) results, Geophys. J. Int , 182, 454– 460. Google Scholar CrossRef Search ADS   van Wijk K., Mikesell T.D., Schulte-Pelkum V., Stachnik J., 2011. Estimating the Rayleigh-wave impulse response between seismic stations with the cross terms of the Green tensor, Geophys. Res. Lett. , 38( 16), L16301. https://doi.org/10.1029/2011GL047442 Google Scholar CrossRef Search ADS   Viegas G., Hutchings L., 2011. Characterization of Induced Seismicity near an Injection Well at the Northwest Geysers Geothermal Field, California, Trans. Geotherm. Resour. Counc. , 35, 1773– 1780. Wang C.Y., Herrmann R.B., 1980. A numerical study of P-, SV- and SH-wave generation in a plane layered medium, Bull. seism. Soc. Am. , 70( 4), 1015– 1036. Wapenaar K., Fokkema J., 2006. Green's function representations for seismic interferometry, Geophysics  71( 4), SI33– SI46, doi:10.1190/1.2213955. Google Scholar CrossRef Search ADS   Weemstra C., Westra W., Snieder R., Boschi L., 2014. On estimating attenuation from the amplitude of the spectrally whitened ambient seismic field, Geophys. J. Int. , 197( 3), 1770– 1788. https://doi.org/10.1093/gji/ggu088 Google Scholar CrossRef Search ADS   Weemstra C., Snieder R., Boschi L., 2015. On the estimation of attenuation from the ambient seismic field: inferences from distributions of isotropic point scatterers, Geophys. J. Int. , 203( 2), 1054– 1071. https://doi.org/10.1093/gji/ggv311 Google Scholar CrossRef Search ADS   Wegler U., Nakahara H., Sens-Schonfelder C., Korn M., Shiomi K., 2009. Sudden drop of seismic velocity after the 2004 M w 6.6 mid-Niigata earthquake, Japan, observed with Passive Image Interferometry, J. geophys. Res. , 114( B6), B06305. https://doi.org/10.1029/2008JB005869 Google Scholar CrossRef Search ADS   Wielandt E., Forbriger T., 1999. Near-field seismic displacement and tilt associated with the explosive activity of Stromboli, Ann. Geofis. , 42( 3), 407– 416. Yao H., van der Hilst R.D., de Hoop M.V., 2006. Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis—I. Phase velocity maps, Geophys. J. Int. , 166, 732– 744. https://doi.org/10.1111/j.1365-246X.2006.03028.x Google Scholar CrossRef Search ADS   Yokoi T., Margaryan S., 2008. Consistency of the spatial autocorrelation method with seismic interferometry and its consequence, Geophys. Prospect. , 56, 435– 451. https://doi.org/10.1111/j.1365-2478.2008.00709.x Google Scholar CrossRef Search ADS   Zhang J., Yang X., 2013. Extracting surface wave attenuation from seismic noise using correlation of the coda of correlation, J. geophys. Res. , 118, 2191– 2205. https://doi.org/10.1002/jgrb.50186 Google Scholar CrossRef Search ADS   Zhan Z., Ni S., 2010. Stationary phase approximation in the ambient noise method revisited, Earthq. Sci. , 23, 425– 431. https://doi.org/10.1007/s11589-010-0741-7 Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1. Filtered waveforms (black, reference station GDXB; grey, other stations) of two earthquakes used to determine the sensor misalignments at The Geysers. The earthquake, year, MW and frequency passband are indicated at the top. The stations used for the two earthquakes (grey traces) are (a) ACR, AL6, CLV, DRK, DVB, FUM, HVC, MNS, TCH, JKB and DRH, and (b) ACR, CLV, DRK, DVB, FUM, HVC, MNS, NEG, TCH and DRH. Original waveforms are in the top panels, while the bottom panels show waveforms after they have been time-shifted and horizontal components have been rotated to increase similarity with GDXB waveforms that remain fixed. Figure S2. VR between waveforms of large earthquakes recorded at three-component stations with corresponding waveforms at USGS broad-band station GDXB as a function of anticlockwise rotation in azimuth of that station. Black solid lines are curves for individual phases of the same or different earthquakes and the red crosses show coordinates of the maximum VR of each curve. The grey solid and dashed lines mark the mean rotation angle (of all the red crosses) and ± 1 standard deviation. The station name, mean angle and standard deviation, number of phases used and the average spread (angle from the VR maxima at which VR drops by 5 per cent) are indicated. Figure S3. Figure explaining the processing of synthetics in this study. (a) Raw TT component synthetic velocity Green's function (GF) (velocity response to an input step function) at interstation distance ∼2 km computed using the REF model. We use this raw velocity time-series as raw synthetic Green's function (RSGF) for amplitude decay analysis in the section Amplitude Decay in the main text. (b) Time-series in (a) integrated to displacement and bandpass-filtered between 0.2 and 1.6 Hz with a zero-phase eight-pole Butterworth filter (black trace) and its time-reversed version (red trace). The acausal filtering is the effect of the frequency domain taper applied to keep cross-correlation spectral amplitudes band-limited (see subsection Cross-correlation analysis in the Methodology section in the main text). (c) Sum of causal and anticausal (time-reversed) traces. Ideally, this is the perfect noise-derived Green's function (NGF) we expect to recover from noise cross-correlation for our choice of parameters. (d) Symmetric component (average of causal and anti-causal sides) of the trace in (c) as an estimate of synthetic Green's function (SGF; black trace). Comparing against the original filtered displacement time-series (red trace), the amplitudes of SGF at times ≲ 3 s are contaminated by the acausal amplitudes of the time-reversed trace (c). Blue and green ticks are theoretical P and S arrival times, respectively, and the grey tick at ∼8.1 s corresponds to 5 s after the empirical coda start time. We choose the waveform segments between the two grey marks for calculating goodness-of-fit measures and spectral amplitude estimates in section Amplitude Decay in the main text. Figure S4. Computational domain for synthetic tests with the stations (black triangles) and ∼2400 random noise source locations (grey stars) generated in 6 min over the earth's surface. Figure S5. Comparison of nine-component SNGF (black) and SGF (red) tensors at various interstation distances. The SNGFs were obtained from cross-correlation of synthetic noise. Each panel (3 × 3 group of subplots) is specific to an interstation distance (shown at the top of the panels). Each subplot in a panel is specific to a component pair (component pair names are at top right corner of each subplot; distance and azimuth are at the top left corner of the TZ subplots). In component pairs (i.e. TR), the first component (T) corresponds to the force direction and second component (R) corresponds to the displacement direction at the receiver. The velocity model used to compute the SGF tensor and the synthetic noise waveforms is indicated at the bottom left corner of the TZ subplot. For the overall tensor comparisons of a station pair, two goodness-of-fit estimates, VR and normalized zero-lag CC, are mentioned at the top left corner of the TR subplots. For all non-zero component pairs, the normalized CC along with the time lag between the SNGF and SGF waveforms is indicated at the bottom right corners of the subplots, if the CC exceeds 0.7. For reference, blue and green marks are theoretical P and S arrival times, respectively, and the grey mark corresponds to 5 s after the empirical coda start time. For a particular station pair, following normalization of amplitudes of the SNGF and SGF tensors, only relative amplitudes are meaningful. Figure S6. Same as Fig. 6 in the main text, but for surface-borehole sensor pairs DVB-DEB and SSR-SRB. Figure S7. (a) Same as Fig. 10(c) in the main text but for RR component SGFs at interstation distances of 14–17 km. (b) Same as (a) but for ZR component SGFs. Figure S8. Same as Fig. 6(b) in the main text, but for all distances and stations pairs with station DRH. Figure S9. (a) Sum of squares of amplitudes on T-[R,Z] components of NGFs normalized by the sum of squares of amplitudes on all components, as a function of interstation distance and interstation azimuth/backazimuth in [0°–180°] range. Coloured symbols represent data points for different distance ranges. Horizontal grey lines represent the azimuth range of the assumed dominant noise source direction [30°, 70°]. Vertical grey line marks φ = –1.0, the quality threshold used by Roux (2009). Colour coded numbers are the mean and standard deviation of φ in different distance ranges. AZ2 and AZ1 refer to NGFs of station pairs in azimuth/backazimuth range [30°–70°] and all other azimuths, respectively. Note that unlike fig. 5 in Roux (2009), this figure shows the data points plotted at actual interstation azimuth/backazimuth. (b) Same as (a) but with reduced φ following application of Optimal Rotation Algorithm. The data points are plotted at original interstation azimuth/backazimuth. (c) Same as Fig. 6(b) in the main text, but for two station pairs with very small amplitudes in the T-[R,Z] components despite being oriented (115° and 137°, respectively) at large angles with respect to the assumed dominant noise-source direction. Figure S10. Figure similar to Figs 16(c) and 16(d) in the main text but for TT component Fourier Spectral Amplitudes (FSAs) at 0.42 Hz instead of 0.72 Hz. (a) TT component NGF and SGF FSAs at 0.42 Hz as a function of interstation distance. The NGF FSAs (grey +) are uniformly scaled up by a constant such that the maximum value is 0. The frequency is indicated at the top left corner of the plots. Black diamonds and error bars represent the mean of NGF amplitudes and their standard deviation ± 1σ in bins for all bins with more than five data points. The bin widths are 1 and 2 km for distances ≤10 and >10 km, respectively. ‘SGF 3D’ and ‘SGF 1D’ are FSAs of equivalent SGFs of the 3-D and 1-D REF velocity models, respectively, plotted for the same station pairs. ‘No Q’ and ‘Q∼15’ indicate synthetics computed with weak and strong attenuation models, respectively. The SGF decay curves are scaled such that their amplitude at ∼1.5 km is similar to the average amplitude of NGF FSAs in the 1–2 km range. Scaling the binned NGF decay curve instead, to directly fit the REF model SGF decay curve by minimizing the L1 norm (black dashed line), doesn't change the absolute amplitudes significantly. The RSGF decay curve (red dashed line) is scaled to fit the REF model SGF decay curve at distances ≳ 8 km. (b) The data (green +) are tangential component FSAs of earthquake records at 0.42 Hz and the red diamonds and error bars represent the mean amplitudes and their standard deviation ± 1 σ in the bins. For comparison, we also plot binned TT component NGF FSAs at 0.42 Hz (black diamonds and error bars; same as [a]) that are scaled to minimize L1 norm of binned FSAs at distances ≳ 4 km. ‘N’ at the top right corner indicates the number of data points. Figure S11. (a) Tangential component velocity waveforms of an MD (duration magnitude) 2.68 earthquake at The Geysers on 2014–11-06, 09:54:57.50 (Northern California Seismic Network Event ID: 72336075). The waveforms are high-pass filtered at 0.1 Hz with a causal two-pole Butterworth filter. The name of the recording station and the peak amplitude in m s−1 are indicated near each waveform. For reference, the blue and green marks are the theoretical P- and S-wave arrival times, respectively. The waveform segments between the two grey marks were used for extracting spectral amplitudes (shown in panel b). There are large differences between the observed S-wave arrival times and the ones predicted by the REF velocity model at large distances. (b) Fourier amplitude spectra of the earthquake ground-motion velocity waveforms (between the grey marks in panel a) along with spectra of pre-event noise windows of the same duration. Black ‘+’ signs indicate the extracted earthquake FSAs at 0.25, 0.42 and 0.72 Hz used in this study. They pass the quality criteria, that is, they are greater than twice the noise FSAs (blue ‘+’ marks) for 12, 18 and 18 out of 18 stations at 0.25, 0.42 and 0.72 Hz, respectively. Figure S12. (a) RR component NGF and SGF/RSGF FSAs for a variety of models at 0.25 Hz (frequency indicated at the top) as a function of distance. Black diamonds and error bars represent the mean of NGF amplitudes and their standard deviation ± 1 σ in bins for all bins with more than 5 data points. The bin widths are 1 and 2 km for distances ≤10 and >10 km, respectively. 3-D and 1-D indicate synthetic FSAs for 3-D and 1-D velocity models, respectively. The synthetic decay curves are uniformly scaled by a constant to minimize the L1 norm of their difference with the binned NGF FSAs. (b) Same as (a) but for FSAs at 0.42 Hz. (c) The data (green +) are radial component FSAs of earthquake records at 0.25 Hz and the red diamonds and error bars represent the mean amplitudes and their standard deviation ± 1 σ in the bins. For comparison, we also plot binned RR component NGF FSAs at 0.25 Hz (black diamonds and error bars; same as panel a) that are scaled to minimize L1 norm of binned FSAs. ‘N’ at the top right corner indicates the number of data points. (d) Same as (c) but for FSAs at 0.42 Hz. Figure S13. RR component Green's function Fourier spectral amplitude |G11(ω)| (eq. 5) as a function of distance r and various terms in the solution at two different frequencies (indicated at the top right corner). The average decay of the near-field term and the far-field term are proportional to r− 2 and r− 1, respectively (indicated by the dashed green and red lines, respectively). The dashed gray box indicates the extent shown in Figs S12(a) and (b). Figure S14. TT component binned NGF FSAs at 0.72 Hz (black diamonds) and REF model SGF FSAs with various degrees of anelastic attenuation are same as in Fig. 16(c) in the main text. Different analytical formulations from Tsai (2011) plotted for distances r ≳ 3.6 km (∼1λ) are E2 = G(r) × (I0(αr))− 1 and $$E3\ = \ G( r ) \times {( {I_0^2( {\alpha r} ) - L_0^2( {\alpha r} )} )^{ - 0.5}}$$ for uniform and one-sided far-field noise source distributions, respectively, and $$E4\ = \ G( r ) \times {( {1 - \frac{{{r^2}}}{{4r{'^2}}}} )^{0.5}}$$ for truncated near-field noise source distribution. For comparison, we also plot simple E1 = G(r) and E5 = G(r)e− αr. Geometrical spreading G (r) = r− 0.5, α ∼ 0.04 km–1 at 0.72 Hz and r' ∼ 12 km. I0 is the modified Bessel function of first kind and order 0 and L0 is the modified Struve function of order 0. The analytical decay curves are scaled to have the same amplitude as the REF model SGF decay curve at ∼3.6 km. Table S1. List of earthquakes used to determine the horizontal orientation of sensors in this study. All but three earthquakes are located at teleseismic distances (>1000 km). FP is the frequency passband used: ‘1’ and ‘2’ are 0.02–0.05 and 0.05–0.1 Hz, respectively; N is the number of phases or waveform segments in the waveforms of each earthquake that were used for measuring alignment angles. Table S2. Mean azimuth of horizontal components of the sensors with respect to the geographic NS–EW reference frame (alignment angle), determined from the analysis of long period earthquake records and the associated uncertainty from multiple measurements. N is the number of measurements, that is, total number of alignment angles measured on separate waveform segments or phases; stations with just one measurement are in italic. The values of uncertainties are likely underestimated in the case of stations with too few measurements. We consider the results of other stations, not listed in this table, to be unreliable owing to large scatter in the measurements on different teleseismic records or to poor quality of the noise cross-correlations. However, the alignment angle at station SQK was verified in the field to be ∼14° (Ramsey Haught, private communication 2016). Table S3. Earthquakes at The Geysers used to study path attenuation of ground motions in section Amplitude decay in the main text. This set of earthquakes is adopted from the Northern California Seismic Network (NCSN) earthquake catalogue. The depths are relative to the Earth's geoid surface (approximately equal to mean sea level). Our reference datum is assumed to be ∼800 m above the mean sea level. MD = duration magnitude. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. APPENDIX A1 Noise cross-correlation methodology We download hourly time-series from the Northern California Earthquake Data Center (NCEDC), discarding any data with gaps in time. The data are demeaned, detrended, tapered, decimated to 20 Hz using an antialiasing filter and corrected for the corresponding instrument phase response to velocity. The waveforms are subsequently bandpass-filtered between 0.08 Hz and 4 Hz using a zero-phase 2-pole Butterworth filter. While large magnitude earthquakes at teleseismic distances are poorly recorded by the short-period sensors at The Geysers, high-rate microseismicity in the geothermal area, characterized by short duration (∼3–5 s), high frequency (≳1.0 Hz) signals, is a major source of non-stationary signals contaminating the ambient noise wavefield. We suppress the local earthquake signals present in the seismic records by applying a non-iterative water-level normalization (Bensen et al.2007) with a threshold of five times the daily mean of absolute values. To each sample exceeding this threshold, we apply an inverse weight equal to 25 times the ratio of mean amplitude measured over a 5 s envelope around that sample. Following Seats et al. (2012), we select 75 per cent overlapping, 150 s duration time windows for computing interstation coherency, applying amplitude response of an eight-pole Butterworth bandpass filter with corners at 0.2 Hz and 1.6 Hz as a band-limiting taper during spectral normalization (e.g. see fig. 7b in Bensen et al.2007). Coherency has been shown to perform better at retrieving phase information than the impulse response function in the presence of highly variable and additive random noise (Nakata et al.2011) and, therefore, is expected to provide more accurate phase information at an industrial site such as The Geysers. We preserve the relative amplitudes between different components of the three-component sensors by using the maximum of the temporal normalization weights of the three individual components and the mean of their spectral amplitudes for temporal and spectral normalization of all components, respectively (Lin et al.2014). The cross-correlations for all windows in one day are stacked to form a daily average and the averages are stacked for all available days to form a final reference stack that we use as our estimate of NGF. Cross-correlations are done in an EW-NS-Z reference frame and are subsequently rotated to a T-R-Z reference frame relative to the source-receiver azimuths of the station pairs (Lin et al.2014). A2 Quality control of NGFs For quality control, we identify time periods when the cross-correlations of various stations with the broad-band reference stations HOPS and MNRC return little or no coherent energy by cross-correlating daily stacks with reference stacks for the overall time period (correlation threshold ≲ 0.3). These daily stacks are subsequently removed from the reference stacks after visual confirmation of all poor coherence time periods that lasted ≳ 10 d. Through this process, ∼45 per cent of all daily noise cross-correlations were removed. While the time period of our study is March 2012–August 2015, the actual duration of usable data for different station pairs varies from ∼25 days to the entire time period. The daily stacks are also corrected for large and obvious (>0.05 s) clock or time stamp errors that were identified by cross-correlating daily stacks with the reference stack (e.g. Sens-Schönfelder 2008). Station pairs with little or no coherent energy in the ballistic wave arrival time window (6.0 to 1.5 km s) in the reference stacks are removed from further analysis. Most pairs involving the 4.5 Hz geophone stations BRP, EPR, JKR, RGP, SB4 and SQK were discarded during the quality control process indicating high sensor self-noise, poor coupling of the sensors to the ground, lack of coherent ambient noise at these sites, or other sensor/digitizer problems. We also observe strong narrow-band spikes in amplitude spectra of cross-correlations at specific frequencies, namely 1.00, 1.76, 1.80, 1.97, 2.28 and 2.57 Hz, similar to the ones observed at high frequencies in other studies (e.g. Wegler et al.2009; Takagi et al.2015). These peaks were also found in spectra of raw data and we speculate that they result from coherent data logger noise, with the peak at 1.00 Hz most likely resulting from GPS time calibration (Takagi et al.2015). We manually reduce the amplitudes of noise spectra by a factor of ∼10−3 in narrow frequency bands of width ∼0.02 to ∼0.07 Hz around these frequencies during spectral whitening. Finally, we are left with NGFs of ∼568 station pairs not including station HOPS or MNRC. A3 Processing of earthquake waveforms We selected ∼121 earthquakes between MD ∼ 2.5–2.7 and depths <4 km from the NCSN earthquake catalogue available at NCEDC (Supporting Information Table S3). While earthquakes with greater magnitudes provide records with better signal-to-noise ratio, we find that the data recorded by the 4.5 Hz geophones are contaminated with long-period transients synchronous with the onset of strong shaking at many stations. For this analysis, we also use the co-located three-component accelerometers at the USGS stations GAXB and GDXB and remove any record with maximum absolute value >106 counts. The EW and NS component data were decimated to 20 Hz, corrected for instrument response to represent velocity, visually checked for quality and rotated to radial and tangential directions. Waveform segments of the same distance-dependent durations were selected and FSAs were extracted using the methodology similar to that adopted for the RSGFs. For the earthquakes, we also extract FSAs of similar duration noise-windows prior to the origin time and use only data, for which signal-to-noise ratio of FSAs is >2. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

Empirical Green's tensor retrieved from ambient noise cross-correlations at The Geysers geothermal field, Northern California

Loading next page...
 
/lp/ou_press/empirical-green-s-tensor-retrieved-from-ambient-noise-cross-HrzVZszfqE
Publisher
The Royal Astronomical Society
Copyright
© The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx534
Publisher site
See Article on Publisher Site

Abstract

Summary We retrieve empirical Green's functions in the frequency range (∼0.2–0.9 Hz) for interstation distances ranging from ∼1 to ∼30 km (∼0.22 to ∼6.5 times the wavelength) at The Geysers geothermal field, Northern California, from coherency of ambient seismic noise being recorded by a variety of sensors (broad-band, short-period surface and borehole sensors, and one accelerometer). The applied methodology preserves the intercomponent relative amplitudes of the nine-component Green's tensor that allows us to directly compare noise-derived Green's functions (NGFs) with normalized displacement waveforms of complete single-force synthetic Green's functions (SGFs) computed with various 1-D and 3-D velocity models using the frequency-wavenumber integration method and a 3-D finite-difference wave propagation method, respectively. These comparisons provide an effective means of evaluating the suitability of different velocity models to different regions of The Geysers, and assessing the quality of the sensors and the NGFs. In the T-Tangential, R-Radial, Z-Vertical reference frame, the TT, RR, RZ, ZR and ZZ components (first component: force direction, second component: response direction) of NGFs show clear surface waves and even body-wave phases for many station pairs. They are also broadly consistent in phase and intercomponent relative amplitudes with SGFs for the known local seismic velocity structure that was derived primarily from body-wave traveltime tomography, even at interstation distances less than one wavelength. We also find anomalous large amplitudes in TR, TZ, RT and ZT components of NGFs at small interstation distances (≲4 km) that can be attributed to ∼10°–30° sensor misalignments at many stations inferred from analysis of longer period teleseismic waveforms. After correcting for sensor misalignments, significant residual amplitudes in these components for some longer interstation distance (≳8 km) paths are better reproduced by the 3-D velocity model than by the 1-D models incorporating known values and fast axis directions of crack-induced VS anisotropy in the geothermal field. We also analyse the decay of Fourier spectral amplitudes of the TT component of NGFs at 0.72 Hz with distance in terms of geometrical spreading and attenuation. While there is considerable scatter in the NGF amplitudes, we find the average decay to be consistent with the decay expected from SGF amplitudes and with the decay of tangential component local-earthquake ground-motion amplitudes with distance at the same frequency. Seismic anisotropy, Seismic interferometry, Seismic noise, Surface waves and free oscillations, Wave propagation, Crustal structure INTRODUCTION The Geysers in the Mayacamas Mountains of Northern California are the largest complex of geothermal power plants in the world today. Spread over an area of ∼115 km2, it currently produces ∼835 MW of electricity (http://www.energy.ca.gov/tour/geysers/, last accessed August 2016). It is a vapour-dominated geothermal reservoir (Allis & Shook 1999) in which approximately 20 million gallons of reclaimed wastewater are injected on a daily basis (http://www.geysers.com/water.aspx, last accessed August 2016) from neighbouring communities (Lake County, since 1998; City of Santa Rosa, since 2004). The water is injected to tap the heat in the volcanic reservoir rocks that are at ∼240 to ∼300° C (Truesdell et al.1991; Lowenstern & Janik 2003; Majer & Peterson 2007). While the injection of wastewater has successfully replenished the resource and the steam production that was in severe decline during the late 1980s, it has also been linked to elevated rates of microseismicity (MW ≤ 3) in the reservoir when compared to regional rates (Ludwin et al.1982). Many studies have demonstrated strong temporal correlation of microseismicity with both steam production and water injection at The Geysers, on both localized and field-wide scales (Oppenheimer 1986; Stark 1991; Majer & Peterson 2007; Johnson et al.2016; Trugman et al.2016). Detailed analysis of microseismicity has also revealed significant temporal variations in earthquake hypocenter distributions, in b-value and in the number of large magnitude (MW > 4) earthquakes (Trugman et al.2016). Investigations of focal mechanisms and moment tensors of these earthquakes using body-wave polarities, amplitude ratios, and low-frequency displacement waveforms have shed some light on various aspects of their source mechanisms (Oppenheimer 1986; Ross et al.1999; Guilhem et al.2014; Boyd et al.2015). However, imperfect knowledge of subsurface seismic velocity structure and assumption of 1-D velocity models are believed to introduce considerable uncertainties in synthetic ray paths (Oppenheimer 1986) and in synthetic Green's functions (SGFs) required for source inversion. 3-D P and S-wave traveltime tomography studies have revealed considerable heterogeneities in the subsurface (Julian et al.1996; Ross et al.1999; Gritto et al.2013a) that must be taken into account in the inversion of polarities and seismic waveforms for recovery of accurate source mechanisms. In addition to the temporal changes in microseismicity, temporal changes in bulk medium properties of the reservoir such as fluid pressure and VP/VS ratio have also been observed from levelling and GPS studies, time-lapse body-wave travel time tomography, etc. (Foulger et al.1997; Mossop & Segall 1999; Gritto et al.2013b). The response of a thermal reservoir to water injection or production, generally indicated by (1) the spatiotemporal behaviour and mechanisms of seismicity, and/or (2) changes in seismic properties of the subsurface, is of great interest to further our understanding of induced seismicity and geothermal systems. This knowledge can help to better characterize and manage the seismic hazard by optimizing fluid injection into the subsurface for energy or storage purposes. In this study, we establish a framework for ambient seismic noise cross-correlation at The Geysers, which will aid in improving our understanding of the seismic velocity structure. The cross-correlations of ambient seismic noise obtained at pairs of seismic sensors have been shown to converge over a period of time to the empirical Green's functions (GFs) or the seismic medium's response at one receiver to a unit force applied at the other receiver (Shapiro & Campillo 2004; Shapiro et al.2005; Bensen et al.2007). In numerous studies, the body and surface waves in the noise-derived Green's functions (NGFs) have been used for tomography at various scales (e.g. Shapiro et al.2005; Lin et al.2008, 2014; Nakata et al.2015). Additionally, the coda of NGFs can be used for spatiotemporal monitoring of small subsurface seismic velocity changes (e.g. Sens-Schönfelder & Wegler 2006; Brenguier et al.2014; Taira et al.2015) that have previously been observed at The Geysers with body-wave tomography. However, the non-uniformity in ambient noise source distribution results in systematic differences between recovered NGFs and the exact GFs in the real Earth (Tsai 2009). This can lead to errors in the estimates of subsurface seismic properties inferred from NGFs and limits their utility and reliability. In this study, we retrieve NGFs, represented by interstation coherency (Prieto et al.2011), in the frequency range of ∼0.2–0.9 Hz for interstation distances from ∼1 to ∼30 km at the Geysers. We use seismic data that are continuously recorded by a wide variety of sensors in and around the geothermal reservoir. The minimum VS is ∼2 km s−1 and the interstation distances range from ∼0.22 λ (near-field) to ∼6.5 λ (far-field) assuming a phase velocity ∼2.3 km s−1 at a frequency of ∼0.5 Hz, where λ is the frequency-dependent surface-wave wavelength and distances ≳1λ–2λ are assumed to be far-field. For pairs of three-component sensors (T: Tangential, R: Radial, Z: Vertical), we retrieve nine components of NGFs (TT, TR, etc. where the first component is direction of the applied force and the 2nd component is the displacement direction). We directly compare NGF waveforms to normalized single-force displacement SGFs, computed using various 1-D velocity models utilized for various monitoring applications at The Geysers. We evaluate the similarity of NGF and SGF waveforms in terms of waveform fits, phase and relative intercomponent amplitudes. NGFs with little or no coherent energy in the ballistic wave arrival time window or NGFs that bear little or no resemblance to SGFs expected for a wide variety of reasonable and expected velocity models are interpreted to be contaminated with errors possibly due to poor sensor coupling, or non-uniformity of ambient noise source distribution. After discarding these erroneous NGFs, we use the remaining NGFs to evaluate the applicability of various 1-D velocity models to different sub-regions of The Geysers as NGFs contain information about real Earth 3-D wave propagation (Ma et al.2008). We examine if NGFs for station pairs in a particular sub-region are systematically better fit by SGFs of a particular velocity model compared to SGFs of other velocity models. We interpret the phase differences between NGFs and SGFs in terms of the deficiencies in the velocity models. While some phase errors in the NGFs are expected due to the non-uniformity of ambient noise source distribution, our interstation paths are well distributed in azimuth, which allows us to draw inferences from spatially coherent phase differences pervasive over multiple azimuths. We show that waveforms, phases and intercomponent relative amplitudes of the primary non-zero components of retrieved NGFs (i.e. TT, RR, RZ, ZR, ZZ) are similar to those of SGFs, even at distances < λ. We compare our inferences on the features of subsurface seismic velocities at The Geysers from the NGF-SGF comparisons with the inferences drawn in many previous body-wave traveltime studies. We also examine if the NGFs compare better to SGFs computed using a 3-D velocity model of the geothermal field than with SGFs from the 1-D velocity models. Our study additionally includes: (1) an analysis of long period teleseismic waveforms to verify sensor orientations, (2) application of the Optimal Rotation Algorithm (Roux 2009) on the NGFs to detect any dominant noise source illumination direction, and (3) comparison with SGFs computed using the 1-D models incorporating known values and fast axis directions of crack-induced VS anisotropy in the geothermal field. We evaluate the contribution of the above-mentioned factors to the significant non-zero amplitudes in TR, TZ, RT and ZT components of NGFs for many stations pairs. We also compare the decay of Fourier spectral amplitudes of the TT component of NGFs at 0.72 Hz as a function of distance with the decay of synthetic spectral amplitudes for velocity models incorporating strong, weak and known values of anelastic attenuation, and with the decay of tangential component local-earthquake ground-motion amplitudes. BACKGROUND AND MOTIVATION Noise cross-correlation tomography studies typically involve phase velocity measurements on the fundamental-mode Rayleigh waves obtained from cross-correlation of ambient noise recorded on vertical components of broad-band stations at far-field interstation distances (e.g. Yao et al.2006; Lin et al.2008, 2014). The theoretical proof for retrieval of ballistic surface waves from cross-correlation of diffused waves at two receivers was first provided by Snieder (2004) under a stationary phase approximation valid at far-field distances (Yokoi & Margaryan 2008; Zhan & Ni 2010). Subsequently, using source-receiver reciprocity, Wapenaar & Fokkema (2006) derived an expression for the velocity GF response at a receiver to a single force applied at another receiver in an inhomogeneous anisotropic loss-less medium from cross-correlation of responses of sources distributed over an enclosing surface. The variations of their basic formulation involve variations in the nature of the boundary and the outside medium, types of sources, types of responses cross-correlated or retrieved, distance of the enclosing source surface from the receivers and whether the cross-correlation operates on the responses of individual sources or the net ground motion field at the receivers. In practice, the application of most temporal and spectral normalization methods removes the absolute amplitude information in component-pair NGFs leading to their ‘empirical’ nature (Bensen et al.2007). However, significant relative amplitude information in addition to the phase information can still be retrieved from NGFs, although many aspects are still under considerable debate as discussed in the Amplitude decay section. Cross-correlation of spectrally whitened noise is equivalent to computing interstation coherency (Bensen et al.2007; Prieto et al.2011). For components i, j of velocity v recorded at stations A, B at positions xA, xB respectively, coherency γij(xA, xB, ω) as a function of frequency ω is defined as   \begin{eqnarray} - {G_{ij}}\left( {{x_A},{x_B},\omega } \right)& \approx & {\gamma _{ij}}\ \left( {{x_A},{x_B},\omega } \right) \nonumber\\ &=& \left\langle {\frac{{v_i^*\left( {{x_A},\omega } \right)\ {v_j}\left( {{x_B},\omega } \right)}}{{\left\{ {\left| {{v_i}\left( {{x_A},\omega } \right)} \right|} \right\}\left\{ {\left| {{v_j}\left( {{x_B},\omega } \right)} \right|} \right\}}}} \right\rangle , \end{eqnarray} (1) where Gij(xA, xB, ω) is the displacement NGF with single force, displacement directions i, j at stations A, B that are acting as source and receiver, respectively; <> implies stacking results for data recorded over multiple time windows known as ensemble averaging, {} implies spectral amplitude smoothing (three-point ∼0.02 Hz moving window average in our study). The directions i, j can be radial (R; positive outwards), transverse (T; positive clockwise 90° from radial) or vertical (Z; up positive). In eq. (1), in the following and throughout this text, we define single-force GF Gij(xA, xB, ω) as the displacement response to an input step force and directly relate it to interstation coherency obtained from the cross-correlation of spectrally whitened ambient noise. This is different than the traditional definition of a GF in which it is defined as the displacement response to an input impulsive force (Aki & Richards 2002). For the traditional definition (i.e. impulsive force as an input), the GF can be related to the time derivative of coherency as in many other studies, that is, Lin et al. (2008), Ma et al. (2008), etc. For a particular component pair, γ has been shown to preserve geometrical spreading and anelastic attenuation information as a function of interstation distance and structure (Prieto et al.2011; Lawrence et al.2013). Additionally, if the same normalization factors are used for all components of a station, relative intercomponent amplitudes for a particular station pair can be preserved enabling extraction of attributes such as the Rayleigh wave ellipticity from the amplitude ratio of R and Z components of NGFs (Lin et al.2014). The amplitude spectrum of the noise at the source station can also be used to normalize the noise spectra at both stations, which yields the impulse response function (Prieto et al.2011; Denolle et al. 2013), and is known as deconvolution interferometry (Nakata et al.2011). The impulse response function has been used to extract site response of basins and 3-D structures at long periods, ∼4–10 s (Prieto & Beroza 2008; Denolle et al.2014). The retrieved phase information is the same for different spectral normalization methods (Prieto et al.2011). While interstation distances ≳λ are considered to be ‘far-field’ for surface waves (Lin et al.2013), phase velocity measurements from NGFs are usually restricted to interstation distances ≳3λ in order to avoid bias at shorter distances caused by inhomogeneous noise source distributions (Lin et al.2008, 2014). Various numerical studies have shown that phase velocity measurements derived from NGFs at interstation distances ≲3λ–5λ can have significant error, even for seismic noise sources that are azimuthally widely distributed (Kimman & Trampert 2010; Zhan & Ni 2010). However, reliable phase velocities can be extracted from NGFs at distances down to ∼1λ using sophisticated methods that are able to reduce biases caused by inhomogeneously distributed noise sources. These methods include measurement of phase velocities by phase-front tracking on dense arrays (Lin et al.2013) and azimuthal averaging of multiple coherency measurements as in spatial autocorrelation (Tsai & Moschetti 2010). The studies on relative amplitude extraction from NGFs, mentioned earlier, have also utilized NGFs at primarily far-field distances assuming Rayleigh waves on RZ, RR, ZR and ZZ components and Love waves on TT components, which are the non-zero components of the GF tensor for 1-D isotropic media. The TR, TZ, RT and ZT components of the GF tensor can be non-zero in the presence significant 3-D structure or anisotropy. In NGFs, these components can also be non-zero due to inhomogeneous noise source distribution [e.g. Parkfield, California (Durand et al.2011)] and have not been analysed in detail in crustal structure studies [e.g. some paths in Southern California (Denolle et al.2013)]. Sensor misalignments or incorrect knowledge of sensor gains and orientations may lead to incorrectly oriented NGFs and anomalous amplitudes in the TR, TZ, RT and ZT components as well. To the best of our knowledge, the retrieval of relative amplitude information in NGFs at distances < λ has not yet been explored. METHODOLOGY Data Fig. 1 shows a map of The Geysers, the seismic stations employed in this study and the approximate outline of the steam field of the reservoir area (Gritto et al.2013a). Most of our data are recorded by the Lawrence Berkeley National Laboratory (LBNL) network at The Geysers consisting of ∼30 short-period (4.5 Hz) three-component geophones (OYO-GS11D), one Nanometrics Titan accelerometer (station DRH) and three short-period (8.0 Hz) three-component borehole sensors (OYO-GS11D8) at depths of ∼35–150 m (stations DEB, JKB and SRB). We include data recorded by short-period (1.0 Hz) vertical component sensors (L4) and by one broad-band STS2 sensor (station GDXB) at stations operated by the United States Geological Survey (USGS) in and around the geothermal field. Data recorded by the two closest broad-band stations (an STS1 at HOPS and an STS2 at MNRC) of the Berkeley Digital Seismic Network (BDSN) are also included to provide a reference for quality control of data recorded at all other stations in our study. Figure 1. View largeDownload slide Map shows the outline of the reservoir area (red polygon), faults (solid black lines) and seismic stations plotted on the gradient of topography at The Geysers. Black triangles: LBNL 4.5 Hz geophones, green squares: LBNL 8.0 Hz borehole geophones, red square: LBNL accelerometer, blue triangles: USGS 1.0 Hz vertical component sensors, and blue square: USGS broad-band sensor. Inset shows the location of The Geysers (red square) and BDSN stations (red triangles) with respect to the San Francisco Bay area in Northern California, USA. Figure 1. View largeDownload slide Map shows the outline of the reservoir area (red polygon), faults (solid black lines) and seismic stations plotted on the gradient of topography at The Geysers. Black triangles: LBNL 4.5 Hz geophones, green squares: LBNL 8.0 Hz borehole geophones, red square: LBNL accelerometer, blue triangles: USGS 1.0 Hz vertical component sensors, and blue square: USGS broad-band sensor. Inset shows the location of The Geysers (red square) and BDSN stations (red triangles) with respect to the San Francisco Bay area in Northern California, USA. Cross-correlation analysis Our data pre-processing and cross-correlation workflows are derived from Bensen et al. (2007), Seats et al. (2012) and Lin et al. (2014). Most importantly, we use the same temporal and spectral normalization factors for all components of a three -component station in order to preserve the relative amplitudes between the components (e.g. Lin et al.2014). Other details of the methodology adopted are provided in  Appendix A1. The details of the subsequent quality control procedure are described in  Appendix A2. After quality control, we average the causal and anticausal sides of the final stacked cross-correlations or coherencies, extracting the symmetric component of NGFs in the T-R-Z reference frame. Fig. 2 shows ZZ and TT components of NGFs in the interstation distance range 0–25 km in which most of our station pairs are located. The waveforms are stacked in 0.5 km interstation distance bins. Rayleigh and Love waves can easily be identified in components ZZ and TT, respectively, with apparent velocities of ∼2.5 and ∼2.8 km s−1, respectively. ZZ components at distances of ∼17–20 km also show faint but coherent body-wave or higher-mode Rayleigh wave energy that is faster than the fundamental-mode Rayleigh waves but slower than direct P-waves, traveling at ∼3.6 km s−1 (composition of waveforms further discussed in the Results section). We also estimate an empirical expression for distance-dependent approximate start time of coda waves following the passage of the surface waves by fitting a fourth degree polynomial to times at which envelopes of NGFs drop below a level of 0.2 times the peak absolute amplitude. This empirical coda-wave start time is used for selecting waveform segments that contain body-wave and surface-wave phases in order to evaluate waveform fits between NGFs and SGFs. Figure 2. View largeDownload slide NGFs filtered between ∼0.2–0.9 Hz and stacked in 0.5 km interstation distance bins. Dashed red curves are reference traveltime curves for various constant apparent velocities (shifted to the left by 0.7 s for clarity). Dashed blue curve shows the empirically estimated approximate start time of coda waves after the passage of surface waves. The NGF waveforms are proportional to displacement response for an input step force. Figure 2. View largeDownload slide NGFs filtered between ∼0.2–0.9 Hz and stacked in 0.5 km interstation distance bins. Dashed red curves are reference traveltime curves for various constant apparent velocities (shifted to the left by 0.7 s for clarity). Dashed blue curve shows the empirically estimated approximate start time of coda waves after the passage of surface waves. The NGF waveforms are proportional to displacement response for an input step force. Velocity models Here we briefly describe the velocity structure of the reservoir area. Gritto et al. (2013a) conducted a tomographic campaign to jointly invert for the 3-D P- and S-wave velocity structures and the hypocentre locations of microseismicity at The Geysers using body-wave arrival times of ∼32,000 earthquakes at the LBNL stations. The starting 3-D velocity model was based on 3-D interfaces in the reservoir (interface of top steam entry and interface of felsite) coupled with P- and S-wave velocity estimates that were derived from a vertical seismic profiling (VSP) experiment in the centre of the reservoir (O'Connell & Johnson 1991). The top of the steam interface is considered the top of the reservoir and separates mélange (a complex assemblage of Franciscan greywacke, greenstone and serpentinite) from the deeper metamorphosed reservoir rocks (Franciscan greywacke and metagraywacke) at depths of ∼1–3 km. The top of the felsite interface demarks the transition from metamorphosed reservoir rock to the underlying granitic pluton at depths greater than ∼3–5 km. On average, VP, VS values increase, respectively, from ∼3.6 and ∼2 km s−1 at the surface, to ∼5.9 and ∼3.3 km s−1 at depth of ∼6 km, with typical values of ∼4.8 and ∼2.8 km s−1 for reservoir rocks, and ∼5.5 and ∼3.0 km s−1 for the underlying felsite (O'Connell & Johnson 1991). The reservoir area is bounded by NNW–SSE trending faults, namely the Mercuryville fault to the southwest and the Collayomi fault to the northeast. There is widespread lateral heterogeneity in and around the reservoir, with up to ∼20 per cent variations in VP at shallow (less than ∼2.5 km) depths (Julian et al.1996; Ross et al.1999). On average, the velocities in the northwestern section of The Geysers have been found to be lower by ∼10 per cent than in the central and southeastern section (Eberhart-Phillips 1986; Julian et al.1996; Gritto et al.2013a). The VP/VS ratio, believed to be relatively insensitive to lithology but quite sensitive to the saturation of rocks (Gritto & Jarpe 2014) and to the compressibility of pore fluids, shows a large negative anomaly (up to ∼–9 per cent) at the centre and southeast of the centre of the reservoir (Julian et al.1996; Foulger et al.1997; Ross et al.1999). This anomaly exists at depths up to ∼2.5 km and is anticorrelated with a positive VS anomaly (Gritto et al.2013a). We select four 1-D velocity models (Fig. 3) at The Geysers to compute single force SGFs for comparison with the NGFs of all station pairs not including stations HOPS and MNRC. The models are, (1) REF, (2) BACK1, (3) VSP0, and (4) AVG1. Model REF is a slightly modified version of the 1-D average of the 3-D model of Julian et al. (1996) for the NW Geysers and was used by Guilhem et al. (2014) for modelling 0.5–2.5 Hz waveforms of MW ∼ 3.3–3.9 earthquakes at distances up to a few km. Model BACK1 is a 1-D depth average of the 3-D starting velocity model used in Gritto et al. (2013a), while model AVG1 is a 1-D depth average of their final 3-D velocity model within the reservoir area. VSP0 is the velocity profile shown in fig. 2 in Gritto et al. (2013a) and is derived from a VSP profile at the felsite peak in the SE Geysers (O'Connell & Johnson 1991). The depth-dependent density and anelastic attenuation quality factors (QP and QS) are borrowed from Guilhem et al. (2014). In our study, we ignore the surface topography at The Geysers and assume our reference datum to be ∼800 m above the mean sea level, which is the average elevation of stations that we are using. For all station pairs that include station HOPS or MNRC, we use the California Central Coast Ranges velocity model GIL7 (Stidham et al.1999) and its modified version GIL7.1 (Fig. 3) to compute SGFs (modifications explained in section NGFs with BDSN stations HOPS and MNRC). Upper crustal QP and QS values in GIL7 and its modified versions are reduced to values comparable to those of REF. Figure 3. View largeDownload slide 1-D velocity models used in our study. Figure 3. View largeDownload slide 1-D velocity models used in our study. Fig. 4 shows a few example depth slices of G3D1, a modified version of the 3-D velocity model of Gritto et al. (2013a) used to compute SGFs accounting for the 3-D structure in the reservoir in our study. G3D1 consists of a 3-D velocity structure surrounded by a 1-D background model (BACK1). The horizontal extent of the 3-D structure is approximately restricted by the outline of the reservoir steam field at the surface (red polygon in Fig. 1) (except towards the northwest where the LBNL network extends beyond the outline; see stations HER, AL2 and DRH in Fig. 1), and decreases with increasing depth. The modifications to the model are explained as follows. For nodes with VP/VS ratio $$ < \sqrt 2 $$, VP and VS values are increased and decreased respectively in small increments, until they satisfy the minimum VP/VS ratio criteria of $$\sqrt 2 $$ of our 3-D wave propagation software. The velocities are tapered at the edges of the model so that they smoothly transition to model BACK1. The model is extended laterally assuming the values of BACK1, and extended in depth assuming a half-space with VP ∼ 5.86 km and VS ∼ 3.35 km at depths ≳ 5.3 km, which are comparable to the deepest velocities in BACK1 and velocities in model GIL7 at the same depths. Figure 4. View largeDownload slide Depth slices of model G3D1 used to compute SGFs. Dashed line is the outline of the geothermal steam field and black triangles are locations of stations used in this study (same as in Fig. 5). Outside this domain, we extrapolate the velocity model to match the 1-D model BACK1. Figure 4. View largeDownload slide Depth slices of model G3D1 used to compute SGFs. Dashed line is the outline of the geothermal steam field and black triangles are locations of stations used in this study (same as in Fig. 5). Outside this domain, we extrapolate the velocity model to match the 1-D model BACK1. Synthetic Green's functions For all 1-D models, single force SGFs were computed using the frequency–wavenumber integration (FKI) method based on Haskell (1964) and Wang & Herrmann (1980), as provided in Herrmann (2013a). The Herrmann (2013a) FKI code computes complete three-component seismograms consisting of all the terms of the elastic wave equation solution (including the near- and intermediate-field terms) and all body-wave and surface-wave phases for isotropic 1-D layered velocity models with anelastic attenuation. We use 10 m as source depth for all stations at the surface (the numerical integration method used requires a non-zero depth difference between the source and the receiver) and depths reported in Northern California Earthquake Data Center (NCEDC) as source or receiver depths for borehole stations. For computation of all SGFs (both 1-D and 3-D) in this study, we use impulse or Dirac source-time functions for numerical stability. The ground motions thus obtained are displacement waveforms for an impulse input or effectively velocity waveforms for a step input. The velocity waveforms are subsequently integrated to displacement SGFs (for an input step force). We use model G3D1 with the 3-D seismic wave propagation code SW4 to compute SGFs accounting for the 3-D structure in the reservoir. SW4 is a 3-D finite difference seismic wave propagation code that solves the elastic wave equation in displacement formulation with fourth order accuracy in space and time (Sjögreen & Petersson 2012; Petersson & Sjögreen 2014a,b). Absorbing supergrid boundary conditions are used at the far-field boundaries and anelastic attenuation is also allowed (Petersson & Sjögreen 2014b). Our grid spacing is 80 m, which allows ∼27 grid points per shortest wavelength for the minimum VS ∼ 2.0 km s−1 in our model and the maximum frequency of our interest ∼0.9 Hz. For the maximum VP ∼ 6.4 km s−1 in the model, the time step is set to ∼1.18 × 10−2 s for numerical stability. Our model dimensions are ∼36.6 km (NS) x ∼41.2 km (EW) and ∼14.0 km (Z), ∼42 million grid points, with supergrid absorbing layer thickness of 75 grid points. Fig. 5 shows a snapshot of an example wavefield computed using SW4 in response to a single force applied at station DRH. The displacement time histories for this wavefield can be compared to NGFs with DRH as the ‘source’ station. The effect of faster than background velocities in the central section of The Geysers is reflected in the bulge in the otherwise circular ∼0.9 Hz wave front along with reduced amplitudes. Note the absence of any significant reflections from the boundaries of the domain showing the effectiveness of the supergrid absorbing boundary conditions. Figure 5. View largeDownload slide The left subplot shows snapshot of the vertical component velocity wavefield at the surface at elapsed time 9.86 s in response to a vertical force (1e + 10 N) applied at location of station DRH (red star). The amplitudes of the wavefield are normalized by the instantaneous peak absolute amplitude (indicated at the bottom right corner) and plotted on a red (-1) to blue ( + 1) colour scale. The domain shown here corresponds to the computational domain used in SW4 and the dotted line marks the supergrid absorbing layer boundary. Explanation of other symbols on the map is the same as in Fig. 4. The right subplot shows the time history of the same wavefield at station GDXB (red square on the map). The red tick is the instant corresponding to the snapshot shown on the map. A Gaussian source time function with fundamental frequency ∼0.9 Hz (freq parameter in SW4 ∼ 5.65) was used for this particular wavefield. Figure 5. View largeDownload slide The left subplot shows snapshot of the vertical component velocity wavefield at the surface at elapsed time 9.86 s in response to a vertical force (1e + 10 N) applied at location of station DRH (red star). The amplitudes of the wavefield are normalized by the instantaneous peak absolute amplitude (indicated at the bottom right corner) and plotted on a red (-1) to blue ( + 1) colour scale. The domain shown here corresponds to the computational domain used in SW4 and the dotted line marks the supergrid absorbing layer boundary. Explanation of other symbols on the map is the same as in Fig. 4. The right subplot shows the time history of the same wavefield at station GDXB (red square on the map). The red tick is the instant corresponding to the snapshot shown on the map. A Gaussian source time function with fundamental frequency ∼0.9 Hz (freq parameter in SW4 ∼ 5.65) was used for this particular wavefield. Sensor orientations Rotation of the NGF tensor to the T-R-Z reference frame with respect to the source-receiver azimuth requires precise knowledge of sensor orientations. Among the stations employed in our study, the exact orientation of the borehole sensors is unknown (http://www.ncedc.org/egs/, last accessed July 2016). Moreover, preliminary comparisons of NGFs and SGFs (described in the following section) indicated possible sensor misalignments, that is, non-zero angles between the horizontal components of the three-component sensors and the geographic NS–EW reference frame. We compare three-component long period waveforms (∼0.02–0.1 Hz) of regional and teleseismic earthquakes at all stations to the corresponding waveforms at the broad-band USGS station GDXB to verify sensor orientations and determine the degree of misalignment, if any (e.g. Grigoli et al.2012). The details of the methodology and results are presented in the Supporting Information (section Sensor orientations, Figs S1 and S2, and Tables S1 and S2). Our results provide the first orientation estimates for the borehole sensors and show surprisingly large alignment angles (>10°) at many of the surface stations as well (stations ACR, AL6, CLV, DRK, DVB, FUM, HVC, LCK, MNS, NEG and TCH) with standard deviations < 5° over multiple measurements. Procedure for comparison between NGFs and SGFs We find that proper filtering of SGFs is a critical step for comparison with NGFs at small interstation distances. As a result of the band-limiting taper applied during spectral whitening and cross-correlation (see Appendix A1), the NGFs are acausal. For positive correlation lag time, NGFs at short interstation distances will have non-zero acausal amplitudes leaking to t < 0 following zero-phase bandpass filtering that contaminate the symmetric time-reversed NGFs at negative correlation time and vice versa for t > 0 (e.g. Cho et al.2007). Therefore, we zero pad the raw SGFs at t < 0, bandpass-filter them between 0.2 and 1.6 Hz with a zero-phase eight-pole Butterworth filter (same as the taper applied during spectral normalization), time-reverse the SGFs and take the symmetric component. This is explained in detail in Supporting Information Fig. S3. All NGFs and SGFs are subsequently bandpass-filtered between 0.2 Hz and 0.9 Hz with a causal two-pole Butterworth filter. For comparison of the NGF and the SGF tensors, we choose waveform segments between 1 s prior to the theoretical P-wave arrival time or t = 0 (whichever is later) and 5 s after the empirical coda-wave start time. The tensors are normalized by the overall L2 norm of all available components, which preserves relative intercomponent amplitudes and makes the scaling factor invariant under the rotation of the station components. However, this also creates a slight bias against the SGFs of 1-D velocity models, because TR, TZ, RT and ZT components of SGFs are identically zero, whereas corresponding components in NGFs are always non-zero. This is manifested by slightly greater amplitudes for the five primary, non-zero components- TT, RR, RZ, ZR and ZZ of the SGF tensors compared to the same for the NGF tensors for many station pairs, especially at greater distances at which the amplitudes in the other NGF components are dominated by noise in the NGF time-series (see section Contribution of 3-D structure and anisotropy). The overall similarity between the two tensors is evaluated by two estimates of goodness-of-fit, namely by the variance reduction (VR; e.g. Guilhem et al.2014), and by the zero-lag normalized correlation coefficient (CC). P- and S-wave arrival times are computed using velocity model REF for all station pairs excluding stations MNRC and HOPS, and models GIL7 or GIL7.1 for all station pairs including stations MNRC or HOPS, respectively, and are shown on all the following figures that contain waveforms, for reference. We also examine the time delay of each component of SGFs with respect to the corresponding NGF component (e.g. Ma et al.2008) if the normalized CC (not zero-lag) for that component exceeds 0.7. Positive time delay implies that the NGF arrives earlier than the SGF and the actual velocities are faster than the model velocities (at these frequencies, waveforms are primarily sensitive to VS). Synthetic tests Prior to comparing NGF and SGF tensors at The Geysers, we perform tests on ‘synthetic noise’ to evaluate the extent to which phase and relative amplitude information can be retrieved from coherency of seismic noise. Our synthetic tests are modified after synthetic tests done by Herrmann (2013b). We design a 100 km × 100 km domain with stations placed at various positions at the centre allowing us to obtain synthetic noise Green's functions (SNGFs), that is, NGFs estimated from coherency of synthetic noise, for interstation distances from 1 to 35 km (Supporting Information Fig. S4). For constructing synthetic noise records at the stations, we sum filtered (0.2–0.9 Hz) three-component velocity waveforms that are generated by randomly oriented force vectors (amplitude range −1 to +1) at random locations (but at least 50 m away from stations) on the surface with 20 sources acting together every three seconds. Different from the synthetic tests of Herrmann (2013b), we use sources overlapping in time as in Kimman & Trampert (2010) and the FKI method described earlier to compute complete single force responses on the REF model instead of surface wave responses computed by modal summation. The results of the synthetic tests are described in detail in the Supporting Information section Synthetic tests and Supporting Information Fig. S5. Based on the synthetic tests, we find that interstation coherency of velocity noise recorded at two stations compares well with normalized displacement SGFs (for an input step force) for the methodology and scenarios (distances and frequencies, layered attenuating medium) in our study (VR ≳ 80 per cent for distances ∼λ and higher VR at shorter and longer distances) under the assumption of an idealized homogenous noise distribution on the Earth's surface. We also observe better recovery of the Rayleigh wave cross-term components, RZ and ZR, compared to the RR and ZZ components, similar to results of observational (van Wijk et al.2011) and numerical studies (Haney et al.2012). Our failure to retrieve the exact GFs (i.e. VR = 100 per cent) can be attributed to the absence of depth sources, the absence of deformation rate or dipole sources, cross-terms of different phases, incomplete cancellation of overlapping sources and ubiquitous distribution of noise sources instead of their clustering near the stationary phase region (Kimman & Trampert 2010). RESULTS In the context of noise cross-correlations, ‘distances’ in the current and following sections implies interstation distances. In the comparisons of NGFs and SGFs, each SGF tensor corresponds to the velocity model that returns the highest VR with NGFs of that specific station pair, unless stated otherwise. First we discuss NGFs involving stations within the steam field or adjacent to it (i.e. all stations excluding HOPS and MNRC). Effects of sensor misorientations on NGFs In our preliminary analysis, we observed that the NGFs of some stations pairs at very small distances (<λ) had large amplitudes on the TR, TZ, RT and ZT components (hereinafter referred to as T-[R,Z] components) that are unlikely the result of 3-D velocity structure or anisotropy. These anomalous amplitudes were restricted to NGFs involving some specific stations that also returned fairly large alignment angles in the analysis of teleseismic waveforms. Fig. 6(a) shows some examples with surface stations FUM and DVB and borehole station JKB. We consider, sensor misalignments as a likely cause of these non-zero amplitudes. We use alignment angles obtained from earthquake waveforms to rotate the NGF tensors to correct orientations. For the NGFs in Fig. 6(a), rotation to the correct orientation reduces the large amplitudes on the T-[R,Z] components and significantly improves the VR between the SGFs and the NGFs (Fig. 6b). On an average, correcting the NGFs for all alignment angles > 10° increases VR between the SGFs and the NGFs, and reduces the mean root mean square amplitude of the T-[R, Z] components for 76.4 per cent and 79.8 per cent of the 407 interstation pairs, respectively. Figure 6. View largeDownload slide (a) Nine-component NGF (black) and SGF (red) tensors for specific station pairs with large anomalous amplitudes on the T-[R,Z] components. Each panel (3 × 3 group of subplots) is specific to a station pair (names of stations shown at the top of the panels). Each subplot in a panel is specific to a component pair (component pair names are at top right corner of each subplot; distance and azimuth are at the top left corner of the TZ subplots). The first station/component corresponds to the force location/direction and the second station/component corresponds to the receiver location/direction. For each station pair, the velocity model used to compute the best fitting SGF tensor is indicated at the bottom left corner of TZ subplot. For the overall tensor comparisons of a station pair, two goodness-of-fit estimates, VR and normalized zero-lag CC are mentioned at the top left corner of TR subplots. For all non-zero component pairs, the normalized CC along with the time lag between the NGF and SGF waveforms is indicated at the bottom right corners of the subplots, if the CC exceeds 0.7. For reference, the blue and green marks are the theoretical P- and S-wave arrival times, respectively, while the grey mark corresponds to 5 s after the empirical coda start time. For a particular station pair, following normalization of amplitudes of the NGF and SGF tensors, only relative amplitudes are meaningful. (b) Same as (a) but after the NGF tensors are corrected for sensor misalignments. The rotation angles for different stations are indicated on the arrows between (a) and (b). NGFs in all the following figures have been similarly corrected for sensor misorientations using alignment angles obtained from the analysis of earthquake waveforms. Figure 6. View largeDownload slide (a) Nine-component NGF (black) and SGF (red) tensors for specific station pairs with large anomalous amplitudes on the T-[R,Z] components. Each panel (3 × 3 group of subplots) is specific to a station pair (names of stations shown at the top of the panels). Each subplot in a panel is specific to a component pair (component pair names are at top right corner of each subplot; distance and azimuth are at the top left corner of the TZ subplots). The first station/component corresponds to the force location/direction and the second station/component corresponds to the receiver location/direction. For each station pair, the velocity model used to compute the best fitting SGF tensor is indicated at the bottom left corner of TZ subplot. For the overall tensor comparisons of a station pair, two goodness-of-fit estimates, VR and normalized zero-lag CC are mentioned at the top left corner of TR subplots. For all non-zero component pairs, the normalized CC along with the time lag between the NGF and SGF waveforms is indicated at the bottom right corners of the subplots, if the CC exceeds 0.7. For reference, the blue and green marks are the theoretical P- and S-wave arrival times, respectively, while the grey mark corresponds to 5 s after the empirical coda start time. For a particular station pair, following normalization of amplitudes of the NGF and SGF tensors, only relative amplitudes are meaningful. (b) Same as (a) but after the NGF tensors are corrected for sensor misalignments. The rotation angles for different stations are indicated on the arrows between (a) and (b). NGFs in all the following figures have been similarly corrected for sensor misorientations using alignment angles obtained from the analysis of earthquake waveforms. We also consider DEB and SRB that were replacement borehole stations for surface stations DVB and SSR at depths ∼150 m and ∼140 m, and at horizontal distances of ∼30 m and ∼70 m, respectively from the surface sites. Both the surface and borehole stations were simultaneously operational for ∼25 days. The NGFs and the SGFs for the two surface-borehole sensor pairs are shown in Supporting Information Fig. S6. At these small distances, GFs are primarily proportional to filtered source-time functions with the largest near-equal amplitudes on the TT, RR and ZZ components (e.g. Wielandt & Forbriger 1999). Orientation corrections to the NGFs significantly reduce the amplitudes in the off-diagonal components of these tensors (Supporting Information Fig. S6b). Although there can be other factors responsible for the off-diagonal terms, utilizing the NGF tensors may be a useful alternative approach for assessing orientations of borehole sensors. The successful validation of angles measured independently from long period teleseismic earthquake waves on higher frequency NGF tensors provides confidence in our results. Alignment angles at three stations (FUM, NEG and TCH) have also been independently obtained at the sites (Ramsey Haught, private communication, 2016), which are reasonably close to our estimates within the uncertainties (Supporting Information Table S2). In the following sections, all NGFs have been corrected using alignment angles obtained from the analysis of earthquake waveforms. NGF versus SGF at various distances Fig. 7 shows some examples of NGFs and SGFs at distances ≲2 km or ∼0.44 λ. At these small distances, we observe simple waveforms, strong NGFs, low relative amplitudes of noise trailing the ballistic phases with respect to the peak amplitudes, good fits at VR ≳ 40 per cent, and consistency in relative amplitudes of different components of the tensors. Small time delays are sometimes observed between NGFs and SGFs (e.g. AL3-AL5 shown Fig. 7); however, these delays are not consistent across all components, and based on the synthetic tests it is unclear whether they are real or due to bias introduced by the noise distribution or by filtering or processing artefacts. For NGFs at distances >2 km and ≲8 km (or >0.44λ and ≲1.7λ) shown in Fig. 8, we observe a decrease in the quality of retrieved NGFs as evidenced by the increase in relative amplitude of trailing noise. The waveforms also grow slightly more complex and the amplitude misfit between NGFs and SGFs increases, with VR decreasing to ∼20 per cent. We start to see clear inadequacies in the velocity models; for example, in HBW-HER in Fig. 8, the TT component NGF is slower than the SGF indicating that the 3-D model, which is not resolved in this region and therefore is essentially the same as the starting average 1-D model (BACK1) in this area, is too fast for this path. While these distances are too small for reliable phase measurements, large reliable time shifts measured on multiple components of the NGFs with good signal-to-noise ratio can be utilized as ‘broad-band phase delays’ in waveform tomography assuming that multiple scattering in the Earth's crust alleviates the ill-effects of inhomogeneous noise source distribution (Fichtner 2014; Lee et al.2014). Figure 7. View largeDownload slide Same as Fig. 6(b) but for station pairs at interstation distances ≤ 2 km. Figure 7. View largeDownload slide Same as Fig. 6(b) but for station pairs at interstation distances ≤ 2 km. Figure 8. View largeDownload slide Same as Fig. 6(b) but for station pairs at interstation distances > 2 km and ≤ 8 km. Figure 8. View largeDownload slide Same as Fig. 6(b) but for station pairs at interstation distances > 2 km and ≤ 8 km. Fig. 9 shows examples for distances >8 km. These are characterized by lower signal-to-noise ratio, larger amplitudes of trailing noise and large amplitudes on the TR, TZ, RT and ZT components. The simple waveforms on the TT components consist of Love waves that depend on VS structure only and show good fits, whereas the waveforms on the RR, RZ, ZR and ZZ components consisting of P–SV and Rayleigh waves depend on both, VP and VS structure and exhibit rapidly increasing misfits with distance. However, the RR components of the NGFs in Fig. 9 compare well to the SGFs. P–SV body-wave energy arriving before the theoretical S-wave arrival time can be clearly identified above the noise floor in the ZR components of AL6-NEG, AL2-STY, AL3-DES and AL2-DVB. Similar to results of Lin et al. (2008), but for a different distance and frequency range, we obtain higher amplitude coherency in the TT component than in the RR and the ZZ components at distances >2λ that compare well to the SGFs. Therefore, Love waves on the TT component NGFs should be used to constrain the shear wave velocity models in addition to Rayleigh waves whenever three-component sensors are available, as the former are generally retrieved with better signal-to-noise ratio. Similar to our synthetic tests and observations in other studies (van Wijk et al.2011), we generally find good fits between ZR and RZ component NGFs and SGFs at all distances. Figure 9. View largeDownload slide Same as Fig. 6(b) but for station pairs at interstation distances > 8 km. Figure 9. View largeDownload slide Same as Fig. 6(b) but for station pairs at interstation distances > 8 km. Compiling overall results, the average peak-to-peak amplitude ratio of the primary components (TT, RR, RZ, ZR, ZZ) of NGFs to the corresponding components of SGFs for all station pairs is ∼0.7. The standard deviation is ∼0.18 in log10 units (i.e. difference of a factor of ∼1.5, bigger or smaller) with respect to the mean amplitude ratio. To select only robust time delays between individual components of NGFs and corresponding components of the best-fitting SGFs for interpretation, we require that VR for a component pair is ≥ 35 per cent following correction for the time delay estimated from normalized cross-correlation. This criterion ensures similarity in both waveform shape and amplitude. Out of 238 station pairs that have at least 2 such components, only ∼19 per cent of station pairs have large contrasting time delays on different components (maximum time delay ≥0.1 s and minimum time delay ≤ –0.1 s) indicating possible errors in the recovered NGFs. For many station pairs such as HBW-HER in Fig. 8, large robust time delays on different components of the tensors are of the same sign, indicating an overall consistency in observations. Contribution of body waves and surface waves To study the contributions of different body-wave and surface-wave phases to the observed NGF waveforms, we compare pure surface-wave SGFs computed using the modal summation method as provided in Herrmann (2013a) to complete FKI SGFs. We use the REF velocity model for this analysis; the synthetic phase velocities and depth-dependent eigenfunctions are computed using codes provided in Herrmann (2013a). Phases that are present in the complete FKI synthetics but absent in the surface-wave synthetics are interpreted to be body-wave phases. Synthetic phase velocity dispersion curves based on the REF model for fundamental- and first higher-mode Rayleigh waves are plotted in Fig. 10(a). Fig. 10(b) shows an enlarged view of the fast arrivals at ∼3.6 km s−1 in the ZZ component NGFs at distances 17–20 km in Fig. 2. Fig. 10(c) shows SGFs at similar distances containing the fundamental-mode (‘Mode 0’) and first higher-mode Rayleigh waves (‘Mode 1’) and their combination (‘Mode 0+1’). For a source at ∼10 m depth, the peak amplitudes of the higher-mode are ∼14 times smaller than peak amplitudes of the fundamental-mode and therefore, the higher mode has a very small contribution to the total waveforms containing both modes. The amplitudes of the second higher-mode Rayleigh wave are even smaller and therefore, it is not discussed here. The recovery of higher-mode surface waves on land is more difficult than the recovery of fundamental-mode surface waves primarily due to the near-surface concentration of ambient noise sources (Kimman & Trampert 2010). However, some studies have reported success in this regard (e.g. Nishida et al.2008; Lin et al.2013). The body-wave energy observed in the FKI synthetics is considerably weaker than energy observed in the ZZ component NGFs between the theoretical P-wave and S-wave arrival times. Compared to the SGFs, the NGFs appear to be depleted in high frequencies (≳0.5 Hz) as well. More detailed analyses including the analysis of particle motions and SGFs computed with other velocity models are required to positively identify this phase in the ZZ component. However, since it has greater amplitudes and is faster (∼3.6 km s−1) than the first higher-mode Rayleigh wave expected for the REF model (phase velocity ∼3.0–3.3 km s−1), it is likely to be body-wave energy. Figure 10. View largeDownload slide (a) Synthetic phase velocity dispersion curves for the fundamental-mode and first higher-order mode Rayleigh waves for the REF velocity model. (b) Record section of the ZZ component NGFs at interstation distances of 17–20 km. For reference, the blue and green solid lines denote the theoretical P- and S-wave arrival times, respectively. The two dashed red lines mark the arrival times for apparent medium velocities of 3.6 and 2.5 km s−1. (c) REF model ZZ component SGFs at interstation distances of 17–20 km for different Rayleigh waves modes (‘Mode 0’ and ‘Mode 1’) computed using the modal summation technique, compared to complete synthetics computed using FKI. The type of modes/method is indicated at the top right corner of each subplot. ‘Mode 0’ and ‘Mode 1’ indicate fundamental and first higher-mode Rayleigh waves, respectively. ‘Mode 0+1’ synthetics include both ‘Mode 0’ and ‘Mode 1’. The range of maximum absolute amplitudes of traces in each subplot is indicated in the bottom left corner. For example, the maximum absolute amplitude of ZZ component SGF computed using FKI varies from 4.7 × 10–4 to 5.4 × 10–4 cm over distances 17–20 km. Figure 10. View largeDownload slide (a) Synthetic phase velocity dispersion curves for the fundamental-mode and first higher-order mode Rayleigh waves for the REF velocity model. (b) Record section of the ZZ component NGFs at interstation distances of 17–20 km. For reference, the blue and green solid lines denote the theoretical P- and S-wave arrival times, respectively. The two dashed red lines mark the arrival times for apparent medium velocities of 3.6 and 2.5 km s−1. (c) REF model ZZ component SGFs at interstation distances of 17–20 km for different Rayleigh waves modes (‘Mode 0’ and ‘Mode 1’) computed using the modal summation technique, compared to complete synthetics computed using FKI. The type of modes/method is indicated at the top right corner of each subplot. ‘Mode 0’ and ‘Mode 1’ indicate fundamental and first higher-mode Rayleigh waves, respectively. ‘Mode 0+1’ synthetics include both ‘Mode 0’ and ‘Mode 1’. The range of maximum absolute amplitudes of traces in each subplot is indicated in the bottom left corner. For example, the maximum absolute amplitude of ZZ component SGF computed using FKI varies from 4.7 × 10–4 to 5.4 × 10–4 cm over distances 17–20 km. Comparing surface-wave (‘Mode 0/1/0+1’) SGFs and complete FKI SGFs at distances of 14–17 km for the ZR and RR components (Supporting Information Fig. S7), we positively identify strong P-wave energy in the RR and ZR component NGFs for many station pairs (e.g. AL3-DES in Fig. 9) in light of the negligible contribution of the first higher-mode Rayleigh wave compared to the fundamental mode. NGFs with USGS 1.0 Hz sensors Fig. 11 shows GF comparisons for station pairs with one station being a USGS vertical component short period (∼1.0 Hz) station, for which only the TZ, RZ and ZZ components are available (with the vertical component station being the receiver). The NGF-SGF comparisons for these station pairs are important as four of these stations (GAXB, GBG, GGPB and GSG) lie well outside the reservoir area (>3 km away from the reservoir outline and the LBNL network) and none of our 1-D velocity models are calibrated to fit these paths. Moreover, a single background 1-D model is assumed in model G3D1 for almost all regions outside the reservoir area and an artificial, unrealistic boundary exists between the 1-D and the 3-D models along the reservoir outline, except to the northwest. Model REF, which is a 1-D average of a 3-D model of NW Geysers, fits most paths to stations GPM and GSG that lie to the NW and NE of The Geysers, respectively. However, large time delays, −0.4 to −0.75 s, are clearly observed at station GSG (see waveforms for ACR-GSG, JKR-GSG, STY-GSG, FUM-GSG in Fig. 11) consistent with body-wave tomography studies that have determined lower velocities to the northeast of The Geysers attributed to rocks of the Clear Lake volcanics and interbedded sedimentary deposits of the Great Valley sequence (Hearn et al.1981; Eberhart-Phillips 1986; Julian et al.1996). The delays are also consistent with a 0.25 s P-wave station correction at GSG with respect to GBG (Eberhart-Phillips & Oppenheimer 1984). The retrieval of robust NGFs from these sensors demonstrates their potential in supplementing temporary deployments of higher quality broad-band sensors for ambient noise Rayleigh wave tomography studies (e.g. Porritt et al.2011). This is especially true for regions like Northern California, where they are a significant part of pre-existing seismic infrastructure in the Northern California Seismic Network (http://www.ncedc.org/ncsn/, last accessed July 2016). Figure 11. View largeDownload slide Same as Fig. 6(b), but for all distances and stations pairs with one USGS vertical component 1.0 Hz sensor. Interstation distance and azimuth, and the overall goodness-of-fit estimates are at the top of TZ subplots. Figure 11. View largeDownload slide Same as Fig. 6(b), but for all distances and stations pairs with one USGS vertical component 1.0 Hz sensor. Interstation distance and azimuth, and the overall goodness-of-fit estimates are at the top of TZ subplots. We were also able to retrieve robust NGFs with the accelerometer at station DRH (Supporting Information Fig. S8) that compare very well to SGFs at distances up to 20 km. The RR component NGFs for station pairs CLV-DRH, HVC-DRH and ACR-DRH also show strong P-waves (see Supporting Information Fig. S7a for comparison). Low gain accelerometers have been successfully used in some ambient noise tomography studies (e.g. Cho et al.2007). CONTRIBUTION OF 3-D STRUCTURE AND ANISOTROPY Even after the orientation corrections, the NGFs for some longer distance paths (≳8 km or ∼1.7 λ) show significant non-zero amplitudes on the T-[R,Z] components that appear robust with respect to the trailing noise (Fig. 12). These waveforms are reproduced well by the SGFs computed using the 3-D velocity model indicating the effects of significant 3-D structure at distances ≳2 λ. We visually identified these interstation paths, which appear to be approximately parallel to the southwestern boundary of the reservoir area along the extent of the high-velocity region in the southeast section of the reservoir. While most of these paths involve stations SSR and SRB, the teleseismic wave analysis described in the Supporting Information indicates that the sensor installed at SSR was correctly oriented. It is possible that off-great-circle wave propagation or multipathing effects along the large velocity contrast cause these large amplitudes on the T-[R,Z] components. However, the 3-D velocity model is poorly resolved outside the reservoir boundary (Gritto et al.2013a). There are other interstation paths for which we see robust non-zero amplitudes on these components but they are not reproduced by the 3-D velocity model. Figure 12. View largeDownload slide Same as Fig. 6(b), but for some station pairs with significant non-zero amplitudes on the T-[R,Z] components of NGFs even after corrections for sensor misalignments. These amplitudes are reproduced well by the 3-D velocity model. Figure 12. View largeDownload slide Same as Fig. 6(b), but for some station pairs with significant non-zero amplitudes on the T-[R,Z] components of NGFs even after corrections for sensor misalignments. These amplitudes are reproduced well by the 3-D velocity model. Assuming far-field noise source incidence only, if the station pairs are not aligned with the dominant noise source incidence direction, non-zero amplitudes can be expected in T-[R,Z] components of NGFs (Roux 2009; Durand et al. 2011). It is possible to minimize the amplitudes in these components by rotating the reference frames at the two stations so that the new reference frames are approximately aligned with the azimuth of dominant noise source incidence direction [the method is referred to as Optimal Rotation Algorithm (ORA); Roux 2009]. We quantify the relative amplitudes of T-[R,Z] components of the NGFs at The Geysers, in terms of the misfit parameter φ, defined as the ratio of sum of squares of amplitudes of the four off-diagonal components to the sum of squares of all components of the NGF tensor (eq. 17 in Saade et al. 2015). We examined φ as a function of azimuth of the station pairs to investigate any systematic relationship. Our results, which are fully described in the Supporting Information section Application of Optimal Rotation Algorithm, indicate that the non-zero signals in T-[R,Z] components of our NGFs at large distances are just leading or trailing noise in the NGF time-series rather than coherent signals resulting from systematic azimuthal deviation of interstation paths from a particular noise source direction. It is possible that seismic anisotropy also contributes to the non-zero amplitudes in these components at larger distances by giving rise to quasi-Love or quasi-Rayleigh waves (Maupin & Park 2007). The shallow subsurface (top ∼3–5 km) at The Geysers can be considered as a horizontal transversely isotropic (HTI) medium with an average ∼4 per cent VS anisotropy (up to ∼11 per cent at some locations) estimated from shear-wave splitting measurements on local earthquakes records (Majer et al.1988; Elkibbi & Rial 2005; Elkibbi et al.2005). The VS anisotropy is defined as $$\frac{{{V_{S,{\rm{fast}}}} - {V_{S,{\rm{slow}}}}}}{{{V_{S,{\rm{fast}}}}}}$$ × 100 as in Elkibbi et al. (2005). The anisotropy at The Geysers is believed to be caused by stress-induced alignment of fractures and cracks. While the fast axes are generally parallel or sub-parallel to the N-to-NE direction of the regional maximum compressive stress, a smaller set of NW–SE oriented fast axes are also observed, primarily in the southeast Geysers and likely related to local fault shearing effects (Elkibbi & Rial 2005; Elkibbi et al.2005). The 3-D model of Gritto et al. (2013a,b) was inverted utilizing a well distributed set of earthquakes and stations, and therefore can be considered as an isotropic average of the velocity structure at The Geysers, with a maximum of 19 per cent ± 7 per cent variation in VS at depths between 0.5 km and 4.0 km. Therefore, seismic wave propagation at The Geysers at large distances is possibly controlled, to first order, by the average 3-D isotropic structure in the reservoir. In order to test the amplitudes caused by anisotropy and expected in the T-[R,Z] components of the NGFs, we also compute SGFs for station pairs at The Geysers incorporating an HTI medium in the 1-D velocity models. The 3-D finite-difference seismic wave propagation code SW4 provides the capability to model seismic waveforms for layered anisotropic but purely elastic media. Therefore, we restrict this analysis to stations pairs at distances <10 km, so that we can ignore the effects of anelastic attenuation (further discussed in section Amplitude decay). We also ignore the considerable scatter in the degree of anisotropy and in the fast axis directions across The Geysers that were observed by Elkibbi & Rial (2005) and Elkibbi et al. (2005). We use the formulation of Hudson (1981, 1982) as provided in eqs (1) to (4) in Crampin (1984) to estimate the five anisotropic elastic constants for an HTI medium—A, C, F, L and N in Love notation (Saade et al.2015). The differences between (A, C) and (L, N) determine the VP and VS anisotropy, respectively. A review on the $$\eta \ = \frac{F}{{A - 2L}}\ $$ parameter can be found in Kawakatsu et al. (2015). The effective anisotropic elastic tensor in a cracked medium can be modelled as the sum of the elastic tensor for the uncracked solid and first and second order perturbation terms that depend on crack density (ε), crack aspect ratio (d), Lame's constants for the uncracked solid, and Lame's constants for the weak crack inclusions (λ', μ'). This formulation assumes a weak distribution of parallel disconnected penny-shaped cracks with dimensions and spacing much smaller than the seismic wavelengths under consideration and is valid for small values of ε < 0.1 (Crampin 1984; Peacock & Hudson 1990; Cheng 1993). The Geysers is characterized by a vapour-dominated geothermal field with a considerable volume of water injected into the reservoir. For this study, we assume that the cracks are filled with water ( VP = 1.5 km s− 1, λ' = 2.25e9 Nm,  μ' = 0). It is important to note that the presence of steam in water can drastically reduce VP of the two-phase mixture, especially if the water and steam phases are in thermodynamic equilibrium. However, this effect is expected to decrease with increasing pressures and temperatures (Liu & Kieffer 2009). Majer et al. (1988) investigated one site at The Geysers and found little or no evidence for VP anisotropy. Therefore, we assume thin cracks, d ∼ 0, which gives A = C. An assumed value of ε ∼ 0.036 leads to an approximate VS anisotropy ∼4 per cent in the layers of the top ∼3.1 km of the 1-D velocity models. We obtain values of η ∼ 0.82–0.89, compared to η = 1.1 as adopted by Saade et al. (2015) based on scaling relationships (Montagner & Anderson 1989). We test two directions of crack normals, N60° W and N125° W, transverse to the average directions of the fast axes determined by Elkibbi & Rial (2005). For the computation in SW4, the crack normal is assumed to be in direction 1 and the station paths are rotated accordingly. Fig. 13(a) shows NGF and SGF comparisons for 2 station pairs at distance < 10 km for which the significant non-zero amplitudes on the T-[R,Z] components of the NGFs are reproduced well by the SGFs computed using the 3-D velocity model. Fig. 13(b) shows the same NGFs as Fig. 13(a) but with SGFs computed with an elastic (no anelastic attenuation) and anisotropic version of the 1-D velocity model VSP0, incorporating ∼4 per cent VS anisotropy in the top ∼3.1 km layers with crack normal direction N125° W for station pairs FUM-SSR and DRK-SSR. The VR and the amplitudes on the T-[R,Z] components for the 1-D anisotropic models are clearly smaller than those for the 3-D model. The 1-D anisotropic model SGFs yield higher VR compared to the isotropic 1-D or 3-D models for only 13 out of ∼270 nine-component station pairs (∼5 per cent) at distances < 10 km; for 6 out of 13 station pairs, the improvement in VR was negligible (<1 per cent). In the case of the SGFs dominated by surface waves in an anisotropic medium, ORA can also be used to minimize the amplitudes in the T-[R,Z] components (Durand et al.2011; Saade et al.2015). We also perform a synthetic test, in which we compute SGFs assuming the anisotropic REF model with one reference station at the centre of a circle and other stations placed at radial-distance intervals of 0.5 km up to a total distance of 10 km and 2.0° azimuthal spacing. For these GFs, the maximum value of φ ranges from ∼4e-4 to ∼1e-2 for distances between 0.5 to 10 km. These values are smaller than φ for our NGFs at The Geysers (Supporting Information Fig. S9a) at distances <10 km. This observation, combined with the failure of ORA to significantly reduce φ in our NGFs (Supporting Information Fig. S9b), indicates that the large residual amplitudes in the T-[R,Z] components of the NGFs for most station pairs (<10 km) at The Geysers correspond to just leading or trailing noise in the NGF time-series. These non-zero signals are unlikely to be real coherent signals resulting from systematic effects of non-uniform illumination from far-field noise sources or currently known subsurface anisotropy in the geothermal field. For at least some station pairs that have robust signals in the T-[R,Z] components, the 3-D model provides better fits than the anisotropic 1-D models. We leave analysis of the effects of anisotropy at distances > 10 km for future studies. Figure 13. View largeDownload slide (a) Same as Fig. 12 but for two different station pairs at interstation distances <10 km. (b) Same as (a) but for SGFs computed using an elastic and anisotropic version of the 1-D velocity model VSP0 with ∼4 per cent VS anisotropy to a depth of ∼3.1 km and crack normal in the direction N125° W. ‘–A’ has been added to the model names to indicate anisotropy. SGFs computed using the 3-D velocity model fit the significant non-zero amplitudes on the T-[R,Z] components of NGFs better than SGFs computed using the 1-D anisotropic velocity model. Figure 13. View largeDownload slide (a) Same as Fig. 12 but for two different station pairs at interstation distances <10 km. (b) Same as (a) but for SGFs computed using an elastic and anisotropic version of the 1-D velocity model VSP0 with ∼4 per cent VS anisotropy to a depth of ∼3.1 km and crack normal in the direction N125° W. ‘–A’ has been added to the model names to indicate anisotropy. SGFs computed using the 3-D velocity model fit the significant non-zero amplitudes on the T-[R,Z] components of NGFs better than SGFs computed using the 1-D anisotropic velocity model. EVALUATION OF VELOCITY MODELS Fig. 14 shows a summary of ∼385 interstation noise coherency paths with the best-fitting VR between the NGFs and the SGFs ≥ –20 per cent, grouped by the best-fitting velocity models, that is, the models that provide the SGFs with the highest VR for the NGFs along these paths. The VR threshold of −20 per cent is a quality check on both, the quality of the NGFs (poor quality NGFs will return poor VR upon comparison with SGFs), and the ability of velocity models to fit the NGFs (paths along which the velocity models are found to be deficient [e.g. paths to station GSG] are filtered out). Out of the four 1-D models tested, REF provides the highest VR for the most number of paths, primarily across NW and central Geysers (Fig. 14a). This can be attributed to its origin from a 3-D model of the NW Geysers region. Similarly, VSP0, which corresponds to a velocity profile in the SE Geysers, is the best-fitting velocity model for most NGFs across the SW- and the SE Geysers (Fig. 14c). The 3-D model G3D1 provides the best-fitting SGFs for a slightly higher number of NGF paths compared to REF (Fig. 14e). To find regions where the 3-D model fits the NGF waveforms significantly better than the 1-D velocity models in a relative sense, we examine paths for which the VR of the 3-D model SGFs exceeds those of the 1-D models by an arbitrarily chosen threshold of 7 per cent or more (Fig. 14f). Most of these paths are across the transition zone between the slower NW section and the faster SE section of The Geysers. The 3-D model also incorporates a prominent low VP/VS ratio anomaly in this area with estimates (down to ∼1.42) considerably lower than those in the 1-D models (Figs 3 and 4). Therefore, it is not surprising that the 3-D model fits long distance paths across this complex region better than the 1-D models. Models AVG1 and BACK1 provide minor improvements over other velocity models for a smaller number of paths in the NW Geysers (Figs 14b and d). As expected, the 1-D models fit most of the paths to the stations outside the reservoir area (stations GAXB, GBG, GGPB, GSG) better than G3D1 (59 out of 62). Figure 14. View largeDownload slide First five subplots (a–e) show interstation paths grouped by best-fitting velocity models. The velocity model and the number of paths are mentioned in the top and bottom left corner, respectively. Subplot (f) ‘3D BEST’ shows paths for which VR of NGFs with the 3-D model SGFs exceed the VR with the 1-D model SGFs by 7 per cent or more (34 out of 124). For all the paths shown in this figure, the best-fitting VR is ≥–20 per cent. Figure 14. View largeDownload slide First five subplots (a–e) show interstation paths grouped by best-fitting velocity models. The velocity model and the number of paths are mentioned in the top and bottom left corner, respectively. Subplot (f) ‘3D BEST’ shows paths for which VR of NGFs with the 3-D model SGFs exceed the VR with the 1-D model SGFs by 7 per cent or more (34 out of 124). For all the paths shown in this figure, the best-fitting VR is ≥–20 per cent. Next we examine phase differences represented by robust time delays (defined in section NGF versus SGF at various distances) between NGFs and SGFs computed with different velocity models for some station pairs. We first consider NGFs and SGFs for intermediate distance (between 2 km and 12 km) interstation paths in ∼NW and Central Geysers (Figs 14a and b), for which REF and BACK1 are the best-fitting velocity models. Barring some large outliers (<–1.0 s or >1 s; ≲2 per cent of measurements removed), the robust time delays for these pairs of NGFs and SGFs are ∼–0.02 ± 0.17 s (mean ± 1 standard deviation; for ∼240 measurements). For NGFs along the same paths, SGFs computed using the faster velocity model VSP0 lead to time delays of ∼–0.12 ± 0.18 s, indicating slightly earlier arrival of SGFs with respect to NGFs. Similarly, time delays between NGFs between SGFs computed with the best-fitting model VSP0 for intermediate distance paths in ∼SE Geysers (Fig. 14c) are ∼+0.12 ± 0.11 s (for ∼75 measurements). The actual velocities could be slightly greater than VSP0 velocities. For NGFs along the same paths, SGFs computed using the slower velocity model REF lead to time delays of ∼+0.26 ± 0.16 s, indicating more delayed arrival of SGFs with respect to NGFs. Our interstation paths in Figs 14(a)–(c) are well distributed in azimuth which is expected to help reduce the effects of phase errors in NGFs due to non-uniform noise source distribution. The meaningful grouping of interstation paths in various sub-regions of The Geysers (Figs 14a–c) indicates an overall consistency of the results. NGFs WITH BDSN STATIONS HOPS AND MNRC We first analyse the simple waveforms of Love waves on the TT components of the NGFs from BDSN stations HOPS and MNRC to stations at the periphery of the reservoir area to gain insight into the shear wave velocity structure in the region to the northwest and to the east of The Geysers (Fig. 15a). The REF model SGFs are too fast compared to the NGFs (Figs 15b and c), which is consistent with body-wave traveltime studies that suggest that the reservoir area has faster than regional seismic velocities (Majer & McEvilly 1979; Eberhart-Phillips 1986). The SGFs of the regional velocity model, GIL7 (Fig. 3) can reasonably fit the Love waves on the paths from the eastern boundary of the reservoir to MNRC (Fig. 15b). GIL7 SGFs are also consistent with the NGFs of paths from some of these stations to station GSG (Fig. 15e), for which the REF model was determined to be too fast (Fig. 11). But we find that GIL7 is too slow for paths from the northwestern boundary of the reservoir area to HOPS (not shown here) with broad-band phase delays from ∼1.3 s to ∼2.1 s. Therefore, we iteratively perturb the VS values in the top ∼5 km of GIL7 to improve the waveform fits of the TT components of the NGFs for these paths. Primarily deceasing the VS values in the top ∼4 km layers of GIL7 by ∼0.1–0.3 km s−1 (model GIL7.1 shown in Fig. 3) produces satisfactory waveform fits (Fig. 15c). In addition to Love waves, clear SH waves can also be identified in the NGFs of some paths to station HOPS (from HBW, RGP, MCL, BRP) that compare well to the SGFs. Figure 15. View largeDownload slide (a) Map showing interstation paths (red and green lines) from stations on the periphery of the reservoir area (white polygon) to BDSN stations HOPS and MNRC (red triangles), respectively. Meaning of other symbols is same as in Fig. 1; (b) TT component NGFs (black traces) and SGFs (red traces) computed using REF (left subplot) and GIL7 (right subplot) velocity models for paths from stations at the eastern boundary of The Geysers to station MNRC arranged in the order of increasing interstation distance. Station names and azimuths are shown near the traces; (c) same as (b) but for paths from stations at the northwestern boundary of The Geysers to HOPS. The velocity models are REF (left subplot) and GIL7.1 (right subplot); (d) Observed Love wave phase velocity dispersion curves measured on TT component NGFs for paths in Figs 15(a)–15(c) plotted against phase velocities predicted for models REF, GIL7.1 and GIL7; (e) same as Fig. 11 but for station pairs including GSG and stations at the northeastern boundary of the reservoir and using velocity model GIL7; (f) same as Fig. 9 but for a station pair including MNRC. Figure 15. View largeDownload slide (a) Map showing interstation paths (red and green lines) from stations on the periphery of the reservoir area (white polygon) to BDSN stations HOPS and MNRC (red triangles), respectively. Meaning of other symbols is same as in Fig. 1; (b) TT component NGFs (black traces) and SGFs (red traces) computed using REF (left subplot) and GIL7 (right subplot) velocity models for paths from stations at the eastern boundary of The Geysers to station MNRC arranged in the order of increasing interstation distance. Station names and azimuths are shown near the traces; (c) same as (b) but for paths from stations at the northwestern boundary of The Geysers to HOPS. The velocity models are REF (left subplot) and GIL7.1 (right subplot); (d) Observed Love wave phase velocity dispersion curves measured on TT component NGFs for paths in Figs 15(a)–15(c) plotted against phase velocities predicted for models REF, GIL7.1 and GIL7; (e) same as Fig. 11 but for station pairs including GSG and stations at the northeastern boundary of the reservoir and using velocity model GIL7; (f) same as Fig. 9 but for a station pair including MNRC. We also compare the Love wave phase velocities extracted from the TT component in the NGFs for these paths (in Figs 15b and c) with synthetic phase velocities for the 1-D velocity models (Fig. 15d). We first apply multiple filter analysis (Dziewonski et al.1969; Herrmann 1973, 2013b) on the NGFs to measure group velocities at periods ∼1.1 to ∼5.0 s (or ∼0.2 to ∼0.9 Hz) and then determine phase velocities using the reference dispersion curves to resolve the 2πN ambiguity, in which N is an integer (Bensen et al.2007; Lin et al.2008). In this frequency passband, observed phase velocities for paths to station HOPS are clearly greater than those for station MNRC by ∼0.1 to ∼0.3 km s−1, and they agree well with phase velocities predicted by GIL7.1 and GIL7 models, respectively. All of these phase velocities are slower than the REF model velocities by ∼0.3 to ∼0.7 km s−1. While it is possible that other velocity models (different from GIL7 and GIL7.1) might provide equal or better fits to these NGFs, it appears that to the first order, shallow crustal shear wave velocities in the region to the northwest of The Geysers are higher than those to the east of The Geysers, and they are both lower than velocities within the reservoir area. Comparing the RR, ZZ, ZR and RZ components suggests that the NGFs at most stations at The Geysers with BDSN stations HOPS and MNRC at distances >30 km are characterized by complicated waveforms, poor fits with synthetic waveforms of GIL7 and its modified versions, and higher relative amplitudes of trailing noise. However, for some stations pairs (e.g. FUM-MNRC), good quality NGFs are retrieved that show remarkable similarity to the SGFs computed with GIL7 with slightly modified values (VS in top 3 layers increased by 100, 50 and 20 m s−1, respectively) in the shallow crust (Fig. 15f). Note the large amplitudes in the TR and RT component of the NGFs in Fig. 15(f), in which FUM has already been corrected for sensor misalignment (Fig. 6). AMPLITUDE DECAY We study the decay of ground motion amplitudes in the reservoir area by analysing the NGF amplitudes as a function of frequency f and interstation distance r between 1 and 25 km. Observations and interpretation of TT component amplitudes First, we focus on the TT component as they consist of relatively simple SH and Love waves, and usually have the largest relative amplitudes among all components of NGFs (see Figs 6–9, 12 and 13, and Supporting Information Figs S8 and S9c). We use the NGF waveform segments in the data window prior to the causal filtering and the tensor amplitude normalization steps described in subsection Procedure for comparison between NGFs and SGFsin the Methodology section. We extract spectral amplitudes from the smoothed Fourier amplitude spectra of the waveform segments at multiple frequencies, f = 0.25 Hz, 0.42 Hz and 0.72 Hz. REF model Love wave group velocities UL(f) at these frequencies are ∼2.44, 2.29 and 2.22 km s−1, respectively. The Fourier Spectral Amplitudes (FSAs) are extracted by interpolating the spectra at 7 points around the frequency of interest with half-width ∼0.1 Hz and averaging the amplitudes with a 7-point Hanning window. We compare the NGF FSAs with the FSAs of the equivalent SGFs of the 3-D and the 1-D REF models that are computed the exact same way. A second group of FSAs of raw synthetic Green's functions (RSGFs) is extracted from the raw unfiltered velocity GFs of the same duration prior to the integration, zero-phase bandpass filtering, time-reversal and symmetric component steps in subsection Procedure for comparison between NGFs and SGFs (Supporting Information Fig. S3). This allows us to evaluate the degree of contamination of the SGF amplitudes with respect to the ‘true’ amplitudes and compare the decay characteristics better to earthquake ground-motion amplitudes. Without the contamination, the distance scaling of RSGF spectral amplitudes should be approximately similar to that of NGF amplitudes in a narrow frequency passband. To distinguish the effects of distance-dependent geometrical spreading, G(r), and total anelastic attenuation, represented by the quality factor Q(f), we also compare the FSAs of SGFs computed for velocity models with weak and strong anelastic attenuation. SW4 allows computation of synthetic seismograms for purely elastic 3-D velocity models. Synthetic seismograms for the weak anelastic attenuation 1-D models were computed by setting QP and QS to high values (1000 and 500, respectively) in all layers. For the strong anelastic attenuation 1-D models, QP and QSare set to ∼15 for depths < 5.8 km and ∼25 for greater depths. In this paper, we primarily focus on FSAs at 0.72 Hz, as the effects of anelastic attenuation are expected to be more prominent at higher frequencies. Fig. 16(a) shows the TT component of the RSGF FSAs at 0.25 Hz as a function of r. The scatter in the amplitudes for the 3-D velocity model is caused by smooth 3-D variations in the seismic velocities. At 0.25 Hz and for r < 25 km, the ground motion decay is very similar for both, low Q and high Q models, implying that the effect of anelastic attenuation can be neglected at these low frequencies and small distances. Directly approximating G(r) by the amplitude decay, we obtain G(r)  ∼  r− 0.82 to G(r)  ∼  r− 1.14, implying dominance of body waves. This decay is stronger than the surface wave geometrical spreading(r− 0.5) expected at distances ≳1λ–2λ (1λ ∼ 11.4 km for Love wave phase velocity cL∼2.85 km s−1 at ∼0.25 Hz) for surface sources. We also plot the RSGF FSAs of the 3-D elastic model and 1-D weak attenuation AVG1 model at 0.72 Hz. AVG1 is an average 1-D model of the 3-D velocity model within the reservoir area. Using various approximations, Menon et al. (2014) have shown that azimuthally averaging the coherence estimates of different station pairs for a 3-D inhomogeneous but elastic velocity structure can introduce apparent anelastic attenuation in the decay of average coherence amplitudes if one were simply fitting an exponential decay model to the amplitudes. The similarity in decay characteristics of the FSAs obtained from the average 1-D and 3-D weak attenuation models at both 0.25 and 0.72 Hz suggests that it might be possible to detect amplitude decay caused by anelastic attenuation following proper averaging of a 3-D velocity structure across all paths (including the ones outside the reservoir), at least for the wavelengths and heterogeneity scales analysed in this study. Figure 16. View largeDownload slide (a) TT component RSGF FSAs for different velocity models and frequencies on log10 scale as a function of distance. The decay curves are scaled so that they have similar amplitudes at ∼1–2 km. 1-D and 3-D indicate FSAs of synthetics computed with 1-D and 3-D velocity models, respectively. ‘No Q’ and ‘Q∼15’ indicate synthetics computed with weak and strong attenuation models, respectively. 1-D RSGF decay curve at 0.72 Hz used the AVG1 model unlike other 1-D model decay curves that use the REF model. (b) TT component NGF and SGF FSAs at 0.25 Hz as a function of interstation distance. The NGF FSAs (grey +) are uniformly scaled up by a constant such that the maximum value is 0. The frequency is indicated at the top left corner of the plots. Black diamonds and error bars represent the mean of NGF amplitudes and their standard deviation ± 1 σ in bins for all bins with more than five data points. The bin widths are 1 and 2 km for distances ≤10 km and >10 km, respectively. ‘SGF 3D’ and ‘SGF 1-D’ are FSAs of equivalent SGFs of the 3-D and 1-D REF velocity models, respectively, plotted for the same station pairs. The FSAs of SGFs are uniformly scaled by a constant to minimize the L1 norm of their difference with the binned NGF FSAs. (c) Similar to (b) but for FSAs at 0.72 Hz. We also plot synthetic FSAs of weak and strong attenuation models. The SGF decay curves are scaled such that their amplitude at ∼1.5 km is similar to the average amplitude of NGF FSAs in the 1–2 km range. Scaling the binned NGF decay curve instead, to directly fit the REF model SGF decay curve by minimizing the L1 norm (black dashed line), doesn't change the absolute amplitudes significantly. The RSGF decay curve (red dashed line) is scaled to fit the REF model SGF decay curve at distances ≳ 6 km. (d) The data (green +) are tangential component FSAs of earthquake records at 0.72 Hz and the red diamonds and error bars represent the mean amplitudes and their standard deviation ± 1 σ in the bins. For comparison, we also plot binned TT component NGF FSAs at 0.72 Hz (black diamonds and error bars; same as [c]) that are scaled to minimize L1 norm of binned FSAs at distances ≳ 3 km (∼1 wavelength). ‘N’ at the top right corner indicates the number of data points. Figure 16. View largeDownload slide (a) TT component RSGF FSAs for different velocity models and frequencies on log10 scale as a function of distance. The decay curves are scaled so that they have similar amplitudes at ∼1–2 km. 1-D and 3-D indicate FSAs of synthetics computed with 1-D and 3-D velocity models, respectively. ‘No Q’ and ‘Q∼15’ indicate synthetics computed with weak and strong attenuation models, respectively. 1-D RSGF decay curve at 0.72 Hz used the AVG1 model unlike other 1-D model decay curves that use the REF model. (b) TT component NGF and SGF FSAs at 0.25 Hz as a function of interstation distance. The NGF FSAs (grey +) are uniformly scaled up by a constant such that the maximum value is 0. The frequency is indicated at the top left corner of the plots. Black diamonds and error bars represent the mean of NGF amplitudes and their standard deviation ± 1 σ in bins for all bins with more than five data points. The bin widths are 1 and 2 km for distances ≤10 km and >10 km, respectively. ‘SGF 3D’ and ‘SGF 1-D’ are FSAs of equivalent SGFs of the 3-D and 1-D REF velocity models, respectively, plotted for the same station pairs. The FSAs of SGFs are uniformly scaled by a constant to minimize the L1 norm of their difference with the binned NGF FSAs. (c) Similar to (b) but for FSAs at 0.72 Hz. We also plot synthetic FSAs of weak and strong attenuation models. The SGF decay curves are scaled such that their amplitude at ∼1.5 km is similar to the average amplitude of NGF FSAs in the 1–2 km range. Scaling the binned NGF decay curve instead, to directly fit the REF model SGF decay curve by minimizing the L1 norm (black dashed line), doesn't change the absolute amplitudes significantly. The RSGF decay curve (red dashed line) is scaled to fit the REF model SGF decay curve at distances ≳ 6 km. (d) The data (green +) are tangential component FSAs of earthquake records at 0.72 Hz and the red diamonds and error bars represent the mean amplitudes and their standard deviation ± 1 σ in the bins. For comparison, we also plot binned TT component NGF FSAs at 0.72 Hz (black diamonds and error bars; same as [c]) that are scaled to minimize L1 norm of binned FSAs at distances ≳ 3 km (∼1 wavelength). ‘N’ at the top right corner indicates the number of data points. Fig 16(b) compares NGF and SGF FSAs at 0.25 Hz. Since the absolute amplitudes of the NGFs are unknown, we scale the FSA decay curves by increasing or decreasing all the data in the log scale by some constant for comparison. The shape of the decay of the SGF amplitudes is different from the decay seen for the RSGF amplitudes, because of the amplitude contamination as explained previously. The degree of contamination will be lower for higher frequencies and for gently decaying band-limiting tapers applied to cross-correlation spectral amplitudes such as a cosine taper. While there is considerable scatter in the NGF amplitudes, the mean amplitudes binned by distance show systematic decay with increasing distance. The variability in the amplitudes, indicated by their standard deviation in the bins (∼0.3 in log10 units or a factor of ∼2), is likely caused by significant differences in time periods over which daily NGFs of different station pairs were stacked (∼25 d to ∼3 yr) or by effects of inhomogeneous noise source distribution or even by variability in the coupling of the sensors with the ground at different stations. For example, in Fig. 16(b), amplitudes of NGFs for pairs with station JKR are systematically lower than mean amplitudes at the same distances. The amplitude decay characteristics of the NGFs provide useful information specifically for surface sources and they are expected to be free from earthquake radiation and source effects (but not from effects of inhomogeneous noise source distributions) that are likely present in earthquake ground motions at low frequencies and small epicentral distances. While the NGF amplitudes generally exhibit similar decay behaviour as the SGF amplitudes, they appear to be systematically higher at r ≳ 10 km. Fig. 16(c) compares SGF FSAs at 0.72 Hz for a variety of 1-D and 3-D models with NGF FSAs. In this figure, the decay curves are scaled such that the mean of the NGF FSAs at 1–2 km distance is similar to the SGF FSAs of different models at ∼1.5 km. The decay of mean NGF amplitudes binned by distance is inside the domain spanned by the SGF amplitudes for 1-D velocity models with weak and strong anelastic attenuation. The presence of NGF amplitudes well outside this domain would have led to the obvious conclusion that amplitude decay characteristics obtained from coherency of the ambient noise at The Geysers are physically unrealistic and incorrect. If we assume that the relative interstation NGF amplitudes are correct and the REF model is a realistic 1-D representation of the velocity structure at The Geysers, the difference between the average NGF FSAs binned by distance and the SGF FSAs for the weak attenuation REF model at f = 0.72 Hz can be attributed to anelastic attenuation. Fitting this difference in amplitudes by a factor of e− αr, where the attenuation coefficient α is defined as $$\alpha \ = \frac{{\pi f}}{{Q( f )U( f )}}\ $$, we obtain α ∼ 0.03 km–1 and Q ∼ 33 at 0.72 Hz (in Herrmann 2013a, the attenuation coefficient's symbol is γ instead of α). The corresponding values of α and Q at 0.42 Hz are ∼0.02 km–1 and ∼27, respectively. By employing synthetic amplitudes from a realistic low-attenuation model as reference, we avoid any assumption regarding the functional form of geometrical spreading which can strongly depend on crustal structure (Burger et al.1987; Bowman & Kennett 1991). Notwithstanding the scatter in the NGF FSAs, we find satisfactory agreement between the average NGF and the REF model SGF FSA decay curves at 0.72 Hz in Fig. 16(c). Scaling the average NGF decay curve to directly fit the REF-model SGF decay curve by minimizing the L1 norm between the amplitudes does not lead to any significant change in amplitudes. We reach similar conclusions for FSAs at 0.42 Hz (Supporting Information Fig. S10a). Measured values of α and Q are also comparable to Love wave α values predicted by the REF-model (∼0.022 and ∼0.04 km–1 at 0.42 and 0.72 Hz, respectively) and constant S-wave Q values in the REF-model (∼25 to ∼40 in the top ∼4.5 km), respectively. The similarity between the RSGF and the SGF FSAs at distances r≳ 3 km shows that the contamination of the NGFs from the acausal amplitudes of the anti-causal component can be ignored at higher frequencies and large distances. For FSAs at 0.42 Hz, they are similar at distances r ≳ 8 km (Supporting Information Fig. S10a). Comparison with earthquake ground motion amplitudes We also examine the path attenuation of horizontal component ground motions of earthquakes at The Geysers. The details of the processing of earthquake records are described in  Appendix A3. Supporting Information Fig. S11 shows examples of waveforms and Fourier amplitude spectra of an earthquake. Fig. 16(d) shows tangential component earthquake FSAs at 0.72 Hz normalized by seismic moments of individual earthquakes calculated using MD–M0 relationships valid for central California (Bakun 1984). At these low frequencies, MW < 3 earthquakes can be treated as point sources and effects of source spectrum shape can be ignored. The earthquake ground motion amplitudes are not theoretically comparable to NGF amplitudes as they are a function of force-couple GFs (unlike single force GFs recovered by ambient noise cross-correlation), earthquake radiation patterns, earthquake depths and station azimuths. At distances ≳ 1λ − 2λ, the seismic phases that dominate waveforms of empirical surface-focus GFs and shallow earthquakes should be similar. Therefore, if the empirical GFs have been correctly retrieved from ambient noise cross-correlations in terms of the relative amplitudes, they should exhibit similar amplitude decay behaviour as shallow-earthquake ground motions at far-field distances (e.g. Zhang & Yang 2013). Scattering at high frequencies and averaging across multiple paths and azimuths alleviates the effects of earthquake source radiation pattern. Most of the earthquakes used in our study are very shallow (94 out of 121 earthquakes are < 3 km deep; all earthquakes < 4 km deep; Supporting Information Table S3) with respect to the shortest wavelength examined in this section (λ ∼ 3.2 km at ∼0.72 Hz). In Fig. 16(d), significant variability (standard deviation ∼0.3 log10 units) can be observed in earthquake FSAs similar to other studies (e.g. Atkinson 2004). It is likely that many factors including effects of 3-D structure, site amplification, and location and magnitude uncertainties (∼0.2 units; Peggy Hellweg, private communication, 2017) contribute to the scatter observed in the amplitudes. However, we are find that the mean earthquake FSAs binned by distance show similar decay as the NGF FSAs at 0.72 for distances ≳ 3 km. This provides strong evidence that our NGFs are recovering realistic path attenuation of ground motions. At closer distances, earthquake and NGF amplitudes are not comparable, because NGF amplitudes are contaminated and earthquake ground motion amplitudes are strongly dependent on depth and the source radiation pattern. Mean FSAs at slightly lower frequency 0.42 Hz seem to agree with each other at distances ≳ 4 km (Supporting Information Fig. S10b). Following the large difference between the shape of the TT component SGF and RSGF attenuation curves at 0.25 Hz (Figs 16a and b), we do not attempt the comparison between the tangential component earthquake and TT component NGF attenuation curves at 0.25 Hz. RR component amplitudes We also compare the RR component SGF and RSGF FSAs with NGF FSAs at 0.25 and 0.42 Hz (see Supporting Information section RR Amplitudes, Supporting Information Figs S12 and S13). The attenuation curves decay some oscillatory features that are likely an artefact of the Fourier Transform of the near-field term in the solution of elastic wave equation for a radial force and radial velocity response (Aki & Richards 2002). While there is considerable scatter in the NGF amplitudes, the average RR component NGF decay curves at 0.25 and 0.42 Hz indicate at least partial recovery of the near-field term. We recommend that future noise cross-correlation studies, especially ones involving broad-band sensors at short interstation distances, should explore any evidence for near-field terms. Bias in NGF amplitudes Accuracy of relative amplitude information retrieved from ambient noise cross-correlations has been widely discussed and debated in both theoretical and numerical studies (e.g. Cupillard & Capdeville 2010; Tsai 2011; Lawrence et al.2013). Commonly, a wave equation solution of the form $${e^{ - \alpha r}}{J_0}( {\frac{{\omega r}}{c}} )$$ is assumed in which J0 is the zero-order Bessel function of the first kind with c being the frequency-dependent phase velocity. This solution is valid for surface waves that are usually retrieved from noise cross-correlation, and consists of an implicit geometrical spreading factor in the Bessel function (∝r− 0.5 in the far-field) and attenuation in the form of e− αr. For intrinsic attenuation in a homogenous medium, Tsai (2011) showed that retrieval of this solution from coherency measurements is possible only under the conditions of a uniform noise source distribution. The ambient noise source distribution at The Geysers is evidently not uniform as none of our NGFs are symmetric. In case that the noise source distribution is not uniform but is ubiquitous (present both inside and outside an array), then azimuthal averaging of coherency over same-distance station pairs may provide correct relative amplitudes (e.g. Lawrence et al.2013; Zhang & Yang 2013; Weemstra et al.2015). Whatever the source may be, the scatter in the observed NGF amplitudes requires some averaging for meaningful interpretation, which is also true for earthquake ground motions on the same spatial scale (Fig. 16d). Otherwise, to avoid loss of local resolution from averaging, one can resort to more sophisticated techniques such as the C3 method for analysing NGF amplitudes at receivers along a line (Zhang & Yang 2013). In case of far-field noise source distribution, which is expected for continental regions in primary and secondary microseismic passbands (Stehly et al.2006), correct intrinsic attenuation cannot be retrieved even if the illumination is uniform from all directions (Tsai 2011; Weemstra et al.2015). In our study, the amplitude decay shown by the TT component RSGFs for the weak attenuation REF model at 0.72 Hz can be assumed to be equivalent to geometrical spreading, which gives us r− 0.53 for 1λ ≲ r ≲ 13 km (λ ∼ 3.5 km for cL ∼2.5 km s−1 at 0.72 Hz) and r− 1.17 for r ≳ 16 km, which are steeper than r− 0.5 (similar to SGF in Fig. 16c). Nevertheless, we compare our results for the TT component NGFs at 0.72 Hz with various conclusions drawn analytically by Tsai (2011) for different noise source distributions, assuming plane-wave incidence (r ≫ c/ω) and weak intrinsic attenuation (Supporting Information Fig. S14). We do not attempt to correct for deviations from the analytical formulations arising from orientations of components and polarizations of incident waves or effects of near-field distances (Aki 1957; Haney et al.2012). Since we observe a decay that is stronger than the expected G(r), we rule out the possibility of a single point noise source distribution that would have led to little or no decay in amplitudes with distance (Cupillard & Capdeville 2010; table 1 in Tsai 2011). In case of uniform or one-sided far-field noise source distribution, a stronger than G(r) decay is expected, which is observed in our NGF amplitudes (Cupillard & Capdeville 2010; eqs 25 and 29 in Tsai 2011). The apparent attenuation factors for the two distributions are (I0(αr))− 1 and $${( {I_0^2( {\alpha r} ) - L_0^2( {\alpha r} )} )^{ - 0.5}}$$, respectively, where I0 is the modified Bessel function of first kind and order 0 and L0 is the modified Struve function of order 0. The apparent attenuation decays are significantly weaker than the observed NGF amplitude decay. Given the scatter in the NGF amplitudes, it is difficult to constrain α, and it is possible that the decay observed in our NGF amplitudes might just be an artefact of a non-uniform far-field noise source distribution. However, 0.72 Hz is higher than the frequency band of secondary microseisms generated near the coast (Stehly et al.2006) and our analysis of T-[R,Z] components of NGFs didn't reveal any preferred noise-source incidence direction. Another possible scenario includes that all background noise is generated from the anthropological activities in the geothermal field. While the long semi-axis of the geothermal field at r' ∼12 km is not significantly smaller than the attenuation distance 1/α ∼ 1/0.04 ∼ 25 km at 0.72 Hz, using eq. 44 in Tsai (2011) for truncated near-field source distribution (apparent attenuation ∼ $${( {1 - \frac{{{r^2}}}{{4r{'^2}}}} )^{0.5}}$$), we obtain a decay that is still weaker than the observed decay for r < 20 km. This noise source distribution also predicts little or no recovery of coherence for stations outside the geothermal field whereas we obtain robust NGFs for many station pairs located at the edge of the reservoir area (e.g. SRB-DRH in Fig. 12). We are unable to discern any obvious stronger-than-actual amplitude decay at distances ≲ 2λ expected from eq. (1) in which ensemble averaging is done after spectral normalization as opposed to ensemble averaging the cross-spectrum and the amplitude spectra separately prior to the normalization (Tsai 2011; Weemstra et al.2014). It has been suggested that the approach followed in eq. (1) might be helpful if the background noise in highly non-stationary (Weemstra et al.2014). Given the difficulty in analytically estimating coherency for realistic noise source distributions, Lawrence et al. (2013) employed numerical tests and showed that azimuthal averaging of same-distance coherency measurements, similar to averaging of FSAs in our study, can yield correct attenuation estimates under a wide range of noise source distributions. We note that while averaging should be performed within spatial dimensions with slowly and smoothly varying background medium properties (Weemstra et al.2015), we average over the entire reservoir area, which doesn't seem to produce any obvious artefact in the average decay of amplitudes. Among other factors, incoherent noise locally observed at stations contributes to the autocorrelations in the denominator in coherency expression (eq. 1) and may play an important role as the contribution of coherent wavefield weakens with distance (e.g. Tsai 2011; Lawrence et al.2013). In recent studies, Weemstra et al. (2015) have shown that scattering attenuation in a non-dissipative medium, illuminated uniformly from far-field sources, can be correctly recovered from cross-spectrums. At The Geysers geothermal field, given the high temperatures (Lowenstern & Janik 2003) and considerable heterogeneities and fractures present in the subsurface (Lockner et al.1982; Thompson 1989; Gunderson 1991; O'Connell & Johnson 1991; Sammis et al.1992; Elkibbi et al.2005; Jeanne et al.2014), both intrinsic (Romanowicz 1995) and scattering attenuation are expected to be high. The QS values adopted in the REF model, ∼25 to ∼40 in the top ∼4.5 km and ≳60 at greater depths, were estimated using the NetMoment method (Hutchings 2001; Viegas & Hutchings 2011). They are considerably lower than QS values expected from standard VS-QS relationships (Brocher 2008), for example, QS ∼150–300 for VS ∼ 2.0–3.1 km s−1 in the top ∼4.5 km in the REF model, implying stronger attenuation, and possible contribution of scattering attenuation that might be better recovered under more realistic noise-source distributions. Scattering outside the study region can also act to homogenize the noise source illumination, which would aid in the better recovery of amplitudes from NGFs. In methodological aspects, some studies prefer using the same spectral normalization factors for all station pairs to obtain more appropriate relative amplitudes (e.g. Denolle et al.2013; Bowden et al.2015). However, studies applying spectral whitening at multiple stages have been successful in retrieving reliable estimates of anelastic attenuation as well (Handel et al.2016). The scatter in the NGF amplitudes observed at all frequencies necessitates a rigorous analysis of uncertainties and trade-offs in the ground motion attenuation parameters extracted from NGF FSAs, which is beyond the scope of this study. An investigation of ambient noise source distribution at The Geysers along with more NGF amplitude data are also required to properly resolve the individual contributions of possible biases, G(r) and anelastic attenuation to the overall observed decay of NGF amplitudes with distance. Notwithstanding the incompleteness of our analysis, the similarity between the average attenuation curves of the NGF FSAs, earthquake ground motion FSAs and the SGF FSAs predicted from an appropriate velocity model is very compelling and indicates that NGF amplitudes at The Geysers are believable at the distances and frequencies examined here. CONCLUSIONS AND FUTURE WORK When three-component sensors are deployed for long durations, ambient noise cross-correlation techniques provide a robust empirical, nine-component, NGF tensor that can be used in a variety of different applications. The main conclusions of our study are summarized below: We were able to retrieve NGFs in the frequency range (∼0.2–0.9 Hz) for a range of interstation distances from ∼1–30 km (∼0.22 λ–6.5 λ) at The Geysers using a variety of sensors in and around the reservoir area. For many station pairs, the NGFs are found to be similar to the single force displacement SGFs computed using pre-existing and revised 1-D velocity models in terms of waveforms, phase and the relative amplitudes of all components of the tensors, even at distances < λ. We identify both body-wave and surface-wave phases in the NGF waveforms. The direct comparison of NGFs with SGFs helps to evaluate the quality of the retrieved NGFs and the suitability of different 1-D velocity models to various sub-regions of The Geysers. SGFs computed with the faster model VSP0 and the slower model REF preferentially provide better waveform fits to the NGFs for interstation paths across SE and NW Geysers, respectively. We are also able to confirm the results of previous body-wave traveltime tomography studies such as the low velocities in the region around station GSG and higher velocities of the reservoir area relative to the surrounding region. Large anomalous amplitudes on the off-diagonal T-[R,Z] components of NGFs helped us detect sensor misalignments for many stations of The Geysers. We confirm the alignment angles estimated from the analysis of long period teleseismic waveforms by a significant reduction in these anomalous amplitudes upon rotation of the NGF tensors to correct orientations. The comparison between NGFs and SGFs computed using 1-D anisotropic velocity models and a 3-D isotropic velocity model of the reservoir suggests that robust amplitudes on the T-[R,Z] components of NGFs for some longer distance paths likely result from wave propagation effects caused by strong heterogeneity in 3-D structure. The 3-D model that was derived from body-wave traveltimes using an infinite frequency geometrical ray tomography approach (Gritto et al.2013a) can fit low frequency NGF (0.2–0.9 Hz) waveforms remarkably well. We were unable to detect any dominant ambient noise source illumination direction applying ORA to our NGFs. While there is considerable scatter in the TT component NGF FSAs, we find their average decay with distance to be similar to the decay expected from SGF amplitudes and with the decay of tangential component local-earthquake ground-motion amplitudes at the same frequencies (∼0.25–0.72 Hz). The flattening of RR component FSAs at distances ∼9–16 km suggests possible recovery of the near-field term. The similarity of the NGF and the SGF waveforms computed with appropriate velocity models in this study indicates that the full nine-component NGF tensor should be used in waveform tomography studies (e.g. Lee et al.2014) whenever multicomponent stations are available in dense networks. As demonstrated, higher amplitude TT components provide strong reliable constraints on VS structure. ZR and RZ component NGFs are possibly less susceptible to effects of directional noise source incidence compared to the commonly used ZZ component. Broad-band phase delays from NGF waveforms as in the case of the NGFs to station GSG can be potentially useful at short interstation distances. Calibration of starting velocity models using NGF waveforms rather than earthquake waveforms can be desirable in regions like The Geysers, where earthquakes might have non-double couple source mechanisms (Guilhem et al.2014). At The Geysers, possibilities of future work include further refinement of the velocity models, monitoring of temporal changes in the coda of NGFs and characterization of background noise source distribution. Given the excellent distribution of earthquakes and stations within the geothermal field at The Geysers, the 3-D velocity model derived from body-wave tomography is well resolved in the reservoir area; therefore, it is generally successful fitting low-frequency NGF waveforms. However, the model is not well resolved at depths ≳ 4.2 km owing to the shallow focal depths of the earthquakes ( ≲ 4.5 km) or at very shallow depths ( ≲ 1.0 km). Long period Rayleigh waves (∼5–8 s; cR ∼ 2.6–2.9 km s−1 for the GIL7 model) derived from noise cross-correlation at interstation distances (∼40–70 km; r ≳ 3λ) can be used to investigate the velocity structure of the deeper reservoir rocks and the underlying felsite (depth sensitivity ∼λ/2 ∼ 6–12 km). We have demonstrated good results with some of the short period (∼1.0 Hz) sensors of the NCSN that are located outside The Geysers. Retrieval of longer period measurements will also benefit from temporary broad-band sensors deployed at The Geysers (Specht et al.2014) that are expected to yield more stable and higher quality NGFs than the ∼4.5 Hz geophones used in this study. We didn't find any significant energy above 1.0 Hz in our NGFs and the presence of multiple spectral peaks (see Appendix A2) between 1.7–2.5 Hz precluded any meaningful interpretation at these frequencies. Therefore, retrieval of high frequency (≳ 1 Hz) empirical GFs from cross-correlation of coda waves of the numerous local earthquakes in the reservoir area can be attempted (e.g. Campillo & Paul 2003) to investigate shallow velocity structure (cR ∼ 2.2 km s−1 at ∼1.0 Hz for REF model; depth sensitivity ∼λ/2 ∼ 1.1 km). A refined version of the 3-D velocity model can subsequently be applied to higher frequency (up to ∼2.5 Hz) earthquake waveforms and should increase confidence in source inversion studies by facilitating the inclusion of more stations at longer distances (e.g. Guilhem et al.2014). For any future study attempting improvement of velocity models at The Geysers, the results obtained in this study provide a framework for appropriate initial models for inversion. Possible temporal changes in the coda of noise cross-correlations can be investigated to look for possible changes in subsurface VS related to micro/macroseismicity or injection/production activities. Variations of the order of few per cent in VS similar to ∼4–6 per cent variations in VP/VS ratio inferred from earthquake body-wave travel times (Foulger et al.1997; Gritto & Jarpe 2014) should be readily detectable in the coda of NGFs. However, it must be noted that many stations don't return stable or robust NGFs during the entire time period. Characterization of ambient noise source distribution and incidence direction in the frequency range ∼0.2–0.9 Hz at The Geysers and its temporal variation using beamforming analysis, etc. (e.g. Menon et al.2014) is an important task. As discussed previously, ambient noise source distribution has important implications for the decay of NGF amplitudes and amplitudes in the T-[R,Z] components of NGFs. It would be interesting to see if effects on NGFs expected from the inferred noise source distribution (Fichtner 2014) can be reconciled with our observations. Acknowledgements A.N. is indebted to R.B. Herrmann (St. Louis University) and F.-C. Lin (University of Utah) for their help and guidance on many noise cross-correlation and surface wave concepts, and for reviewing the first draft of this paper. We thank an anonymous reviewer, D.C. Bowden (Caltech), M.A. Denolle (Harvard University) and Editor L. Métivier for carefully reviewing this paper, and for their comments and suggestions, which helped to significantly improve this study. We are thankful to V.H. Lai (Caltech), P. Hellweg [Berkeley Seismological Lab (BSL)], H. Rademacher (BSL) and V.C. Tsai (Caltech) for helpful comments and/or discussions. We also thank Ramsey Haught for visiting the stations and checking the sensor orientations. A.N. thanks V. Maupin (University of Oslo) for pointing out the reference for the elastic tensor for crack-induced anisotropy. We are grateful for permission to use the Baribu cluster at BSL for running SW4. This study was supported by National Science Foundation award EAR-1053211 and partially supported by the France-Berkeley Fund 2014-0051. Data used in this study come from Berkeley Digital Seismic Network (BDSN; dx.doi.org/10.7932/BDSN), operated by the University of California Berkeley Seismological Laboratory, Northern California Seismic Network (NCSN), and Lawrence Berkeley National Laboratory (LBNL) Short Period Network at The Geysers, which are all archived at the Northern California Earthquake Data Center (NCEDC; dx.doi.org/10.7932/NCEDC; http://www.ncedc.org, last accessed July 2016). Codes for computing SGFs for 1-D velocity models by FKI and the modal summation method, codes for extracting Love wave group and phase velocities from NGFs and codes for forward modelling the same using 1-D velocity models are all available in the software Computer Programs in Seismology available at http://www.eas.slu.edu/eqc/eqccps.html, last accessed June 2016. SW4 is hosted by the Computational Infrastructure for Geodynamics (CIG) which is supported by the National Science Foundation award NSF-0949446. SW4 can be downloaded from CIG website https://geodynamics.org/cig/software/sw4/, last accessed June 2016. REFERENCES Aki K., 1957. Space and time spectra of stationary stochastic waves, with special reference to microtremors, Bull. Earthq. Res. Inst. Univ. Tokyo , 25, 415– 457. Aki K., Richards P.G., 2002. Quantitative Seismology , 2nd edn, University Science Books, pp. 27– 72. Allis R.G., Shook G.M., 1999. An alternative mechanism for the formation of The Geysers vapor-dominated reservoir, in Proc. 24th Workshop on Geotherm. Reservoir Engg. , Stanford University, California. Atkinson G., 2004. Empirical attenuation of ground-motion spectral amplitudes in Southeastern Canada and the Northeastern United States, Bull. seism. Soc. Am. , 94( 3), 1079– 1095. https://doi.org/10.1785/0120030175 Google Scholar CrossRef Search ADS   Bakun W.H., 1984. Seismic moments, local magnitudes and coda-duration magnitudes for earthquakes in Central California, Bull. seism. Soc. Am. , 74( 2), 439– 458. Bensen G.D., Ritzwoller M.H., Barmin M.P., Levshin A.L., Lin F.-C., Moschetti M.P., Shapiro N.M., Yang Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, Geophys. J. Int. , 169( 3), 1239– 1260. https://doi.org/10.1111/j.1365-246X.2007.03374.x Google Scholar CrossRef Search ADS   Bowden D.C., Tsai V.C., Lin F.-C., 2015. Site amplification, attenuation, and scattering from noise correlation amplitudes across a dense array in Long Beach, CA, Geophys. Res. Lett. , 42( 5), 1360– 1367. https://doi.org/10.1002/2014GL062662 Google Scholar CrossRef Search ADS   Bowman J.R., Kennett B.L.N., 1991. Propagation of Lg waves in the North Australian Craton: influence of crustal velocity gradients, Bull. seism. Soc. Am. , 81( 2), 592– 610. Boyd O.S., Dreger D.S., Lai V.H., Gritto R., 2015. A systematic analysis of seismic moment tensor at The Geysers geothermal field, California, California, Bul. seism. Soc. Am. , 105( 6), 2969– 2986. https://doi.org/10.1785/0120140285 Google Scholar CrossRef Search ADS   Brenguier F., Campillo M., Takeda T., Aoki Y., Shapiro N.M., Briand X., Emoto K., Miyake H., 2014. Mapping pressurized volcanic fluids from induced crustal seismic velocity drops, Science , 345( 6192), 80– 82. https://doi.org/10.1126/science.1254073 Google Scholar CrossRef Search ADS PubMed  Brocher T.M., 2008. Key elements of regional seismic velocity models for long period ground motion simulations, J. Seismol. , 12( 2), 217– 221. https://doi.org/10.1007/s10950-007-9061-3 Google Scholar CrossRef Search ADS   Burger R.W., Somerville P.G., Barker J.S., Herrmann R.B., Helmberger D.V., 1987. The effect of crustal structure on strong ground motion attenuation relations in Eastern North America, Bull. seism. Soc. Am. , 77( 2), 420– 439. Campillo M., Paul A., 2003. Long-range correlations in the diffuse seismic coda, Science , 299( 5606), 547– 549. https://doi.org/10.1126/science.1078551 Google Scholar CrossRef Search ADS PubMed  Cheng C.H., 1993. Crack models for a transversely isotropic medium, J. geophys. Res. , 98( B1), 675– 684. https://doi.org/10.1029/92JB02118 Google Scholar CrossRef Search ADS   Cho K.H., Herrmann R.B., Ammon C.J., Lee K., 2007. Imaging the upper crust of the Korean peninsula by surface-wave tomography, Bull. seism. Soc. Am. , 97( 1B), 198– 207. https://doi.org/10.1785/0120060096 Google Scholar CrossRef Search ADS   Crampin S., 1984. Effective anisotropic elastic constants for wave propagation through cracked solids, Geophys. J. Int. , 76( 1), 135– 145. https://doi.org/10.1111/j.1365-246X.1984.tb05029.x Google Scholar CrossRef Search ADS   Cupillard P., Capdeville Y., 2010. On the amplitude of surface waves obtained by noise correlation and the capability to recover the attenuation: a numerical approach, Geophys. J. Int. , 181, 1687– 1700. Denolle M.A., Dunham E.M., Prieto G.A., Beroza G.C., 2013. Ground motion prediction of realistic earthquake sources using the ambient seismic field, J. geophys. Res. , 118( 5), 2102– 2118. https://doi.org/10.1029/2012JB009603 Google Scholar CrossRef Search ADS   Denolle M.A., Miyake H., Nakagawa S., Hirata N., Beroza G.C., 2014. Long-period seismic amplification in the Kanto Basin from the ambient seismic field, Geophys. Res. Lett. , 41( 7), 2319– 2325. https://doi.org/10.1002/2014GL059425 Google Scholar CrossRef Search ADS   Durand S., Montagner J.P., Roux P., Brenguier F., Nadeau R.M., Ricard Y., 2011. Passive monitoring of anisotropy change associated with the Parkfield 2004 earthquake, Geophys. Res. Lett. , 38( 13), n/a-n/a. https://doi.org/10.1029/2011GL047875 Google Scholar CrossRef Search ADS   Dziewonski A.M., Bloch S., Landisman M., 1969. A technique for the analysis of transient seismic signals, Bull. seism. Soc. Am. , 59, 427– 444. Eberhart-Phillips D., 1986. Three-dimensional velocity structure in northern California Coast Ranges from inversion of local earthquake arrival times, Bull. seism. Soc. Am. , 76( 4), 1025– 1052. Eberhart-Phillips D., Oppenheimer D.H., 1984. Induced seismicity in The Geysers geothermal area, California, J. geophys. Res. , 82( B2), 1191– 1207. https://doi.org/10.1029/JB089iB02p01191 Google Scholar CrossRef Search ADS   Elkibbi M., Rial J.A., 2005. The Geysers geothermal field: results from shear-wave splitting analysis in a fractured reservoir, Geophys. J. Int. , 162( 3), 1024– 1035. https://doi.org/10.1111/j.1365-246X.2005.02698.x Google Scholar CrossRef Search ADS   Elkibbi M., Yang M., Rial J.A., 2005. Crack-induced anisotropy models in The Geysers geothermal field, Geophys. J. Int. , 162( 3), 1036– 1048. https://doi.org/10.1111/j.1365-246X.2005.02697.x Google Scholar CrossRef Search ADS   Fichtner A., 2014. Source and processing effects on noise correlations, Geophys. J. Int. , 197( 3), 1527– 1531. https://doi.org/10.1093/gji/ggu093 Google Scholar CrossRef Search ADS   Foulger G.R., Grant C.C., Ross A., Julian B.R., 1997. Industrially induced changes in earth structure at The Geysers geothermal area, California, Geophys. Res. Lett. , 24( 2), 135– 137. https://doi.org/10.1029/96GL03152 Google Scholar CrossRef Search ADS   Grigoli F., Cesca S., Dahm T., Krieger L., 2012. A complex linear least-squares method to derive relative and absolute orientations of seismic sensors, Geophys. J. Int. , 188( 3), 1243– 1254. https://doi.org/10.1111/j.1365-246X.2011.05316.x Google Scholar CrossRef Search ADS   Gritto R., Yoo S.-H., Jarpe S.P., 2013a. 3D seismic tomography at The Geysers geothermal field, CA, USA, in Proc. Thirty-Eighth Workshop on Geothermal Reservoir Engineering , Stanford University, California, 11–13 February 2013. Gritto R., Yoo S.-H., Jarpe S., 2013b. Seismic Imaging of reservoir structure at The Geysers geothermal reservoir, in AGU Fall Meeting 2013 S33D-2460 , San Francisco, USA, 9–13 December 2013. Gritto R., Jarpe S.P., 2014. Temporal variations of VP/VS-ratio at The Geysers geothermal field, USA, Geothermics , 52, 112– 119. https://doi.org/10.1016/j.geothermics.2014.01.012 Google Scholar CrossRef Search ADS   Guilhem A., Hutchings L., Dreger D.S., Johnson L.R., 2014. Moment tensor inversions of M ∼ 3 earthquakes in the Geysers geothermal fields, California, J. geophys. Res. , 119( 3), 2121– 2137. https://doi.org/10.1002/2013JB010271 Google Scholar CrossRef Search ADS   Gunderson R.P., 1991. Porosity of reservoir greywacke at The Geysers, in Monograph on The Geysers Geothermal Field , Geothermal Resources Council Spl. Rept. No. 17, pp. 89– 93. Handel A., Ohrnberger M., Krüger F., 2016. Extracting near-surface QL between 1–4 Hz from higher-order noise correlations in the Euroseistest area, Greece, Geophys. J. Int. , 207( 2), 655– 666. https://doi.org/10.1093/gji/ggw295 Google Scholar CrossRef Search ADS   Haney M., Mikesell T.D., van Wijk K., Nakahara H., 2012. Extension of the spatial autocorrelation (SPAC) method to mixed-component correlations of surface waves, Geophys. J. Int. , 191( 1), 189– 206. https://doi.org/10.1111/j.1365-246X.2012.05597.x Google Scholar CrossRef Search ADS   Haskell N.A., 1964. Radiation pattern of surface waves from point sources in a multi-layered medium, Bull. seism. Soc. Am. , 54( 1), 377– 393. Hearn B.C.Jr, Donnelly-Nolan J.M, Goff F.E., 1981. The Clear Lake volcanics: Tectonic setting and magma sources, U.S. Geol. Surv. Prof. Pap.  1141, 25– 45. Herrmann R.B., 1973. Some aspects of band-pass filtering of surface waves, Bull. seism. Soc. Am. , 63( 2), 663– 671. Herrmann R.B., 2013a. Computer programs in seismology: An evolving tool for instruction and research, Seismol. Res. Lett. , 84( 6), 1081– 1088. https://doi.org/10.1785/0220110096 Google Scholar CrossRef Search ADS   Herrmann R.B., 2013b. Update to do_mft for the determination of phase velocities from empirical Green's functions from noise cross-correlation, http://www.eas.slu.edu/eqc/eqc_cps/TUTORIAL/EMPIRICAL_GREEN/index.html, last accessed June 2016. Hudson J.A., 1981. Wave speeds and attenuation of elastic waves in material containing cracks, Geophys. J. Int. , 64( 1), 133– 150. https://doi.org/10.1111/j.1365-246X.1981.tb02662.x Google Scholar CrossRef Search ADS   Hudson J.A., 1980. Overall properties of a cracked solid, Math. Proc. Camb. Phil. Soc. , 88( 02), 371– 384. https://doi.org/10.1017/S0305004100057674 Google Scholar CrossRef Search ADS   Hutchings L., 2001. Program NetMoment; simultaneous calculation of moment, source corner frequency, and site specific t* from network recordings, Lawrence Livermore National Laboratory Report  UCRL-ID-135693-REV-1. Jeanne P., Rutqvist J., Hartline C., Garcia J., Dobson P.F., Walter M., 2014. Reservoir structure and properties from geomechanical modeling and microseismicity analyses associated with an enhanced geothermal system at The Geysers, California, Geothermics , 51, 460– 469. https://doi.org/10.1016/j.geothermics.2014.02.003 Google Scholar CrossRef Search ADS   Johnson C.W., Totten E.J., Bürgmann R., 2016. Depth migration of seasonally induced seismicity at The Geysers geothermal field, Geophys. Res. Lett. , 43( 12), 6196– 6204 https://doi.org/10.1002/2016GL069546 Google Scholar CrossRef Search ADS   Julian B.R., Ross A., Foulger G.R., Evans J.R., 1996. Three-dimensional seismic image of a geothermal reservoir: The Geysers, California, Geophys. Res. Lett. , 23( 6), 685– 688. https://doi.org/10.1029/95GL03321 Google Scholar CrossRef Search ADS   Kawakatsu H., Montagner J.-P., Song T.-R.A., 2015. On DLA's η, in The Interdisciplinary Earth: A Volume in Honor of Don L. Anderson , Geol. Soc. Am. Spl. Pap. 514 and Amer. Geophys. Union Spl. Publ. 71, pp. 33– 38, eds Foulger G.R., Lustrino M, King S.D, Geological Society of America. Google Scholar CrossRef Search ADS   Kimman W.P., Trampert J., 2010. Approximations in seismic interferometry and their effects on surface waves, Geophys. J. Int. , 182, 461– 476. Lawrence J.F., Denolle M.A., Seats K.J., Prieto G.A., 2013. A numeric evaluation of attenuation from ambient noise correlation functions, J. geophys. Res. , 118( 12), 6134– 6145. https://doi.org/10.1002/2012JB009513 Google Scholar CrossRef Search ADS   Lee E.-J., Chen P., Jordan T.H., Maechling P.B., Denolle M.A., Beroza G.C., 2014. Full-3-D tomography for crustal structure in Southern California based on the scattering-integral and the adjoint-wavefield methods, J. geophys. Res. , 119( 8), 6421– 6451. https://doi.org/10.1002/2014JB011346 Google Scholar CrossRef Search ADS   Lin F.-C., Moschetti M.P., Ritzwoller M.H., 2008. Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps, Geophys. J. Int. , 173( 1), 281– 298. https://doi.org/10.1111/j.1365-246X.2008.03720.x Google Scholar CrossRef Search ADS   Lin F.-C., Li D., Clayton R.W., Hollis D., 2013. High-resolution 3D shallow crustal structure in Long Beach, California: application of ambient noise tomography on a dense seismic array, Geophysics , 78( 4), Q45– Q56. https://doi.org/10.1190/geo2012-0453.1 Google Scholar CrossRef Search ADS   Lin F.-C., Tsai V.C., Schmandt B., 2014. 3-D crustal structure of the western United States: application of Rayleigh-wave ellipticity extracted from noise cross-correlations, Geophys. J. Int.  198( 2), 656– 670. https://doi.org/10.1093/gji/ggu160 Google Scholar CrossRef Search ADS   Liu X., Kieffer S.W., 2009. Thermodynamics and Mass Transport in Multicomponent, Multiphase H 2 O Systems of Planetary Interest, Annu. Rev. Earth Planet. Sci. , 37( 1), 449– 477. https://doi.org/10.1146/annurev.earth.031208.100109 Google Scholar CrossRef Search ADS   Lockner D.A., Summers R., Moore D., Byerlee J.D., 1982. Laboratory Measurements of Reservoir Rock from the Geysers Geothermal Field, California, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. , 19( 2), 65– 80. https://doi.org/10.1016/0148-9062(82)91632-1 Google Scholar CrossRef Search ADS   Lowenstern J.B., Janik C.J., 2003. The origins of reservoir liquids and vapors from The Geysers geothermal field, California (USA), in Volcanic, Geothermal and Ore-forming Fluids: Rulers and Witnesses of Processes within the Earth , pp. 181– 195, eds Simmons S.F., Graham I., Society of Economic Geologists Special Publication 10. Ludwin R.S., Cagnetti V., Bufe C.G., 1982. Comparison of seismicity in The Geysers geothermal area with the surrounding region, Bull. seism. Soc. Am. , 72( 3), 863– 871. Ma S., Prieto G.A., Beroza G.C., 2008. Testing community velocity models for southern California using the ambient seismic field, Bull. seism. Soc. Am. , 98( 6), 2694– 2714. https://doi.org/10.1785/0120080947 Google Scholar CrossRef Search ADS   Majer E.L., McEvilly T.V., 1979. Seismological investigations at The Geysers geothermal field, Geophysics , 44( 2), 246– 269. https://doi.org/10.1190/1.1440965 Google Scholar CrossRef Search ADS   Majer E.L., Peterson J.E., 2007. The impact of injection on seismicity at The Geysers, California geothermal field, Int. J. Rock Mech. Min. Sci. , 44( 8), 1079– 1090. https://doi.org/10.1016/j.ijrmms.2007.07.023 Google Scholar CrossRef Search ADS   Majer E.L., McEvilly T.V., Eastwood F.S., Myer L.R., 1988. Fracture detection using P-wave and S-wave vertical seismic profiling at The Geysers, Geophysics , 53., 76– 84. https://doi.org/10.1190/1.1442401 Google Scholar CrossRef Search ADS   Maupin V., Park J., 2007. Theory and observations—wave propagation in anisotropic media, in Treatise on Geophysics, Vol. 1: Seismology and the Structure of the Earth , . 289– 321, eds Romanowicz B., Dziewonski A., Elsevier. Menon R., Gerstoft P., Hodgkiss W.S., 2014. On the apparent attenuation in the spatial coherence estimated from seismic arrays, J. geophys. Res. , 119( 4), 3115– 3132. https://doi.org/10.1002/2013JB010835 Google Scholar CrossRef Search ADS   Montagner J.P., Anderson D.L., 1989. Petrological constraints on seismic anisotropy, Phys. Earth planet. Inter. , 54( 1–2), 82– 105. https://doi.org/10.1016/0031-9201(89)90189-1 Google Scholar CrossRef Search ADS   Mossop A., Segall P., 1999. Volume strain within The Geysers geothermal field, J. geophys. Res. , 104( B12), 29 113–29 131. https://doi.org/10.1029/1999JB900284 Google Scholar CrossRef Search ADS   Nakata N., Snieder R., Tsuji T., Larner K., Matsuoka T., 2011. Shear wave imaging from traffic noise using seismic interferometry by cross-coherence, Geophysics , 76( 6), SA97– SA106. https://doi.org/10.1190/geo2010-0188.1 Google Scholar CrossRef Search ADS   Nakata N., Chang J., Lawrence J.F., Boué P., 2015. Body wave extraction and tomography at Long Beach, California, with ambient-noise interferometry, J. geophys. Res. , 120( 2), 1159– 1173. https://doi.org/10.1002/2015JB011870 Google Scholar CrossRef Search ADS   Nishida K., Kawakatsu H., Obara K., 2008. Three-dimensional crustal S wave velocity structure in Japan using microseismic data recorded by Hi-net tiltmeters, J. geophys. Res. , 113( B10), B10302. https://doi.org/10.1029/2007JB005395 Google Scholar CrossRef Search ADS   O'Connell D.R.H., Johnson L.R., 1991. Progressive inversion for hypocenters and P wave and S wave velocity structure: Application to the Geysers, California, Geothermal Field, J. geophys. Res. , 96( B4), 6223– 6236. https://doi.org/10.1029/91JB00154 Google Scholar CrossRef Search ADS   Oppenheimer D.H., 1986. Extensional tectonics at The Geysers geothermal area, California, J. geophys. Res. , 91( B11), 11 463–11 476. https://doi.org/10.1029/JB091iB11p11463 Google Scholar CrossRef Search ADS   Peacock S., Hudson J.A., 1990. Seismic properties of rocks with distributions of small cracks, Geophys. J. Int. , 102( 2), 471– 484. https://doi.org/10.1111/j.1365-246X.1990.tb04479.x Google Scholar CrossRef Search ADS   Petersson N.A., Sjögreen B., 2014a. SW4-v1.1 Users Guide , https://geodynamics.org/cig/software/sw4/, last accessed, December 2015. Petersson N.A., Sjögreen B., 2014b. Super-grid modeling of the elastic wave equation in semi-bounded domains, Commun. Commut. Phys. , 16( 04), 913– 955. https://doi.org/10.4208/cicp.290113.220514a Google Scholar CrossRef Search ADS   Porritt R.W., Allen R.M., Boyarko D.C., Brudzinski M.R., 2011. Investigation of Cascadia segmentation with ambient noise tomography, Earth planet. Sci. Lett. , 309(1–2), 67– 76. https://doi.org/10.1016/j.epsl.2011.06.026 Google Scholar CrossRef Search ADS   Prieto G.A., Beroza G.C., 2008. Earthquake ground motion prediction using the ambient seismic field, Geophys. Res. Lett. , 35( 14), L14304. https://doi.org/10.1029/2008GL034428 Google Scholar CrossRef Search ADS   Prieto G.A., Denolle M.A., Lawrence J.F., Beroza G.C., 2011. On amplitude information carried by the ambient seismic field, C. R. Geosci. , 343( 8–9), 600– 614. https://doi.org/10.1016/j.crte.2011.03.006 Google Scholar CrossRef Search ADS   Romanowicz B., 1995. A global tomographic model of shear attenuation in the upper mantle, J. geophys. Res. , 100( B7), 12 375–12 394. https://doi.org/10.1029/95JB00957 Google Scholar CrossRef Search ADS   Ross A., Foulger G.R., Julian B.R., 1999. Source processes of industrially-induced earthquakes at The Geysers geothermal area, California, Geophysics  64( 4), 1877–1889. Roux P., 2009. Passive seismic imaging with directive ambient noise: application to surface waves and the San Andreas Fault in Parkfield, CA, Geophys. J. Int. , 179( 1), 367– 373. https://doi.org/10.1111/j.1365-246X.2009.04282.x Google Scholar CrossRef Search ADS   Saade M., Montagner J.P., Roux P., Cupillard P., Durand S., Brenguier F., 2015. Influence of seismic anisotropy on the cross correlation tensor: numerical investigations, Geophys. J. Int. , 201( 2), 595– 604. https://doi.org/10.1093/gji/ggu470 Google Scholar CrossRef Search ADS   Sammis C.G., An L., Ershaghi I., 1992. Determining the 3-D fracture structure in the Geysers geothermal reservoir, in Proc. Seventeenth Workshop on Geothermal Reservoir Engineering , Stanford University, California, 29–31 January 1992. Seats K.J., Lawrence J.F., Prieto G.A., 2012. Improved ambient noise correlation functions using Welch's method, Geophys. J. Int. , 188( 2), 513– 523. https://doi.org/10.1111/j.1365-246X.2011.05263.x Google Scholar CrossRef Search ADS   Sens-Schönfelder C., 2008. Synchronizing seismic networks with ambient noise, Geophys. J. Int. , 174( 3), 966– 970. https://doi.org/10.1111/j.1365-246X.2008.03842.x Google Scholar CrossRef Search ADS   Sens-Schönfelder C., Wegler U., 2006. Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia, Geophys. Res. Lett. , 33( 21), L21302. https://doi.org/10.1029/2006GL027797 Google Scholar CrossRef Search ADS   Shapiro N.M., Campillo M., 2004. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise, Geophys. Res. Lett. , 31( 7), L07614. https://doi.org/10.1029/2004GL019491 Google Scholar CrossRef Search ADS   Shapiro N.M., Campillo M., Stehly L., Ritzwoller M.H., 2005. High-resolution surface-wave tomography from ambient seismic noise, Science , 307( 5715), 1615– 1618. https://doi.org/10.1126/science.1108339 Google Scholar CrossRef Search ADS PubMed  Sjögreen B., Petersson N.A., 2012. A fourth order accurate finite difference scheme for the elastic wave equation in second order formulation, J. Sci. Comput. , 52( 1), 17– 48. https://doi.org/10.1007/s10915-011-9531-1 Google Scholar CrossRef Search ADS   Snieder R., 2004. Extracting the Green's function from the correlation of coda waves: a derivation based on stationary phase, Phys. Rev. E , 69( 4), 046610. https://doi.org/10.1103/PhysRevE.69.046610 Google Scholar CrossRef Search ADS   Specht S., Jousset P., Zang A., Gritto R., Bruhn D., 2014. Velocity structure of The Geysers geothermal area (California) from ambient noise cross-correlation, in AGU Fall Meeting 2014 S51A-4425 , San Francisco, USA, 15–19 December 2014. Stark M.A., 1991. Microearthquakes—A tool to track injected water in The Geysers reservoir, in Monograph on The Geysers Geothermal Field , Geothermal Resources Council Spl. Rept. No. 17, pp. 111– 117. Stehly L., Campillo M., Shapiro N.M., 2006. A study of the seismic noise from its long-range correlation properties, J. geophys. Res. , 111( B10), B10306. https://doi.org/10.1029/2005JB004237 Google Scholar CrossRef Search ADS   Stidham C., Antolik M., Dreger D.S., Larsen S., Romanowicz B., 1999. Three-dimensional structure influences on the strong motion wavefield of the 1989 Loma Prieta earthquake, Bull. seism. Soc. Am. , 89, 1184– 1202. Takagi R., Nishida K., Aoki Y., Maeda T., Masuda K., Takeo M., Obara K., Shiomi K., Sato M., Saito K., 2015. A single bit matters: coherent noise of seismic data loggers, Seismol. Res. Lett. , 86( 3), 901– 907. https://doi.org/10.1785/0220150030 Google Scholar CrossRef Search ADS   Taira T., Brenguier F., Kong Q., 2015. Ambient noise-based monitoring of seismic velocity changes associated with the 2014 M w 6.0 South Napa earthquake, Geophys. Res. Lett. , 42( 17), 6997– 7004. https://doi.org/10.1002/2015GL065308 Google Scholar CrossRef Search ADS   Thompson R.C., 1989. Structural stratigraphy and intrusive rocks at The Geysers geothermal field, Trans. Geotherm. Resour. Counc. , 13, 481– 485. Truesdell A.H., Haizlip J.R., Box W.TJr., Amore F.D’, 1991. A geochemical overview of The Geysers geothermal reservoir, in Monograph on The Geysers Geothermal Field , Geothermal Resources Council Spl. Rept. No. 17, 121– 132. Trugman D.T., Shearer P.M., Borsa A.A., Fialko Y., 2016. A comparison of long-term changes in seismicity at The Geysers, Salton Sea, and Coso geothermal fields, J. geophys. Res. , 121( 1), 225– 247. https://doi.org/10.1002/2015JB012510 Google Scholar CrossRef Search ADS   Tsai V.C., 2009. On establishing the accuracy of noise tomography travel-time measurements in a realistic medium, Geophys. J. Int. , 178( 3), 1555– 1564. https://doi.org/10.1111/j.1365-246X.2009.04239.x Google Scholar CrossRef Search ADS   Tsai V.C., 2011. Understanding the amplitudes of noise correlation measurements, J. geophys. Res. , 116( B9), B09311, doi:10.1029/2011JB008483. https://doi.org/10.1029/2011JB008483 Google Scholar CrossRef Search ADS   Tsai V.C., Moschetti M.P., 2010. An explicit relationship between time-domain noise correlation and spatial autocorrelation (SPAC) results, Geophys. J. Int , 182, 454– 460. Google Scholar CrossRef Search ADS   van Wijk K., Mikesell T.D., Schulte-Pelkum V., Stachnik J., 2011. Estimating the Rayleigh-wave impulse response between seismic stations with the cross terms of the Green tensor, Geophys. Res. Lett. , 38( 16), L16301. https://doi.org/10.1029/2011GL047442 Google Scholar CrossRef Search ADS   Viegas G., Hutchings L., 2011. Characterization of Induced Seismicity near an Injection Well at the Northwest Geysers Geothermal Field, California, Trans. Geotherm. Resour. Counc. , 35, 1773– 1780. Wang C.Y., Herrmann R.B., 1980. A numerical study of P-, SV- and SH-wave generation in a plane layered medium, Bull. seism. Soc. Am. , 70( 4), 1015– 1036. Wapenaar K., Fokkema J., 2006. Green's function representations for seismic interferometry, Geophysics  71( 4), SI33– SI46, doi:10.1190/1.2213955. Google Scholar CrossRef Search ADS   Weemstra C., Westra W., Snieder R., Boschi L., 2014. On estimating attenuation from the amplitude of the spectrally whitened ambient seismic field, Geophys. J. Int. , 197( 3), 1770– 1788. https://doi.org/10.1093/gji/ggu088 Google Scholar CrossRef Search ADS   Weemstra C., Snieder R., Boschi L., 2015. On the estimation of attenuation from the ambient seismic field: inferences from distributions of isotropic point scatterers, Geophys. J. Int. , 203( 2), 1054– 1071. https://doi.org/10.1093/gji/ggv311 Google Scholar CrossRef Search ADS   Wegler U., Nakahara H., Sens-Schonfelder C., Korn M., Shiomi K., 2009. Sudden drop of seismic velocity after the 2004 M w 6.6 mid-Niigata earthquake, Japan, observed with Passive Image Interferometry, J. geophys. Res. , 114( B6), B06305. https://doi.org/10.1029/2008JB005869 Google Scholar CrossRef Search ADS   Wielandt E., Forbriger T., 1999. Near-field seismic displacement and tilt associated with the explosive activity of Stromboli, Ann. Geofis. , 42( 3), 407– 416. Yao H., van der Hilst R.D., de Hoop M.V., 2006. Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis—I. Phase velocity maps, Geophys. J. Int. , 166, 732– 744. https://doi.org/10.1111/j.1365-246X.2006.03028.x Google Scholar CrossRef Search ADS   Yokoi T., Margaryan S., 2008. Consistency of the spatial autocorrelation method with seismic interferometry and its consequence, Geophys. Prospect. , 56, 435– 451. https:/