Fault zone reverberations from cross-correlations of earthquake waveforms and seismic noise

SUMMARY Seismic waveﬁelds interact with low-velocity fault damage zones. Waveforms of ballistic fault zone head waves, trapped waves, reﬂected waves and signatures of trapped noise can provide important information on structural and mechanical fault zone properties. Here we extend the class of observable fault zone waves and reconstruct in-fault reverberations or multiples in a strike-slip faulting environment. Manifestations of the reverberations are signiﬁcant, consistent wave fronts in the coda of cross-correlation functions that are obtained from scattered earthquake waveforms and seismic noise recorded by a linear fault zone array. The physical reconstruction of Green’s functions is evident from the high similarity between the signals obtained from the two different scattered waveﬁelds. Modal partitioning of the reverberation waveﬁeld can be tuned using different data normalization techniques. The results imply that fault zones create their own ambiance, and that the here reconstructed reverberations are a key seismic signature of wear zones. Using synthetic waveform modelling we show that reverberations can be used for the imaging of structural units by estimating the location, extend and magnitude of lateral velocity contrasts. The robust reconstruction of the reverberations from noise records suggests the possibility to resolve the response of the damage zone material to various external and internal loading mechanisms.


I N T RO D U C T I O N
The elastic Earth's response or Green's function includes seismic phases associated with main boundaries characterized by significant impedance contrasts.Examples on the global scale include the free surface with Rayleigh waves, the crust-mantle boundary with Moho reflections and refracted head waves, and core-mantle boundary reflections.The associated phases in a seismogram provide constraints on the location of these discontinuities, and therefore on the structure of the Earth.Multiple reflections, also referred to as multiples or reverberations, form characteristic arrival time patterns that are governed by the properties of the waveguide which is typically associated with horizontal layering (Fukao et al. 1983;Verschuur et al. 1992).Vertically oriented waveguides constitute mechanically important zones of strike-slip earthquake faults.The waveguiding wear zones evolve during the slip accumulation process and consist of comminuted and fractured material with lower rigidity and seismic velocity compared to the more competent surrounding rocks (Ben-Zion & Sammis 2003, and references therein).Here we reconstruct reverberations inside a low-velocity damage zone in a strike-slip faulting environment (Fig. 1) and investigate the coupling of both earthquake wavefields and the ambient field with the fault structure.The passive approach establishes new observables that can help to improve fault zone characterization in environments that sustain multiple in-fault reflections.
Ballistic fault zone waves propagate along systematic seismic velocity contrasts across and within fault zones.The term ballistic refers to the general computability of ray paths and waveforms, in contrast to the statistical description of energy propagation associated with multiple scattered or diffuse wavefields.The class of fault zone waves includes, first, head waves that are refracted along velocity contrasts resulting from the cumulative displacement across the slip interface.Second, the interference pattern of critically reflected energy propagating along and within sufficiently persistent and uniform low-velocity damage zones is referred to as trapped or guided waves (Ben-Zion & Aki 1990;Li et al. 1990Li et al. , 1994;;Li & Leary 1990;Ben-Zion 1998).Third, body wave reflections at fault zone boundaries form typical P-and S-wave moveout patterns (Li et al. 2007;Yang et al. 2011).The sensitivity of ballistic and scattered waves to velocity variations alike leads to the noise wavefield that is trapped in the cavity/waveguide-like structure of the damage low-velocity zone (Hillers et al. 2014).Manifestations of trapped noise are an increased randomization of propagation directions through in-fault scattering and subcritical reflections at the fault zone boundaries, which leads to a reduction of wavefield coherence across the region of the velocity gradient.The The configuration represents a ray pattern that differs from the associated traveltime pattern (if the along-strike direction is replaced by traveltime as in Fig. 10a).characteristics of the fault zone wavefields facilitate high-resolution imaging of the waveguide within which they propagate (Lewis et al. 2005;Hillers et al. 2014).This is essential to constrain fault zone dynamics, because the mechanical and rheological properties of the wear zone can have strong effects on co-and post-seismic processes such as rupture propagation regimes, ground motion patterns and healing (e.g.Lewis & Ben-Zion 2010;Yang et al. 2011, and references therein).
The Earth's response between two stations can be estimated from the cross-correlation function (also briefly referred to as correlation) of scattered wavefield records.The wavefield can consist of earthquake waveforms (Campillo & Paul 2003), or it can be the ambient seismic wavefield or noise (Shapiro & Campillo 2004).For high-fidelity estimates of empirical Green's functions the energy propagation in the original wavefield is ideally described by a diffusion process approaching equipartition (Hennino et al. 2001;Lobkis & Weaver 2001;Paul et al. 2005;Campillo 2006).In the multiple scattering regime, this refers to the equipartition of energy in the phase space associated with uncorrelated waves.Properties of seismological wavefields are generally not compatible with this limit regime, but source-and/or time-averaging promote the convergence towards estimates of sufficient accuracy.The robustness of the passive Green's function retrieval is documented by numerous surface wave tomographies (Sabra et al. 2005;Shapiro et al. 2005) and by the reconstruction of body waves including the abovementioned reflections at the Moho (Zhan et al. 2010;Poli et al. 2012a) and the core-mantle boundary (Boué et al. 2013;Lin et al. 2013b).
Detailed images of the hierarchical fault zone architecture (Ben-Zion 2008) are key for an improved understanding of fault dynamics and for progress in seismic hazard assessment.Considering this importance and the many noise-based monitoring studies reporting broadly distributed velocity changes induced by moderate and large earthquakes (e.g.Wegler & Sens-Schönfelder 2007;Brenguier et al. 2008;Hobiger et al. 2012), the number of applications based on the passive Green's function retrieval in fault zone environments appears somewhat limited.To date, investigations of fault zone processes include noise-based imaging studies (Roux 2009;Lin et al. 2013a;Nakata et al. 2015;Zigone et al. 2015), wavefield analyses (Hillers et al. 2014;Zhang & Gerstoft 2014), ground motion estimates (Denolle et al. 2013) and damage zone monitoring using earthquake waveforms (Roux & Ben-Zion 2014).
Extending the scope of the passive approach for the study of fault zone properties we here reconstruct multiples within a strike-slip fault damage zone from cross-correlation of scattered wavefields.We analyse the noise and regional earthquake signals recorded by a linear seismic array installed across the Calico fault, southern California (Section 2).We demonstrate the effect of signal processing choices on the properties of the reconstructed in-fault reverberations (Section 3).The reverberation amplitudes that emerge in the early part of the coda can be tuned relative to the amplitudes of the direct surface wave.A very high similarity (correlation coefficient 0.95) between reverberation waveforms obtained from noise and earthquake records confirms that the observed reverberations are part of the Green's function associated with significant velocity discontinuities.The partition of body and surface wave energy in the reverberation wavefield varies with position relative to the lowvelocity zone and with propagation or lapse time, and depends on original wavefield properties and the applied processing steps (Section 4).We constrain the fault zone velocity structure (Section 5) by reproducing relevant moveout features of the reverberation patterns with synthetic waveforms.The values for the effective fault zone width (∼1 km) and velocity contrast (∼50 per cent) agree with previously obtained estimates and demonstrate that multiples extend the range of fault zone waves that can be used for fault imaging.

