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

Gregor Hillers, Michel Campillo

► To cite this version:

Gregor Hillers, Michel Campillo. Fault zone reverberations from cross-correlations of earthquake waveforms and seismic noise. Geophysical Journal International, 2016, 204, pp.1503-1517. 10.1093/gji/ggv515. insu-03596380

HAL Id: insu-03596380

https://insu.hal.science/insu-03596380

Submitted on 3 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Gregor Hillers and Michel Campillo

Institut des Sciences de la Terre, Université Joseph Fourier, CNRS, F-38400 Grenoble, France. E-mail: gregor.hillers@uff-grenoble.fr

Accepted 2015 November 26. Received 2015 November 24; in original form 2015 June 30

## SUMMARY

Seismic wavefields interact with low-velocity fault damage zones. Waveforms of ballistic fault zone head waves, trapped waves, reflected 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 significant, 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 wavefields. Modal partitioning of the reverberation wavefield 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.

**Key words:** Interferometry; Fault zone rheology; Guided waves; Coda waves; Wave propagation.

## 1 INTRODUCTION

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. 1990, 1994; Li & Leary 1990; Ben-Zion 1998). Third, body wave reflections at fault zone boundaries form typical P- and S-wave move-out 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
<figure><img src="image_3.png" /><figcaption>**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).</figcaption></figure>

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 above-mentioned 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 low-velocity 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.

# 2 METHOD

## 2.1 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 low-velocity 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.
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

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.

## 2.2 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,

<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>

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 $\tau \approx 10\bar{T}$, where $\tau$ denotes the lapse time in the correlation function and $\bar{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_{i0}$ and
<table><caption>**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.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Before correlation</th><th rowspan="2">After correlation<br>Post-filtering</th></tr><tr><th>Pre-filtering</th><th>f-domain normalization</th><th>t-domain normalization</th></tr></thead><tbody><tr><td>$P_0$</td><td>(Highpass filter >0.01 Hz)</td><td>Whiten 1–2 Hz</td><td>Threshold clipping</td><td>—</td></tr><tr><td>$P_1$</td><td>Bandpass filter 0.03–5 Hz</td><td>—</td><td>1-bit normalization</td><td>Bandpass filter 1–2 Hz</td></tr><tr><td>$P_2$</td><td>(Highpass filter >0.01 Hz)</td><td>Whiten 1–2 Hz</td><td>1-bit normalization</td><td>—</td></tr></tbody></table>

$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 in-fault reverberations, and contributions to the converging correlation stacks are systematically investigated.

# 3 RESULTS

## 3.1 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 $\tau$-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).

## 3.2 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$ noise-based 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 signal-to-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><img src="image_6.png" /><figcaption>**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 $\tau$, 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.</figcaption></figure>

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 low-frequency limit of the initial bandpass filter (Table 1). We found that increasing $f_{10}$ 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_{10}$ 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).