Data
We focus on continuous vertical and radial component records of a dense line array consisting of 2-Hz Sercel L-22 sensors across the Calico fault in the Eastern California Shear Zone (Fig. 2a).This is a subset of the array data analysed by Cochran et al. (2009), Yang et al. (2011), Zhang & Gerstoft (2014) and Hillers et al. (2014) for the investigation of fault zone-wavefield interactions.To distinguish potential reverberations excited by earthquakes from those associated with the background wavefield we separate earthquake waveforms and noise.Noise refers to the stationary ambient wavefield that cannot be associated with a specific source, in contrast to transients that are typically identified by their larger amplitudes.We use the term noise for the in-fault and for the fault zone surrounding fields, and in-fault refers to the approximately 1 km wide region (Fig. 2b) around the largest body and surface wave velocity reduction (Cochran et al. 2009;Hillers et al. 2014).The position x = 0 km corresponds to the mapped fault surface trace.
Properties of the ambient wavefield in the Calico fault zone environment depend strongly on frequency and the position relative to the low-velocity zone.Surface wave energy between 0.25 Hz and 2 Hz is pre-dominantly excited west to the site, and propagates approximately perpendicular across the fault.The 1-2 km wide compliant damage zone (Fialko et al. 2002;Cochran et al. 2009;Yang et al. 2011) is relatively transparent for the low frequency wavefield.Around 0.5 Hz ballistic and scattered waves interact with the lowvelocity channel, and above 1 Hz wavefields in-and outside this zone differ significantly (Hillers et al. 2014).The colocated 2-D array (Zhang et al. 2009) can further be used to beamform 1-Hz body waves that are excited off the coast.In summary, the high-frequency wavefield consists of multiple scattered and refracted components, and different data processing and analysis techniques can be used to enhance or mute their respective signatures.To separate regional earthquake signals from noise we search the Southern California Earthquake Data Center (SCEDC 2013) for hypocentre information in the period between 2006 August 25 and November 26.A first earthquake waveform list defines windows of duration w + that begin at the minute t 0 in which M > 3 events occurred.We choose events within a radius R of 250 km around the station located on the mapped fault trace (Fig. 2c).This list contains 25 windows associated with 25 earthquakes that occurred on 19 different days.A second earthquake exclusion list defines long noise windows from these 19 days.The windows are punctuated by short gaps of duration w − = 15 min that follow M > 2 earthquakes occurring up to 400 km away.

Processing
Pre-processing or data conditioning (e.g.Bensen et al. 2007;Groos et al. 2012, and references therein) determines the quality of noise correlation functions and are important for the mitigation of biasing effects associated with non-stationary wavefield fluctuations.Processing steps and parameters have to be adapted to the spatial, temporal and spectral properties of the target problem.The spectral amplitude of a broad-band noise record is typically normalized (whitening) to reduce the effect of inhomogeneously distributed power inherited from frequency dependent excitation, attenuation and amplification.Time domain recipes to decrease the influence of transients include averaging moving-window strategies, threshold clipping, 1-bit normalization and the complete removal of windows containing large amplitudes.Effects of normalization and clipping deserve special attention if amplitude estimates are inverted for attenuation properties (e.g.Cupillard et al. 2011;Prieto et al. 2011).Smart choices of data window lengths can similarly decrease the influence of non-stationary episodes in the resulting correlation stack (Poli et al. 2012b;Seats et al. 2012;Stehly & Hanasoge 2012).Notwithstanding this large parameter space the phase properties of empirical Green's functions are generally well recovered from noise correlations, that is, a reliable reconstruction of the surface wave usually emerges under a broad range of conditions.These considerations apply also to the pre-processing of earthquake waveforms (Campillo & Paul 2003;Paul et al. 2005;Roux & Ben-Zion 2014).The multiply scattered coda wavefields (Fig. 3) govern the faster convergence of the associated correlation functions compared to noise correlations (Gouédard et al. 2008;Chaput et al. 2015).
We focus on the effect of different processing strategies on the reconstruction of in-fault reverberations.These waveforms emerge in the early coda part of the correlations up to τ ≈ 10 T , where τ denotes the lapse time in the correlation function and T is the mean period of a frequency band.Their significance is deduced from the high coherency across the dense line array.Part of all processing sequences (Table 1) is the removal of the instrument response, detrending, highpass filtering, glitch removal and tapering.The reference processing P 0 used by Hillers et al. (2014) includes whitening in the target frequency band.It is fixed to 1-2 Hz throughout this work, following the observation that the trapped noise signature is best expressed in this band.The time domain amplitude is clipped at 3.5 times the standard deviation of the amplitude distribution in each two-hour processing window (threshold clipping).Signatures of trapped noise were similarly obtained using a broad whitening passband, threshold clipping, correlation and subsequent target band filtering (Hillers et al. 2014).This is the most efficient approach for the study of dispersion, and is therefore commonly applied in surface wave imaging.
Here we deviate from this approach.For processing P 1 , records are bandpass filtered in a broad frequency band between f lo and (Highpass filter >0.01 Hz) Whiten 1-2 Hz 1-bit normalization f hi , for example, between 0.03 and 5 Hz (Fig. 3), 1-bit normalized, correlated and filtered in the 1-2 Hz target frequency range.Finally, P 2 processing means records are highpass filtered, whitened between 1-2 Hz, 1-bit normalized and correlated (Table 1).P 1 stands out because of the omitted spectral normalization; P 0 and P 2 differ only in the time domain normalization.All methods are applied to noise data.The decaying envelope of earthquake coda prohibits the application of P 0 including threshold clipping to earthquake waveforms.Therefore, 1-bit clipping is used in P 1 and P 2 .
When we process the noise windows that are chosen from the earthquake catalogue-based exclusion-list, we perform additional statistical tests prior to whitening and/or time-domain normalization to remove spurious large-amplitude transients (Poli et al. 2012b) that are not associated with earthquakes listed in the catalogue.In the following we discuss features of correlation coda wavefields that are obtained with these processing schemes to demonstrate the tuning effect of the different choices.After a brief summary of trapped noise features that are obtained with P 0 , we discuss the very different results that emerge in response to the P 1 -processing of noise and earthquake data.In contrast, noise and earthquake correlation patterns obtained with P 2 processing are much more similar.Both correlation wavefields exhibit strong signals that correspond to infault reverberations, and contributions to the converging correlation stacks are systematically investigated.

Response to P 0
We summarize manifestations of trapped noise in the correlation wavefield obtained from P 0 -processed noise records (Hillers et al. 2014).In-fault 1-2 Hz vertical-vertical (ZZ) noise correlation patterns (Fig. 4a) are characterized by a far-field surface wave that collapses into the large amplitude near-field focal spot (Catheline et al. 2008).The larger surface wave τ -symmetry of in-fault correlations compared to correlations between station pairs outside the fault or across the fault boundaries indicates that the original wavefield is trapped within the low-velocity waveguide (Hillers et al. 2014).Amplitude patterns of correlations between stations located inside the low-velocity channel exhibit coherent wave fronts behind the surface waves (considering absolute lapse times).These wave fronts likely indicate multiples or reverberations, but their amplitudes are small compared to the dominating surface wave.In Figs 4(a)-(d), the correlation reference station S ref is located at x = −0.5 km at the western edge of the low-velocity zone (Figs 2b and c).This representation is chosen to highlight the relevant in-fault features.The loss in surface wave coherence across the velocity gradient can be better visualized using different positions of S ref (Hillers et al. 2014).

Response to P 1
The correlation pattern constructed from 19 d of earthquake-free noise processed with P 1 (Fig. 4b) exhibits an arrival structure that is similar to move-out patterns of low-frequency waves for which the fault environment is relatively transparent (Hillers et al. 2014).The impinging field is generated Southwest of the site, and waves propagate approximately perpendicular across the fault, that is, parallel to the line array.Hillers et al. (2014) showed that these components also exist at 1-2 Hz.Omitting spectral whitening (P 1 , Fig. 4b) enhances those contributions in the correlations that are associated with the directional incident wavefield.In contrast, signatures of the trapped, more isotropic in-fault wavefield emerge using whitening and threshold clipping (P 0 ).Comparison of P 0 , P 1 and P 2 results presented below suggests that the strong differences in the associated patterns are controlled by spectral normalization, and not by the replacement of threshold clipping by 1-bit clipping.It suggests that the spectral amplitude of the incident SW→NE propagating component is relatively large and thus dominates the correlation wavefield if whitening is not applied.
Next we consider short noise windows that are recorded immediately before and after the 25 earthquake signals.We correlate and stack data from two sets.Each consists of 25 5 min long noise windows.In the first (second) set the 5 min long windows begin 6 min before (5 min after) t 0 .The correlation wavefields resulting from both sets are very similar.Fig. 4(c) shows the pattern constructed from the second set.The wavefield is characterized by the very same wave that dominates the discussed 19-d stack of noise correlations (Fig. 4b).Again, this wave traverses the fault in the fault-normal direction.Only the signal-to-noise ratio is smaller compared to the results from the longer averaging times (Sabra et al. 2005;Larose et al. 2007).The similarity between the patterns based on the 25 short post-earthquake coda windows and those based on 19 d of noise implies that energy released by the earthquake is dissipated after 5 min.
This inference is supported by the comparison of the P 1 noisebased patterns (Figs 4a-c) to correlations constructed from earthquake records (Fig. 4d) of the complete wavefield (direct and scattered arrivals).The in-fault amplitude pattern obtained from a stack of 25 earthquake correlations (M > 3, R < 250 km, w + = 5 min; Fig. 4d) also exhibits the focal spot obtained with P 0 noise processing.However, the amplitude symmetry of the surface wave at negative and positive lag times is significantly lower compared to the P 0 noise case in Fig. 4(a).Larger amplitudes emerge at lag times that correspond to SW→NE propagation.This asymmetry emerges from the prevailing flux (Paul et al. 2005) associated with the inhomogeneous hypocentre distribution (Fig. 2a).The good signalto-noise ratio of the arrivals in the correlations of P 1 processed earthquake signals is associated with a relatively short convergence time (25 × 5 min), indicating that scattered earthquake wavefields carry higher coherent energy compared to the background noise.
Importantly, the in-fault correlation coda pattern shows significant amplitudes following the direct far-field surface wave (Fig. 4d).Especially for S ref located at |x| < 200 m, these later arrivals exhibit a typical criss-crossing pattern (shown below) associated with multiple reflected waves within the low-velocity zone.The back-and-forth propagation of energy along the fault-normal x-direction governs the observed high-amplitude symmetry of the Figure 4. ZZ cross-correlation wavefields in the 1-2 Hz range obtained with P 0 and P 1 noise and earthquake processing.Amplitudes are scaled by the peak value in each panel.(a) Correlation wavefield obtained with standard noise processing P 0 .It exhibits a far-field surface wave with relatively symmetric amplitudes at negative and positive τ , and a large-amplitude focal spot.Small-amplitude multiples or reverberations emerge in correlations between stations located inside the low-velocity zone.(b) Processing the same data with the P 1 scheme emphasizes the constituent of the incident ambient seismic field that propagates approximately perpendicular across the fault.(c) The same pattern emerges with lower signal-to-noise ratio from 25 5 min long noise windows that begin 5 min after the earthquake origin times.(d) The correlation pattern obtained from 5 min long earthquake waveforms has more similarities with panel (a) than panel (b) or panel (c); note the increased reverberation amplitudes.reflected phases at negative and positive lapse times.The relative amplitudes of these multiples are considerably larger compared to those from the P 0 processed noise (Fig. 4a).We find that correlating 3, 5 or 10 min long windows, where all windows begin at t 0 , has a small effect on the reverberation amplitudes, although the windowing affects the correlation similarity with noise-based correlations (Section 3.4).It indicates that the reconstructed multiple reflections are primarily controlled by early arrivals, even though 1-bit clipping kills the relative amplitude information.
The quality of the reconstructed reflections depends on the lowfrequency limit of the initial bandpass filter (Table 1).We found that increasing f lo from 0.03 Hz to 0.5 Hz degrades these signals.This can be understood in the context of Stehly & Hanasoge (2012), who showed that the level of random noise added to seismic records prior to 1-bit normalization can affect the partitioning of the reconstructed wave types and symmetry properties of the correlation functions.Contributions from the noise at frequencies outside the 1-2 Hz target band correspond in this view to the added random noise, without, however, making the assumption or implying that long period noise is random.An increased f lo threshold is hence synonymous with reducing the amplitudes of the added noise (Fig. 3b).This is consistent with the observation that spectral noise amplitudes generally decrease from the lower frequencies around the secondary microseisms peak towards our target high-frequency range (Peterson 1993).