### 3.3 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_{\text{ref}}$ at $x = -0.5$ km and stations across the zone of reduced velocities (Fig. 5a) exhibit a nearly amplitude-symmetric surface
<figure><img src="image_7.png" /><figcaption>**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_{\text{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.</figcaption></figure>

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_{\text{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 incipient 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_{\text{ref}}$ and stations around $x = 0$ km. These arrivals can be clearly identified at least up to $|\tau| = 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_{\text{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_{\text{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 $|\tau| = 8$ s. In conclusion, the strongest signature of the reverberations can be extracted from both wavefields with the $P_2$ processing (target-band whitening, 1-bit clipping).
Correlations from earthquake waveforms

<figure><img src="image_8.png" /></figure>

Noise correlations

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**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*₀ + 1 to *t*₀ + 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*₂ processed.</figcaption></figure>

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.

## 3.4 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><img src="image_11.png" /><figcaption>**Figure 7.** Lapse time or τ-averaged correlation coda from *P*₁ and *P*₂ 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.</figcaption></figure>

$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 ($\tau = -20$ s to $\tau = 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.

## 4 ANATOMY OF THE REVERBERATIONS

We estimate the phase $\phi$ 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 $\phi(\tau)$ 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 $\phi(\tau)$ variations for
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Histograms of phase estimates between ZZ and ZR correlation coda waves obtained with $P_1$ and $P_2$ processed earthquake waveforms and noise. Panels in the top (bottom) row correspond to the centre-5 (outer-6) group of stations (Fig. 2). Phases around 90° (45° < $\phi_{surface}$ < 135°) are associated Rayleigh surface wave polarization, and phases around 0° / 180° ($\phi_{body}$ < 45° and $\phi_{body}$ > 135°) indicate longitudinally polarized body waves. The centre-5 (outer-6) histograms for each 2 s window are compiled from 21 × 10 (15) = 210 (315) individual measurements.</figcaption></figure>

these groups and discuss the systematic lapse time dependence of $\phi$ distributions obtained in three 2 s long windows between $\tau = 1.5$ s and $\tau = 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; $\phi$ 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° < $\phi$ < 135°) and body waves (0° and 180° or $\phi$ < 45° and $\phi$ > 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 wavefields 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
<figure><img src="image_13.png" /><figcaption>**Figure 9.** (a) Phase partition pattern between all stations for the 3–5 s coda window obtained from $P_2$ processed earthquake waveforms (see also Figs 8b 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).</figcaption></figure>

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 $\tau = 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 $\tau$ 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).

# 5 REVERBERATION-BASED IMAGING

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
<figure><img src="image_14.png" /><figcaption>**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.4$ km 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.</figcaption></figure>

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 adjusted the free model parameters (*W*, $c_{LVZ}$, $c_{bulk}$) to reproduce the position of the reflection points at $x \approx \pm 0.4$ km (Fig. 10b), to match the slope of the low-velocity ‘⟨⟩⟨⟩’ pattern, the position of the second reflected crossing at $\tau \approx 3$ s and the slope of the transmitted wave fronts. The parameter combination $W \approx 800$ m, $c_{LVZ} \approx 800$ m s⁻¹, $c_{bulk} \approx 1600$ m s⁻¹ 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 low-velocity 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 \approx 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.

## 6 DISCUSSION

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 *P*-wave 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.

## 7 CONCLUSIONS

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 cross-correlation 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.

## ACKNOWLEDGEMENTS

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.

## REFERENCES

Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., 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**, 1239–1260.

Ben-Zion, Y., 1998. Properties of seismic fault zone waves and their utility for imaging low-velocity structure, *J. geophys. Res.*, **103**(B6), 12 567–12 585.

Ben-Zion, Y., 2008. Collective behavior of earthquakes and faults: continuum-discrete transitions, progressive evolutionary changes and different dynamic regimes, *Rev. Geophys.*, **46**(4), doi:10.1029/2008RG000260.

Ben-Zion, Y. & Aki, K., 1990. Seismic radiation from an SH line source in a laterally heterogeneous planar fault zone, *Bull. seism. Soc. Am.*, **80**(4), 971–994.

Ben-Zion, Y. & Sammis, C.G., 2003. Characterization of fault zones, *Pure appl. Geophys.*, **160**, 677–715.

Ben-Zion, Y. *et al.*, 2015. Basic data features and results from a spatially-dense seismic array on the San Jacinto fault zone, *Geophys. J. Int.*, **202**, 370–380.

Boué, P., Poli, P., Campillo, M., Pedersen, H., Briand, X. & Roux, P., 2013. Tele seismic correlations of ambient seismic noise for deep global imaging of the Earth, *Geophys. J. Int.*, **194**, 844–848.

Brenguier, F., Campillo, M., Hadzioioannou, C., Shapiro, N.M., Nadeau, R.M. & Larose, E., 2008. Postseismic relaxation along the San Andreas Fault at Parkfield from continuous seismological observations, *Science*, **321**(5895), 1478–1481.

Campillo, M., 2006. Phase and correlation of ‘random’ seismic fields and the reconstruction of the Green function, *Pure appl. Geophys.*, **163**, 475–502.

Campillo, M. & Paul, A., 2003. Long-range correlations in the diffuse seismic coda, *Science*, **299**(5606), 547–549.

Campillo, M., Sato, H., Shapiro, N.M. & van der Hilst, R.D., 2011. New developments on imaging and monitoring with seismic noise, *C. R. Geosci.*, **343**, 487–495.

Catheline, S., Benech, N., Brum, J. & Negreira, C., 2008. Time reversal of elastic waves in soft solids, *Phys. Rev. Lett.*, **100**(6), 064301, doi:10.1103/PhysRevLett.100.064301.

Chaput, J., Campillo, M., Aster, R.C., Roux, P., Kyle, P.R., Knox, H. & Czoski, P., 2015. Multiple scattering from icequakes at Erebus volcano, Antarctica; implications for imaging at glaciated volcanoes, *J. geophys. Res.*, **120**, 1129–1141.

Cochran, E.S., Li, Y.-G., Shearer, P.M., Barbot, S., Fialko, Y. & Vidale, J.E., 2009. Seismic and geodetic evidence for extensive, long-lived fault damage zones, *Geology*, **37**(4), 315–318.

Cupillard, P., Stehly, L. & Romanowicz, B., 2011. The one-bit noise correlation: a theory based on the concepts of coherent and incoherent noise, *Geophys. J. Int.*, **184**, 1397–1414.

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**, 2102–2118.

Draganov, D., Wapenaar, K., Mulder, W., Singer, J. & Verdel, A., 2007. Retrieval of reflections from seismic background-noise measurements, *Geophys. Res. Lett.*, **34**, doi:10.1029/2006GL028735.

Fialko, Y., Sandwell, D., Agnew, D., Simons, M., Shearer, P. & Minster, B., 2002. Deformation on Nearby Faults Induced by the 1999 Hector Mine Earthquake, *Science*, **297**, 1858–1862.

Finzi, Y., Hearn, E.H., Ben-Zion, Y. & Lyakhovsky, V., 2009. Structural properties and deformation patterns of evolving strike-slip
faults: numerical simulations incorporating damage rheology, *Pure appl. Geophys.*, **166**, 1537–1573.

Fukao, Y., Hori, S. & Ukawa, M., 1983. A seismological constraint on the depth of basalt-eclogite transition in a subducting oceanic crust, *Nature*, **303**, 413–415.

Gomberg, J. & Johnson, P., 2005. Dynamic triggering of earthquakes, *Nature*, **437**, 830, doi:10.1038/437830a.

Gouédard, P., Roux, P., Campillo, M. & Verdel, A., 2008. Convergence of the two-point correlation function toward the Green's function in the context of a seismic-prospecting data set, *Geophysics*, **73**(6), 47–53.

Groos, J.C., Bussat, S. & Ritter, J.R.R., 2012. Performance of different processing schemes in seismic noise cross-correlations, *Geophys. J. Int.*, **188**, 498–512.

Hand, E., 2014. A boom in boomless seismology, *Science*, **345**, 720–721.

Hennino, R., Trégourès, N., Shapiro, N.M., Margerin, L., Campillo, M., van Tiggelen, B.A. & Weaver, R.L., 2001. Observation of equipartitioning of seismic waves, *Phys. Rev. Lett.*, **86**(15), 3447–3450.

Hillers, G., Campillo, M., Ben-Zion, Y. & Roux, P., 2014. Seismic fault zone trapped noise, *J. geophys. Res.*, **119**(7), 5786–5799.

Hillers, G., Retailleau, L., Campillo, M., Inbal, A., Ampuero, J.-P. & Nishimura, T., 2015. In situ observations of velocity changes in response to tidal deformation from analysis of the high-frequency ambient wave-field, *J. geophys. Res.*, **120**, 210–225.

Hobiger, M., Wegler, U., Shiomi, K. & Nakahara, H., 2012. Coseismic and postseismic elastic wave variations caused by the 2008 Iwate-Miyagi Nairiku earthquake, Japan, *J. geophys. Res.*, **117**, doi:10.1029/2012JB009402.

Johnson, P.A., Savage, H., Knuth, M., Gomberg, J. & Marone, C., 2008. Effects of acoustic waves on stick-slip in granular media and implications for earthquakes, *Nature*, **451**, 57–60.

Kaneko, Y. & Fialko, Y., 2011. Shallow slip deficit due to large strike-slip earthquakes in dynamic rupture simulations with elasto-plastic off-fault response, *Geophys. J. Int.*, **186**, 1389–1403.

Landès, M., Hubans, F., Shapiro, N.M., Paul, A. & Campillo, M., 2010. Origin of deep ocean microseisms by using teleseismic body waves, *J. geophys. Res.*, **115**, doi:10.1029/2009JB006918.

Larose, E., Roux, P. & Campillo, M., 2007. Reconstruction of Rayleigh-Lamb dispersion spectrum based on noise obtained from an air-jet forcing, *J. acoust. Soc. Am.*, **122**(6), 3437–3444.

Lewis, M.A. & Ben-Zion, Y., 2010. Diversity of fault zone damage and trapping structures in the Parkfield section of the San Andreas Fault from comprehensive analysis of near fault seismograms, *Geophys. J. Int.*, **183**, 1579–1595.

Lewis, M.A., Peng, Z., Ben-Zion, Y. & Vernon, F.L., 2005. Shallow seismic trapping structure in the San Jacinto fault zone near Anza, California, *Geophys. J. Int.*, **162**, 867–881.

Li, Y.-G. & Leary, P., 1990. Fault zone trapped seismic waves, *Bull. seism. Soc. Am.*, **80**(5), 1245–1271.

Li, H., Zhu, L. & Yang, H., 2007. High-resolution structures of the Landers fault zone inferred from aftershock waveform data, *Geophys. J. Int.*, **171**, 1295–1307.

Li, Y.-G., Leary, P., Aki, K. & Malin, P., 1990. Seismic trapped modes in the Oroville and San Andreas fault zones, *Science*, **249**(4970), 763–766.

Li, Y.-G., Aki, K., Adams, D., Hasemi, A. & Lee, W.H.K., 1994. Seismic guided waves trapped in the fault zone of the Landers, California, earthquake of 1992, *J. geophys. Res.*, **99**(B6), 11 705–11 722.

Lin, F.-C., Li, D., Clayton, R.W. & Hollis, D., 2013a. 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.

Lin, F.-C., Tsai, V.C., Schmandt, B., Duputel, Z. & Zhan, Z., 2013b. Extracting seismic core phases with array interferometry, *Geophys. Res. Lett.*, **40**, 1049–1053.

Lobkis, O.I. & Weaver, R.L., 2001. On the emergence of the Green's function in the correlations of a diffuse field, *J. acoust. Soc. Am.*, **110**(6), 3011–3017.

Ma, S., 2008. A physical model for widespread near-surface and fault zone damage induced by earthquakes, *Geochem. Geophys. Geosyst.*, **9**(11), doi:10.1029/2008GC002231.

Margerin, L., Campillo, M., van Tiggelen, B.A. & Hennino, R., 2009. Energy partition of seismic coda waves in layered media: theory and application to Pinyon Flats Observatory, *Geophys. J. Int.*, **177**, 571–585.

Nakata, N., Chang, J.P., Lawrence, J.F. & Boué, P., 2015. Body wave extraction and tomography at Long Beach, California, with ambient-noise interferometry, *J. geophys. Res.*, **120**, 1159–1173.

Obermann, A., Planès, T., Larose, E., Sens-Schönfelder, C. & Campillo, M., 2013. Depth sensitivity of seismic coda waves to velocity perturbations in an elastic heterogeneous medium, *Geophys. J. Int.*, **194**(1), 372–382.

Paul, A., Campillo, M., Margerin, L., Larose, E. & Derode, A., 2005. Empirical synthesis of time-asymmetrical Green function from the correlation of coda waves, *J. geophys. Res.*, **110**, doi:10.1029/2004JB003521.

Peng, Z., Ben-Zion, Y., Michael, A.J. & Zhu, L., 2003. Quantitative analysis of seismic fault zone waves in the rupture zone of the 1992 Landers, California, earthquake: evidence for a shallow trapping structure, *Geophys. J. Int.*, **155**, 1021–1041.

Peterson, J., 1993. Observations and modeling of seismic background noise, Open file rep., U.S. Geol. Surv.

Poli, P., Campillo, M., Pedersen, H.A. & the LAPNET Working Group, 2012a. Body-wave imaging of Earth's mantle discontinuities from ambient seismic noise, *Science*, **338**, 1063–1065.

Poli, P., Pedersen, H.A., Campillo, M. & the POLENET/LAPNET Working Group, 2012b. Emergence of body waves from cross-correlation of short period seismic noise, *Geophys. J. Int.*, **188**, 549–558.

Prieto, G.A., Denolle, M., Lawrence, J.F. & Beroza, G.C., 2011. On amplitude information carried by the ambient seismic field, *C. R. Geosci.*, **343**, 600–614.

Richter, T., Sens-Schönfelder, C., Kind, R. & Asch, G., 2014. Comprehensive observation and modeling of earthquake and temperature-related seismic velocity changes in northern Chile with passive image interferometry, *J. geophys. Res.*, **119**, 4747–4765.

Roux, P., 2009. Passive seismic imaging with direct ambient noise: application to surface waves and the San Andreas Fault in Parkfield, CA, *Geophys. J. Int.*, **179**, 367–373.

Roux, P. & Ben-Zion, Y., 2014. Monitoring fault zone environments with correlations of earthquake waveforms, *Geophys. J. Int.*, **196**, 1073–1081.

Sabra, K.G., Gerstoft, P., Roux, P., Kuperman, W.A. & Fehler, M.C., 2005. Surface wave tomography from microseisms in Southern California, *Geophys. Res. Lett.*, **32**(14), doi:10.1029/2005GL023155.

Sánchez-Sesma, F.J. et al., 2011. A theory for microtremor H/V spectral ratio: Application for a layered medium, *Geophys. J. Int.*, **186**, 221–225.

SCEDC, 2013. Southern California Earthquake Center. Caltech.Dataset. doi:10.7909/C3WD3xH1.

Seats, K.J., Lawrence, J.F. & Prieto, G.A., 2012. Improved ambient noise correlation functions using Welch's method, *Geophys. J. Int.*, **188**, 513–523.

Shapiro, N.M. & Campillo, M., 2004. Emergence of broadband Rayleigh waves from correlations of the seismic ambient noise, *Geophys. Res. Lett.*, **31**, doi:10.1029/2004GL019491.

Shapiro, N.M., Campillo, M., Stehly, L. & Ritzwoller, M.H., 2005. High-resolution surface-wave tomography from ambient seismic noise, *Science*, **307**, 1615–1618.

Snieder, R. & Larose, E., 2013. Extracting Earth's elastic wave response from noise measurements, *Annu. Rev. Earth Planet. Sci.*, **41**, 183–206.

Stehly, L. & Hanasoge, S., 2012. Processing noise records with non-linear operations to select the distribution of noise sources contributing to the correlations, *EOS, Trans. Am. geophys. Un.*, AGU Fall Meet. Suppl., Abstract S51G-08.

Sylvester, A.G., 1988. Strike-slip faults, *Bull. geol. Soc. Am.*, **100**, 1666–1703.

Tanimoto, T., Ishimaru, S. & Alvizuri, C., 2006. Seasonality of particle motion of microseisms, *Geophys. J. Int.*, **166**, 253–266.
Verschuur, D.J., Berkhout, A.J. & Wapenaar, C.P.A., 1992. Adaptive surface-
related multiple elimination, *Geophysics*, **57**(9), 1166–1177.

Weaver, R.L., 1982. On diffuse waves in solid media, *J. acoust. Soc. Am.*, **71**(6), 1608–1609.

Wegler, U. & Sens-Schönfelder, C., 2007. Fault zone monitoring with passive image interferometry, *Geophys. J. Int.*, **168**, 1029–1033.

Wessel, P. & Smith, W.H.F., 1998. New, improved version of generic mapping tools released, *EOS, Trans. Am. geophys. Un.*, **79**, doi:10.1029/98EO00426.

Yang, H., Zhu, L. & Cochran, E.S., 2011. Seismic structures of the Calico fault zone inferred from local earthquake travel time modelling, *Geophys. J. Int.*, **186**, 760–770.

Zhan, Z., Ni, S., Helmberger, D.V. & Clayton, R.W., 2010. Retrieval of Moho-reflected shear wave arrivals from ambient seismic noise, *Geophys. J. Int.*, **182**(1), 408–420.

Zhang, J. & Gerstoft, P., 2014. Local-scale cross-correlation of seismic noise from the Calico fault experiment, *Earthq. Sci.*, doi:10.1007/s11589-014-0074-z.

Zhang, J., Gerstoft, P. & Shearer, P.M., 2009. High-frequency P-wave seismic noise driven by ocean winds, *Geophys. Res. Lett.*, **36**(9), doi:10.1029/2009GL037761.

Zigone, D., Ben-Zion, Y., Campillo, M. & Roux, P., 2015. Seismic tomography of the Southern California plate boundary region from noise-based Rayleigh and Love Waves, *Pure appl. Geophys.*, **172**, 1007–1032.