Response to P 2
In the previous sections we demonstrated the tuning effect of the different processing schemes P 0 and P 1 on the composition of the correlation wavefields.If now whitening is applied and 1-bit normalization is retained as for P 2 , correlation functions between S ref at x = −0.5 km and stations across the zone of reduced velocities (Fig. 5a) exhibit a nearly amplitude-symmetric surface wave pattern-a signature of trapped noise, since this indicates the increased randomization of propagation directions in the fault zone waveguide (Hillers et al. 2014).As for Fig. 4, S ref at x = −0.5 km is located inside the trapping structure to best visualize the important components.Note that the location of sensors at x > 1 km is not in line with the array (Fig. 2c), which obscures the distance dependent variability compared to the change across the contrast around x = −1 km.The amplitude symmetry in Fig. 5(a) is in strong contrast to the associated pattern shown in Fig. 4(b) that is dominated by the fault-normal propagating component of the incident wavefield.Whitening tends to equilibrate the contribution of wavefield constituents characterized by more isotropic propagation directions.The result is a better reconstruction of the medium response through the reduction of directivity effects associated with unevenly distributed sources.
The most important result of this study is the reconstruction of clean, significant wave fronts that emerge in the coda of the correlations between S ref and stations around x = 0 km.These arrivals can be clearly identified at least up to |τ | = 8 s.The amplitude of these reverberations increases relative to the direct surface wave, that is, the excitation of reverberations is stronger if the virtual source S ref moves from the edge (Fig. 5a) towards the centre of the low-velocity zone (Fig. 5b).This behaviour is similar to the trend observed with P 1 processed earthquake data (Fig. 4d).Recall that P 0 (whitening, threshold clipping) and P 2 (whitening, 1-bit clipping) differ just in the time domain normalization (Table 1).Whereas threshold clipping is sufficient to reconstruct signatures of the isotropic trapped noise wavefield (Fig. 4a), the nonlinear 1-bit clipping enhances reverberations in correlations obtained from noise (Fig. 5a).This trend is again stronger for cases in which S ref is located at small |x|.
The effect of whitening on the reverberation pattern constructed from earthquake waveforms is less striking but still significant.This is inferred by comparing the P 1 and P 2 earthquake based amplitude patterns in Figs 4(d), 5(c) and (d), where clean wave fronts associated with coherent reflections are also visible at least up to |τ | = 8 s.In conclusion, the strongest signature of the reverberations can be extracted from both wavefields with the P 2 processing (targetband whitening, 1-bit clipping).We emphasize the high similarity of the noise and the earthquake correlation results: The patterns in Figs 5(b) and (c) exhibit the same strong multiple reflections.The almost perfect agreement between the correlation waveforms obtained from the two independent scattered wavefields is a strong indicator that both reconstructions constitute physical estimates of the Green's functions.Again, this statement is not limited to the reconstruction of the direct arrival, but includes explicitly the highly τ -symmetric correlation coda arrivals.These arrivals are thus definitely associated with the velocity variations of the medium, in contrast to the often poorly reconstructed coda waves in long-period correlation functions.The high similarity highlights the efficiency of the applied normalization techniques to homogenize the two original wavefields, which are characterized by different intensity distributions.Despite the general similarity secondary differences exist.Whereas the noise-based patterns appear more regular, earthquake-based cross-correlation wavefields are characterized by higher amplitudes at later lapse times.Earthquake-based patterns show a sensitivity to the correlated earthquake coda window (Figs 5c and d).We attribute these differences to incomplete averaging.

Convergence
Considering the reverberations that are reconstructed from earthquake waveforms: How do these signals emerge compared to the coherency build-up in noise correlations?It is evident from Figs 6(a) and (b) that waveforms of individual correlation functions obtained from P 2 processed earthquake seismograms are highly variable.At |τ | > 1 s, there are no coherent features from one trace to the next, and the wavefield clearly lacks a consistent arrival structure.Additional manifestations of incomplete averaging are fluctuations in the otherwise well reconstructed main arrival around τ = 0 s (Roux & Ben-Zion 2014).However, the similarity of the coda waveforms between the stack composed from 25 earthquake-based correlations and the 19 d noise-correlation stack (Fig. 6c, top, and Fig. 7) is remarkable.The similarity is expressed by an average 0.95 correlation coefficient between earthquake-and noise-based cross-correlation waveforms in a window between τ = 1 s and τ = 6 s.These values are obtained from waveforms that are 'τ -averaged' over negative and positive lag times.Such high values are found for station pairs located within a few hundred metres around x = 0 km.Average similarity values drop to the 0.6 level at locations x < −1 km.This implies that the strong lateral velocity contrast controls the rapid convergence between the noise-based and just 25 earthquake-based correlation functions.This argument rests on the interpretation of the noise-correlation stack as a sufficiently converged 'reference' estimate, which follows from the stability of the daily correlations (Fig. 6c).Yet here, too, variations in the signal-to-noise ratio and timing of the multiples indicate fluctuations in the original noise wavefield.It is understood that the remaining difference between noise-and earthquake-based reflection patterns are governed by incomplete averaging.These results confirm that averaging contributions from (many) different events is imperative to obtain sufficiently converged Green's function estimates (Campillo & Paul 2003;Paul et al. 2005;Gouédard et al. 2008).Again, Figs 6(a) and (b) suggest that full convergence is not achieved for a small number of events or with short earthquake data windows.
Correlation waveforms constructed from different earthquake coda windows differ substantially.The position and length of the data window can be thus tuning parameters for improved Green's function estimates (Chaput et al. 2015).For in-fault station pairs we find that correlations constructed from early windows (1-3 min after Figure 7. Lapse time or τ -averaged correlation coda from P 1 and P 2 processed earthquake waveforms and noise.Panels in the top (bottom) row show ZZ (ZR) correlations.The images focus on the 2 km wide region around the mapped fault trace.In the right column, we plot seismograms from the middle column in grey to illustrate the waveform similarity.t 0 ) are more similar to the complete earthquake waveform correlation (0-5 min) compared to correlations constructed from late windows (3-5 min).Supporting our conclusion from Section 3.2, this implies that the properties of earthquake wavefields at early lapse times-which include the direct arrivals-dominate over late coda, that is, the coherency and directivity associated with a more ballistic regime dominate over the diffusive part.In addition, the frequency content changes as a function of the propagation time.Average spectral amplitudes of windowed (τ = −20 s to τ = 20 s) correlations constructed from the early 1-3 min window are significantly larger between 1.2 and 2 Hz compared to spectral amplitudes of correlations obtained from the late 3-5 min window.Green's function estimates reconstructed from the late window are thus depleted at higher frequencies.This is explained by faster high-frequency energy dissipation, visualized by reduced earthquake coda amplitudes after t 0 + 3 min (Fig. 3b).The complete earthquake waveform correlation stack has the highest similarity with the noise-stack and has therefore been used in our analysis.
We investigated factors that can potentially control the similarity between windowed correlation waveforms associated with a single earthquake and the 25-stack.We considered relations between similarity and magnitude, distance, depth, azimuth, correlation peak amplitude and a proxy for earthquake coda energy.The latter is estimated from the envelope-integral of the narrow band filtered velocity seismograms in a 5 min window after t 0 .The results do not indicate any dependence, except for a linear trend between cor-relation coda similarity and correlation peak amplitude, which are both markers of the signal-to-noise ratio.The lack of systematic interrelations implies that other factors such as spatially variable heterogeneity govern the scattered wavefield that excites reverberations inside the Calico fault waveguide.

A N AT O M Y O F T H E R E V E R B E R AT I O N S
We estimate the phase φ between ZZ and ZR correlation coda waves (Fig. 7) to study the polarization state of the wavefield.So far, we focused on the ZZ component of the 3 × 3 correlation tensor, which is the response in the vertical direction to a vertical force.The ZR component captures the empirical response in the radial direction to a vertical force.Radial is here the fault-normal direction x along the instrumented line.Particle motion in the ZZ-ZR plane is indicative of the wave polarization and can thus help to distinguish between body waves and Rayleigh surface wave components (Tanimoto et al. 2006;Landès et al. 2010).We consider the φ(τ ) evolution associated with two different station groups, i.e. we compare phase estimates from all-to-all correlations between the stations in each group.One group consists of five stations that are located around the mapped fault trace at x = 0 km ('centre-5'; Fig. 2).The other group consists of six stations ('outer-6'), where three each are located on both sides to the centre-5 group.We find different φ(τ ) variations for these groups and discuss the systematic lapse time dependence of φ distributions obtained in three 2 s long windows between τ = 1.5 s and τ = 7.5 s (Fig. 8).We neglect estimates in the near-field because of the longitudinal polarization of the focal spot (Catheline et al. 2008).In this analysis we do not account for the systematic moveouts of the direct surface wave (Fig. 7), but we discuss this in detail in Section 5.For each 2 s window we estimate the phase in a shorter 1 s window (±0.5 s) advancing in 0.1 s increments; φ is estimated in the range 0 • -180 • .We apply this to earthquake correlations obtained with P 1 (recall that P 1 noise processing does not reconstruct reverberations), and to earthquake and noise correlations obtained with P 2 (Fig. 7).
We find one main and four secondary trends in the resulting histograms displayed in Fig. 8.The interpretation is based on the phase discrimination for Rayleigh surface waves (90 • or 45 • < φ < 135 • ) and body waves (0 • and 180 • or φ < 45 • and φ > 135 • ) (Tanimoto et al. 2006).The receiver geometry captures energy travelling perpendicular to the fault strike that is polarized in a plane defined by the vertical and radial direction.Amplitudes in the ZZ correlation coda are thus associated with vertically polarized shear waves (SV) and ZR motion is controlled by P-waves propagating along the line.Our phase estimates distinguish between surface and body waves, but a slowness analysis would be required to separate P-from SV-waves.The large ZZ amplitudes compared to the smaller ZR amplitudes in the patterns in Fig. 7 and the high apparent velocities compared to the direct surface wave indicate that SV-waves dominate the reverberation wavefield in the region around x = 0 km.The patterns highlight again the almost perfect agreement between earthquake coda-and noise-based Green's function estimates from which we obtained the 95 per cent similarity.
The main observation is that the properties of the wavefields evolve with lapse time.The compositions at early lapse times differ depending on the original wavefield and on the pre-processing scheme.However, scattering, reflections and the associated mode conversions homogenize the wavefields and lead to the uniform distributions observed at later lapse times (hatched histograms in Fig. 8).The four secondary results are: (1) Body waves dominate all center-5 correlation coda in the central window (Figs 8a-c, open/white histograms).( 2) The centre-5 correlation group based on P 2 processed earthquake data is dominated by surface waves in the first window (Fig. 8b, grey versus open histograms).
(3) The central window body wave dominance (result 1) is for the P 2 processed earthquake waveforms also observed in the outer-6 group (Figs 8b and e, open histograms).( 4) For the P 1 processed earthquake data, the dominating wave type in the first window differs between the two station groups.Body waves tend to dominate the reverberations of the centre-5 group, and surface waves those of the outer-6 group (Figs 8a and d, grey histograms).
For the outer-6 results we considered correlations between each three stations on either side of the centre-5 group, and all pairs crossing the central part.A separation does not change the results.When the position of the centre-5 group shifts by some 100 m in both fault-normal directions, the body wave dominance (central window, result 1) vanishes.It only remains persistent for the P 2 processed earthquake data, which supports again the notion that 1-bit clipped earthquake waveforms yield the strongest body wave reverberations.Body waves in the fault centre can be reconstructed from all processed data sets in the central window.It suggests that they attenuate less than the surface wave components that lead to more homogeneous phase distributions in the early window.The rate at which this mixing occurs can be a signature of the scattering length scales (Obermann et al. 2013).For the outer-6 group, body waves can still be recovered from earthquake waveforms, which implies, yet again, that earthquake signals are more efficient in exciting fault zone reverberations compared to the noise field.
The comparison of results from the centre-5 and outer-6 station groups reveals that important differences in the correlation coda wavefield are controlled by the relative position of the stations normal to the fault.In Fig. 9, we further investigate this by plotting the average fraction of body wave phases as a function of distance from the mapped fault trace and as a function of lapse time.Here we use correlations based on P 2 processed earthquake waveforms.These estimates are obtained from matrices (Fig. 9a) that contain the fraction of body wave phases per 2 s analysis window (the 0.1 s increment yields 21 estimates per window) per station pair.Each datum in Fig. 9(b) is then the average fraction along a row (or column) of this matrix associated with a station and hence position.The average is found by taking the mean over all values within a distance from each reference location.This distance is half the length of the line.We neglect again time windows too close to τ = 0 s to avoid near-field biases.Ideally, the reference lapse time would be the distance dependent main arrival time instead of the here used constant τ value.However, the obtained profiles are not very sensitive to distance or quality dependent weighted averaging.That is, the main arrival potentially captured in early lapse time windows between distant stations is not biasing our conclusions.
In Fig. 9(b), the indicated fraction is generally larger than 50 per cent, which suggests that the correlation coda field consists of more body wave polarization states compared to surface waves.However, this inference has to be moderated, because it is not clear if the propagation direction is always parallel to the fault-normal ZR component.The profiles show a relatively constant body wave fraction at negative distances away from the fault.This level is lower compared to average estimates around the fault zone, but does not indicate pure surface wave propagation (body wave fraction <0.5).We verified that our phase estimator returns the expected 90 degrees phase shift if a window is centred on a surface wave train.However, the values quickly deviate if the window is not centred on the main arrival.At early lapse times, correlations between stations located around the velocity gradients at small negative and positive distances exhibit the largest fraction of body wave phases (two dark-grey peaks in Fig. 9b), compatible with Figs 8(b) and (e).These two peaks merge with increasing lapse time into one peak around x = 0 km, where the peak value decreases from 75 per cent to around 70 per cent.Overall, Fig. 9(b) reveals a complex spatiotemporal evolution of the wavefield composition.The distributions become more homogeneous towards later lapse times.These results, too, are consistent with our interpretations based on Figs 7 and  8.The partitioning at early lapse times reflects the lateral velocity variations and hence fault structure.Mode conversion and scattering eventually homogenize the wavefield (Paul et al. 2005).At later lapse times this leads to partition coefficients that are governed by the local velocity structure.The stabilization of H/V ratios (where H and V denote energy estimates from horizontal and vertical component data, respectively) to a constant value is a characteristics of modal equipartitioning (Weaver 1982;Hennino et al. 2001).These ratios can be used to constrain properties of the propagation medium (Margerin et al. 2009;Sánchez-Sesma et al. 2011).Our conclusion that the fault structure affects the spectral ratio is supported by the fault normal H 2 /V 2 profiles that Hillers et al. (2014) obtained from the raw noise data of the Calico line array.These profiles exhibit the strongest variations in the regions of the velocity gradient, and are thus similar to the discussed early-window results (Fig. 9b).

R E V E R B E R AT I O N -B A S E D I M A G I N G
We compare synthetic and observed Green's functions to estimate the width of the trapping low-velocity zone (W, Fig. 1a), the wave speed within that zone (c LVZ ) and in the surrounding medium (c bulk ) and thus the velocity contrast.This comparison illustrates the possibility to constrain these fault zone properties using the reverberations.We aim to reproduce only a limited set of features of the complete yet complex reconstructed wavefield, that is, we focus on features in the arrival time structure of the ZR patterns that are associated with lateral surface wave propagation.We use the 2-D analytical solution of Ben-Zion & Aki (1990) and Ben-Zion (1998) to compute the scalar wavefield in a vertical low-velocity zone in response to an SH line source.This tool is particularly selected to efficiently replicate the main features of the arrival structure that are, again, associated with fault normal propagation.The polarizations states, that is, the difference between the SH excitation and the target ZR and ZZ fields, do not matter in this kinematic approach.
We neglect the effect of attenuation on the arrival time pattern, which is here limited to the amplitude decay rate with traveltime.This is in contrast to waveform fits of trapped waves interference patterns propagating along a waveguide, for which the Q values in this model are an important tuning parameter (Peng et al. 2003).The medium is uniform in the along-strike direction and with depth.We place the source at the surface and at x = 0 km.We sample the 1-2 Hz filtered synthetic wavefield in a 100 m interval (Fig. 10a).To construct similar images from the P 2 processed earthquake data we binned the ZZ and ZR correlation functions associated with a reference station at the origin (Figs 10b and c).To optimize the images, we used a 60 m sampling interval that accounts for the variable interstation distances (Fig. 2c).
The synthetic arrival time structure in Fig. 10(a) matches relevant features of the ZR pattern in Fig. 10(b) that we associate with lateral propagation.We highlight these features with the dotted lines in the upper half of Fig. 10(b), that is, the criss-cross or ' '-shaped pattern inside the low-velocity zone, and the branches of transmitted energy that propagate away from that zone at higher speeds.The ' ' epitomizes the dotted template.We manually the free model parameters (W, c LVZ , c bulk ) to reproduce the position of the reflection points at x ≈ ±0.4 km (Fig. 10b), to match the slope of the low-velocity ' ' pattern, the position of the second reflected crossing at τ ≈ 3 s and the slope of the transmitted wave fronts.The parameter combination W ≈ 800 m, c LVZ ≈ 800 m s −1 , c bulk ≈ 1600 m s −1 leads to the synthetic wavefield shown in Fig. 10(a), which duplicates fairly good the ZR main features.This is a qualitative statement, but we find that changing any of these values by more than 10-20 per cent produces patterns that are clearly different from the dotted main features in Fig. 10(b).Since we focus on lateral surface wave propagation in the ZZ-ZR plane c corresponds to Rayleigh wave speed.The obtained values are compatible with the estimates based on a traveltime inversion (Hillers et al. 2014).To be exact, the empirical time-domain Green's function is proportional to the time derivative of the cross-correlation function.We compared the correlation patterns (Figs 10b and c) with those of the associated time derivatives and found that the differences including the small phase shift are negligible and do not affect the above W, c LVZ and c bulk estimates.
The synthetic wavefield differs in many ways from the observations.The reversed polarity of the computed waveforms across x = 0 km is not found in the ZR and ZZ observations.This is governed by the different source mechanisms: The synthetics correspond to a SH dislocation, whereas the empirical Green's function estimates are the response to a vertical point force.We iterate that this does not impair the here adopted exploratory approach.A more informative mismatch is that the synthetic wavefield contains less features, for example, it does not reproduce the additional transmitted wave front in the ZR data that emerges between the highlighted branches that depart at the reflection points.In general, there is a larger discrepancy between the synthetics and the ZZ data (Fig. 10c) compared to the ZR data.The ZZ arrival times do not exhibit the ' ' structure, the strongest features have a ' ' geometry instead.The upper parts at x > 0 km are again highlighted by dotted lines in Fig. 10(c).The associated slopes become steeper with increasing lapse time, that is, later arriving wave fronts have a higher apparent velocity.This effect is likely governed by velocity changes in deeper parts of the structure, as waves can interact with the bottom of the low-velocity zone (Cochran et al. 2009;Yang et al. 2011).
We conclude that the employed simplistic scalar model allows us to constrain key average structural properties of the trapping lowvelocity zone by focusing on a limited set of distinct features in the ZR arrival time structure that are associated with fault-normal propagation.The obtained effective width W ≈ 0.8-1 km and the ∼50 per cent velocity contrast indicated in Fig. 9(c) (black) are compatible with the estimates based on different methods discussed below.The approach illustrates the possibility to constrain the lateral velocity variations in a fault zone environment from properties of the reverberation wavefield.However, the model cannot reproduce the full wavefield complexity.The observed richer behaviour is likely governed by fault-normal and depth dependent velocity variations of the possibly wedge-shaped low-velocity zone.A better reconstruction of the wavefield constituents for improved fault zone imaging is expected from the application of a complete wavefield solver.This can potentially be extended into a reverberation full waveform inversion approach considering the full network response consisting of the 3 × 3 correlation tensor between all station pairs.

D I S C U S S I O N
Several conclusions on fault zone-wavefield interactions can be drawn from our analysis of cross-correlation wavefields obtained from a dense line array in a strike-slip faulting environment.Multiple surface and body wave reflections emerge in the early coda part of the Green's function estimates that are constructed from earthquake waveforms and the ambient noise (Figs 4 and 5).The composition of the reverberation wavefield depends strongly on the location of station pairs relative to the low-velocity waveguide and on the traveltime (Figs 8 and 9).Mode partitioning depends also on the wavefield from which the reverberations are constructed (earthquake or noise) and on the employed signal processing.
The robustness of these phases is highlighted by the dense spatial sampling.The interstation distance is sufficiently small compared to the wavelengths at 1-2 Hz, which are in turn smaller than the Calico damage zone width of 1-2 km.Different spectral and temporal normalization choices can tune the partitioning of the reconstructed wavefield components in the early part of the correlation coda.The strongest surface wave signatures are found if the original data are not whitened.The application of whitening and especially 1-bit clipping enhances strong multiple reflections in the correlation waveforms.Processing earthquake data tends to yield stronger signatures of bulk waves compared to surface wave modes.This is compatible with the understanding that the earthquake source process excites more body waves compared to the noise wavefield, which is usually dominated by surface waves.
Fault-normal profiles of the correlation coda composition indicate systematic spatial variations of local wavefield properties (Fig. 9).The overall higher level of longitudinal body wave polarization in the region around the largest velocity reduction is synonymous with the reduction of surface wave propagation modes that is governed by scattering, mode conversion and absorption within or at the boundaries of the compliant damage zone.Variable scattering and attenuation of SV-and P-waves also control the lapse time dependent phase difference between ZZ and ZR motion.The strongest variability of the wavefield anatomy with propagation time is again observed at stations located in the central part of the damage zone.The zoning of this trend implies that this behaviour is governed by evolving properties of the reverberating field.
The emergence of multiples supports the interpretation that the low-velocity zone acts as a reverberant cavity and waveguide (Hillers et al. 2014).The here discussed configuration constitutes an open system, which does not, however, undermine the usefulness of the cavity analogy.The actual source is outside the trapping structure, and in-fault modes are excited trough scattering and mode conversion that do not contain signatures of the earthquake source mechanism.Reverberations excited in the waveguide by the impinging earthquake wavefield are transient.The good agreement between stacked correlation waveforms from earthquake coda and noise data demonstrates the reconstruction of phases that are indeed part of the Green's function.We can conclude that the correlation functions constructed from properly processed earthquake signals-and surely from noise-represent the Earth's response, where correlation coda phases are associated with multiple reflections inside the low-velocity zone.
The waveguiding potential of the Calico fault was first suggested by the space geodesy-based observation of a 1-2 km wide compliant zone (Fialko et al. 2002) and then confirmed with the observation of ballistic trapped waves (Cochran et al. 2009).The estimated Pwave velocity structure (Fig. 9c, solid grey line) shows a ∼1 km wide zone of 30-40 per cent velocity reduction that narrows with depth.A traveltime analysis of high-frequency body waves implied a ∼1.3 km wide dipping low-velocity zone of 40-50 per cent wave speed reduction that extends to 3 km depth (Yang et al. 2011).The surface wave velocity profile obtained from an ambient noise analysis (Hillers et al. 2014) shows a zone of reduced velocities (Fig. 9c, dashed grey line) that is broader compared to the estimates based on body waves.The average lateral velocity structure obtained here from properties of the reverberation wavefield (Fig. 9c, black line) shows again a more narrow zone (∼1 km) of reduced wave speeds.
Together, these results indicate a zone of reduced velocities that likely narrows with depth, and that has a finite depth extend H (Fig. 1a).This is commonly observed in strike-slip environments (Sylvester 1988).The wedge-shape of these 'flower structures' is mainly controlled by the depth dependent increase in the confining pressure (Ma 2008;Finzi et al. 2009).The low-rigidity material in this zone accommodates distributed deformation during an earthquake rupture and may thus partially explain the frequently observed shallow slip deficit (Kaneko & Fialko 2011).The generally depth dependent velocity increase possibly controls the wave fronts in the ZZ correlation patterns at late lapse times.These arrivals can be associated with vertical propagation, and the moveout patterns can be partially controlled by subhorizontal impedance contrasts (Draganov et al. 2007) in addition to the lateral variations discussed here.
Similar to signatures of trapped noise (Hillers et al. 2014), reverberations are expected to emerge along fault segments that do not sustain ballistic trapped waves due to the irregularity of the waveguide.Our results demonstrate the possibility to constrain properties of the damage low-velocity zone from reconstructed multiples.At the same time the results expose the limits of the employed scalar model for the reproduction of the complex wavefield.Meeting or superseding the resolution of previous results requires a synthetic vectorial wavefield together with the application of inversion techniques to phase arrival times or to full reverberation waveforms.Since the reverberation wavefield constitutes a typical signature of the fault, it extends the range of fault zone waves that can be used for the analysis of fault structure.The reconstruction of the multiples from stacks of earthquake waveform correlations shows that earthquake signals can contribute to fault zone imaging even if they do not contain clean signatures of head, trapped or reflected waves.The reconstruction of the multiples from the ambient noise field makes it possible to characterize faults in regions of low background seismicity.This is important for improved seismic hazard assessment and strong ground motion estimates along active faults that are not frequently illuminated by earthquake waves.
The robust noise-based reconstruction of multiple reflections within the damage zone is also important for probing the state of the wear zone.Noise-based monitoring can resolve variations of the elastic and structural properties of a heterogeneous medium (e.g.Campillo et al. 2011;Snieder & Larose 2013, and references therein).Our results imply that high-quality velocity change measurements can be obtained from the analysis of the fault zone reverberations.Importantly, this wavefield is only sensitive to changes in the elastic properties of the damage zone material.The body wave component in the reverberation wavefields suggests that these phases can resolve material changes below the surface wave penetration depth, although more research is necessary to understand the vertical propagation history of the observed wave fronts.
Passive monitoring of crustal material properties in response to a variety of external and internal loading mechanisms is well established, and the here presented results constitute important groundwork for the application of these methods to fault zone rocks.Considering the desert environment and the 6-month duration of the Calico fault zone experiment, the data set is suitable for noise-based investigations of the damage zone response to diurnal temperature changes (Richter et al. 2014) and tidal deformation (Hillers et al. 2015).The resolution of damage zone responses to strain variations characterized by variable amplitude and duration will further benefit from modern seismic array designs (Hand 2014;Ben-Zion et al. 2015).Continuous dense deployments across fault zones can provide important observations for the study of potential triggering mechanisms (Gomberg & Johnson 2005;Johnson et al. 2008).As shown here, the resolution of such phenomena appears most promising in fault zone environments that sustain reverberations.

C O N C L U S I O N S
Using data from a densely instrumented linear array across the Calico fault zone, southern California, we studied signatures of multiple reflections in cross-correlation functions reconstructed from scattered wavefields.These reverberations emerge within the cavity/waveguide-like structure of the low-velocity channel of the fault damage zone and consist of a mix of body wave and surface wave propagation modes.The reverberation arrivals in the crosscorrelation coda can be obtained from earthquake waveforms and from noise records.The consistency of the associated correlation coda wavefields is a strong indicator for the physical reconstruction of the Green's function.Modal partitioning of the emergent reverberations can be tuned using different normalization choices that enhance or mute different components of the original wavefields.We constrained key properties of the damage low-velocity zone (width, velocity contrast) from the arrival structure of the reverberation wavefield that are compatible with estimates obtained with more traditional techniques.Multiples extend the range of fault zone waves that can be used for imaging and the study of time-dependent fault zone properties.The emergence from noise facilitates the characterization of faults in regions of low seismicity.Monitoring the properties of signals reconstructed from waves that are multiply reflected inside the low velocity zone can be used to probe the damage zone state.

A C K N O W L E D G E M E N T S
Discussions with P. Roux, J. Chaput, L. Stehly, P. Boué, F. Brenguier and A. Colombi helped to organize our results.We thank Y. Ben-Zion for discussions and for providing the numerical tool for the computation of the synthetic waveforms.We thank E. Cochran for providing material associated with the Calico field experiment, and X. Briand for computational assistance.Comments by the editor I. Grevemeyer, reviewer S. de Ridder and two anonymous reviewers helped to improve the manuscript.GH acknowledges support from a Heisenberg fellowship of the German Research Foundation (HI 1714/1-2).This work was supported by the European Research Council (Advanced grant Whisper L27507).Figures were made using GMT (Wessel & Smith 1998).The facilities of IRIS Data Services, and specifically the IRIS Data Management Center, were used for access to waveforms used in this study.IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience and EarthScope (SAGE) Proposal of the National Science Foundation under Cooperative Agreement EAR-1261681.

Figure 1 .
Figure 1.(a) Conceptual model of a low-velocity damage zone associated with strike-slip faulting.The width W, depth H, dip angle and the impedance contrast control signatures of ballistic and scattered wavefields in a fault zone environment.W commonly decreases with depth.For the Calico fault zone, W is on the order of 1-2 km; x indicates the fault-normal direction.(b) Map view of the situation in panel (a).Solid arrows indicate the lateral propagation pattern of the reverberating wavefield.The configuration represents a ray pattern that differs from the associated traveltime pattern (if the along-strike direction is replaced by traveltime as in Fig.10a).

Figure 2 .
Figure 2. (a) Map of the study area.The circles indicate the locations of the 25 earthquakes used in this study.Grey lines show geologically mapped faults.The cross indicates the location of the fault zone array.(b) The solid and dashed lines show P-wave and surface wave phase velocity profiles from Cochran et al. (2009) and Hillers et al. (2014), respectively.(c) Station geometry of the line array.'Centre-5' and 'outer-6' label groups of stations that are used in the phase analysis (Section 4).The origin is indicated by the cross in the map.The arrow indicates North.

Figure 3 .
Figure 3. (a) Bandpass filtered seismograms of a regional M3.5 earthquake.(b) The envelopes of the filtered seismograms indicate an increased noise level at low frequencies.Coda energy is dissipated after 3-5 min in the 1-2 Hz target frequency range.

Figure 5 .
Figure 5. ZZ cross-correlation wavefields constructed with the P 2 processing scheme in the 1-2 Hz range.Amplitudes are scaled by the peak value in each panel.(a) The same data and configuration as in Fig. 4(a), but now the reverberations have larger amplitudes.(b) The same results as in panel (a) for the reference station S ref located closer to the mapped fault centre.Amplitudes of later reverberation arrivals are increased relative to the surface wave.(c) The in-fault correlation wavefield constructed from 25 5 min long earthquake waveforms is very similar to the noise result in panel (b), but the reverberation amplitudes are larger.(d) The same situation as in panel (c) for a different earthquake coda window.

Figure 6 .
Figure 6.Contributions of individual earthquake and noise ZZ correlations to the resulting stacks for a station pair located at x = 0 m and x = 125 m.(a) Correlations constructed from 2 min long earthquake coda windows (t 0 + 1 to t 0 + 3 min), and the resulting stack (top, black waveform).For comparison, we plot the stack associated with a later coda window (grey waveform).(b) The correlations and the stack correspond to the original 5 min window.The grey stack is the black stack from panel (a).(c) 19 daily noise correlations and the corresponding stack (top, black waveform).The grey stack is the black stack from panel (b).Distance d in kilometre, azimuth az in degree and magnitude M of the 25 earthquakes used in panels (a) and (b) are indicated to the left.In panels (a) and (b), data are ordered from bottom to top with increasing similarity to the average stack based on the data in panel (b).All data are P 2 processed.

Figure 9 .
Figure9.(a) Phase partition pattern between all stations for the 3-5 s coda window obtained from P 2 processed earthquake waveforms (see also Figs8b and e).Symbol size and grey scale indicate the fraction of body wave phase estimates.Large, dark (small, bright) symbols represent a large (small) fraction.(b) Average fraction of body wave polarization as a function of distance from the fault obtained from overlapping 2 s long analysis windows as in panel (a).The timing of the window-onset is indicated.(c) Velocity profiles from Fig.2(b) (grey lines) and the estimate obtained from our reverberations-based imaging approach in Section 5 (black line).

Figure 10 .
Figure 10.Synthetic and observed Green's functions.(a) Scalar wavefield computation based on the solutions of Ben-Zion & Aki (1990).The imposed low-velocity zone extends from x = −0.4km to x = 0.4 km.(b) ZR and (c) ZZ correlations from P 2 processed earthquake waveforms.Amplitudes are scaled by the peak value in each panel, and the grey scales are adjusted in each figure.Equidistant binning in panels (b) and (c) averages over the variable interstation distances (Fig. 2c).Dotted lines in panels (b) and (c) illustrate relevant, coherent wave fronts.

Table 1 .
Summary of the processing sequences.The initial highpass filter for P 0 and P 2 is put in parentheses since it is not pivotal.