- Split View
-
Views
-
Cite
Cite
Philippe Roux, Yehuda Ben-Zion, Monitoring fault zone environments with correlations of earthquake waveforms, Geophysical Journal International, Volume 196, Issue 2, February, 2014, Pages 1073–1081, https://doi.org/10.1093/gji/ggt441
- Share Icon Share
Abstract
We develop a new technique for monitoring temporal changes in fault zone environments based on cross-correlation of earthquake waveforms recorded by pairs of stations. The method is applied to waveforms of ∼10 000 earthquakes observed during 100 d around the 1999 M 7.1 Duzce mainshock by a station located in the core damage zone of the North Anatolian Fault and a nearby station. To overcome clock problems, the correlation functions are realigned on a dominant peak. Consequently, the analysis focuses on measurements of coherency rather than traveltimes, and is associated with correlation coefficient of groups of events with a reference wavelet. Examination of coherency in different frequency bands reveals clear changes at a narrow band centred around 0.8 Hz. The results show a rapid drop of ∼1–2 per cent of the coherency at the time of the Duzce event followed by gradual recovery with several prominent oscillations over 4 d. The observed changes likely reflect evolution of permeability and fluid motion in the core damage zone of the North Anatolian Fault. Compared to noise correlation processing, our analysis of earthquake waveform correlation (i) benefits from high level of coherence with short duration recorded signals, (ii) has considerably finer temporal sampling of fault dynamics after mainshocks than is possible with noise correlation, (iii) uses the coherence level to track property variations, which may be more robust than traveltime fluctuations in the coda of noise correlations. Studies utilizing both earthquake and noise waveforms at multiple pairs of stations across fault damage zones can improve significantly the understanding of fault zone processes.
1 INTRODUCTION
Natural fault zones include narrow layers (tens to hundreds metres wide) of highly damaged rocks and bimaterial interfaces that contain important information on past and future earthquake activity (e.g. Ben-Zion & Sammis 2003, and references therein). Large earthquakes produce abrupt reductions of seismic velocities (brittle rock damage) near the rupture zones and in broad shallow surrounding areas. Such changes are seen by analysis of repeating events (e.g. Peng & Ben-Zion 2006; Rubinstein et al. 2007), deconvolution of borehole and surface data (e.g. Wu et al.2010; Nakata & Snieder 2011) and other techniques. The coseismic changes in the immediate vicinity of fault zones can approach tens of per cent, and are followed by logarithmic recovery (healing) in time (e.g. Sawazaki et al. 2006; Wu et al. 2009). Similar large reductions of seismic properties during brittle failures are seen in laboratory fracturing experiments (e.g. Gupta 1973; Lockner et al. 1977). High-resolution laboratory and seismic studies indicate that the healing processes already operate at high rates 1 s after the failures, and may reach asymptotic saturation at ∼104–105 s (e.g. Johnson & Jia 2005; Nakatani & Scholz 2006; Wu et al. 2009; Ben-David et al. 2010).
In recent years, noise-based imaging techniques have been developed considerably and used to derive properties of various structures and processes in the Earth (e.g. Shapiro et al. 2005; Brenguier et al. 2008; Roux 2009). This started with the study of Campillo & Paul (2003) on retrieving the Rayleigh wave tensor between a set of sensors in Mexico from the coda of seismic events. Subsequently, cross-correlation techniques rapidly moved to using ambient seismic noise, available at longer time intervals, instead of earthquake data. This created possibilities for imaging regions with limited seismicity and was proven successful over periods ranging from 0.1 to 100 s and corresponding spatial scales from 100 to 1000 km (Campillo et al. 2011, and references herein).
There are several attractive features of noise-based imaging, including, significantly, the ubiquitous operation of sources in essentially all space–time domains. However, since the noise sources are relatively weak and diffuse, they are easily masked by earthquakes and other large episodic sources, which have to be removed in standard noise-processing techniques. The large numbers of aftershocks in the days to weeks following large mainshocks create space–time windows near large earthquake rupture zones where noise data cannot be used. The same may hold during time windows associated with accelerating foreshock sequences. Current noise-based techniques that focus on temporal changes in the coda of noise correlations require stacking of relatively quiet periods (e.g. Campillo 2006; Bensen et al. 2007). This, and the available frequency range of typical noise-based data, precludes imaging properties of structures and dynamic processes that are highly localized in space and time (e.g. sub 100 m and sub 1 d). The high near-fault seismicity rates, combined with strong attenuation and various structural complexities, may explain why applications of noise-based monitoring to fault zone environments have been limited so far (e.g. Sens-Schönfelder & Wegler 2006, Brenguier et al. 2008; Roux et al. 2011; Hillers et al. 2013).
On the other hand, the numerous aftershocks in the early periods following the occurrence of mainshocks generate many waveforms that are ideal for studying rapid fault zone processes, in combination with noise-based correlation waveforms that are adequate for studying structures and processes at other time intervals. Utilizing both types of data is clearly highly desirable. Analysing waveforms that include both earthquake and noise signals with similar processing techniques can provide a self-consistent methodology that would avoid the need to merge results obtained by different methods. As a step towards this goal, we analyse in this paper correlated waveforms generated by earthquakes using noise-related processing techniques.
The study employs seismic waveform data generated by aftershocks of the 1999 Mw 7.4 İzmit and Mw 7.1 Düzce earthquakes along and around the Karadere-Düzce branches of the North Anatolian Fault (NAF). We develop and apply a new correlation-based methodology for monitoring fault zone properties, which preserves amplitude information and focuses on evolution of coherency. The examined earthquake waveforms are recorded by a pair of stations, one located within the core damage zone of the NAF and the other ∼400 m away from the fault. The correlation functions are realigned to overcome clock problems that may commonly exist in rapid deployment data. The results show strong drop of coherency at the time of the Duzce mainshock, followed by oscillatory recovery over 4 d likely associated with fluid motion in the highly damaged fault zone rocks.
2 TECTONIC ENVIRONMENT AND DATA
The NAF system is the largest geological structure in Turkey, extending for over 1500 km from the Karliova triple junction where it intersects with the East Anatolia Fault. The fault defines the northern boundary of the Anatolian block that is extruded westward under the combined effect of Arabia–Asia convergence to the east and southwestward extension related to the Aegean subduction rollback in the west (e.g. Jackson & McKenzie 1988). This regime has been active for ∼10 Myr, has accumulated ∼80 km of right-lateral offset and has a current rate of ∼2.4 cm yr–1 in northwestern Turkey (e.g. Reilinger et al. 1997; Armijo et al. 1999; Sengör et al. 2005). During the past century, the NAF had a sequence of large ruptures that unzipped nearly the entire fault from east to west (e.g. Stein et al. 1997; Sengör et al. 2005). The sequence culminated with the 1999 August 17, Mw 7.4 İzmit earthquake and the following 1999 November 12, Mw 7.1 Duzce earthquake.
A week after the Mw 7.4 İzmit earthquake, a 10-station seismic network was deployed along and around the Karadere-Duzce branch of the NAF (Seeber et al. 2000; Ben-Zion et al. 2003). This network operated for ∼6 months and recorded approximately 26 000 events, including aftershocks of the İzmit earthquake, the Mw 7.1 Düzce earthquake and both foreshocks and immediate aftershocks of the Duzce event. The seismic stations had REFTEK recorders and three-component L22 short-period sensors with a sampling frequency of 100 Hz. The horizontal location errors of the earthquakes are generally less than 1 km around the centre of the network and 1–2 km near the margins. The vertical location errors are somewhat larger. Additional details on the seismic experiment and basic properties of the data are given by Seeber et al. (2000) and Ben-Zion et al. (2003). The recorded seismograms were analysed previously in a set of papers focusing on various seismic signatures of rock damage (Ben-Zion et al. 2003; Peng & Ben-Zion 2004, 2005), temporal changes of material properties (Peng & Ben-Zion 2006; Wu et al. 2009) and earthquake source properties (Lewis & Ben-Zion 2007; Yang et al. 2009; Wu et al. 2014).
In this paper, we use seismic waveforms recorded by a station pair VO and FP having short spacing near the fault (Fig. 1). Station VO is located within the rupture zone of the İzmit earthquake on a shutter ridge with high relief composed of gouge and slope debris. Station FP is ∼400 m off the fault in soil near bedrock with medium relief.
3 A NEW MONITORING TECHNIQUE BASED ON EARTHQUAKE WAVEFORMS
As mentioned in the introduction, the increasing seismic activity immediately after (and sometimes also before) moderate and large earthquakes prevents the use of ambient noise for a considerable period of time (e.g. one or a few days before/after the event). During this time interval, seismic traces are dominated by waveforms generated by a large number of weak/strong earthquakes that require a specific processing when used in cross-correlation algorithms. An alternative method developed here is to use the waveforms generated by each earthquake to construct cross-correlations between various station pairs in a seismic network. Unlike coda wave interferometry (Poupinet et al. 1984; Snieder et al. 2002) that uses only the coda of similar events observed at a given station, the developed method involves computing two-point correlation functions of entire waveforms, including both ballistic and scattered waves, generated by all events that are recorded at a given pair of stations.
The station separation distance should be short since the frequency content of earthquakes is high and the spatial coherence decreases fast with distance and frequency. The spatial distribution of the seismic events or the scattering properties of the medium should be sufficiently diffused to allow extraction of information from the cross-correlation processing.
Fig. 1 displays the locations of the two seismic stations (white letters) and epicentres of ∼10 000 earthquakes used in this study (coloured dots and small stars). Station VO is located within a vertical seismic trapping structure (waveguide) that is ∼75 m wide and has ∼50 per cent velocity reduction relative to the surrounding rock (Ben-Zion et al. 2003). The employed events occurred in a 100-d period that includes the 1999 November 12 [Julian Day (JD) 316], Mw 7.1 Duzce earthquake (focal mechanism beach ball). Stations VO and FP are 410 m apart and both recorded all the seismic events during the examined period. These stations were selected because their strongly contrasting local environments (within and outside the core damage zone) and relatively short separation can be used to monitor changes of medium properties before and after the Duzce earthquake (Peng & Ben-Zion 2006; Wu et al. 2009).
Fig. 2(a) shows two waveforms at stations VO and FP generated by a seismic event with local magnitude ML = 2.8 occurring on JD 280. The signals have similar amplitudes and frequency content and are delayed in time in relation to the distances between the source and receivers. Figs 2(b)–(d) show corresponding bandpass-filtered waveforms in different frequency bandwidths around 3, 1.5 and 0.75 Hz. As expected, the resemblance between the signals decreases with increasing frequency, reflecting the high structural heterogeneity of the fault zone environment at smaller scales.
The collection of correlation functions computed in the (0.5–1 Hz) frequency band for each of the ∼10 000 events between JD 280 and 380 is plotted in Fig. 3(a). All correlation functions exhibit a dominant peak whose arrival time should be associated with the traveltime difference between the pre-processed signals. However, a close inspection reveals jumps and drifts between the clocks at stations V0 and FP that can mask changes in the properties of the seismic events (locations, magnitude and focal mechanism) or changes in the properties of the medium. When each correlation function is realigned with respect to its maximum at time t = 0 (Fig. 3b), the pattern appears very stable and uniform, suggesting that the correlations between the two nearby sensors are largely independent of the event characteristics (location, size and mechanism).
The comparison between Figs 3(a) and (b) point to clock uncertainties between the two stations over the course of the recordings. Such problems prevent the use of time information as classically done in passive monitoring from continuous records (e.g. Brenguier et al. 2008). On the other hand, the amplitude information of the correlation functions may still carry signatures of physical changes in the propagation medium. The purpose of this study is to monitor fault zone properties through amplitude rather than time information extracted from a collection of events-driven correlation functions. The average correlation function of the realigned individual correlations resembles an autocorrelation function (Fig. 3c). This is consistent with the short distance between the two stations (∼400 m) and the relatively low frequency content of the pre-processed signals. Since station VO is in the core damage zone of the NAF, one expects the cross-correlation function to be highly sensitive to local damages at this location, with station FP serving as a (relatively unperturbed) reference. This is consistent with previous results on evolving seismic properties at these locations (Peng & Ben-Zion 2006; Wu et al. 2009).
The new methodology, geared to analysis of earthquake waveforms recorded during aftershock sequences of large events, should account for various clock problems expected in rapid deployments, while capitalizing on the high level of coherence observed on the correlation function for each earthquake recorded by different stations (Figs 3a–c). Fig. 4 shows the distributions of coherence values (calculated as the maximum of the correlation functions) computed in different frequency bandwidths for the entire set of used events. As expected, the average coherence decreases with increasing frequency. The shapes of the distributions follow generally the Poisson law, indicating an overall statistical independence of results associated with different events. However, the detailed analysis done below reveals clear evolution of properties and memory effects that start at the time of the Duzce mainshock. To maintain a high level of coherence for a large set of events, we focus the amplitude analysis of the correlation functions at relatively low frequency around 0.75 Hz. This choice implies an average wavelength of a few hundred metres that is the typical length scale at which the fault zone is monitored in this work.
The monitoring of the fault structure is investigated with data of stations VO and FP during 45 d before and after the JD 315 when the Duzce mainshock occurred. Fig. 5 exhibits the temporal evolution of the number of events recorded at stations VO and FP per 2-h sampling interval. The large increase in the number of events at day 315 represents the beginning of the aftershock sequence of the Duzce mainshock. The time recovery follows overall a classical Omori's law (e.g. Omori 1894; Utsu et al. 1995; Zaliapin & Ben-Zion 2013) with rapid decay over ∼15 d. Note that more than 100 events are recorded for every 2-hr interval in the few days that follow the earthquake. Taking each event as an independent and stable source, as seen from the correlation pattern in Fig. 3(b), this represents a very dense temporal sampling of the physical processes that occur in the fault zone just after the mainshock. The smaller peak a few days earlier is associated with a cluster of events not very close to the Duzce mainshock hypocentres (Wu et al. 2014).
Fig. 6 shows the time evolution of the coherence level for different average periods and different frequency bands. In practice, the collection of correlation functions is calculated for every event over a 0.5-Hz frequency bandwidth and a set of 16 central frequencies between 0.55 and 1.15 Hz. The realigned correlation functions are then averaged over N = 25, 50 or 100 consecutive events for smoothing purposes (Figs 6a–c). The coherence levels of the obtained average correlations are measured as a function of time as a first attempt to monitor the fault damage zone. The time/date associated with a given set of N events is the average time/date of the events. The running average over N consecutive events produces a sampling that is not uniform with time; the period just after the mainshock where the largest changes are expected is monitored with shorter time intervals. A strong decrease of the coherence level is observed at low frequency at the time of the Duzce earthquake (JD 315), with possible recovery over ∼4 d. However, the natural decrease of the coherence level with increasing frequencies and numerous temporal fluctuations over the scanned period do not facilitate detailed interpretation of results.
To obtain more refined observations, we compare each of the N-averaged correlations to the reference wavelet in Fig. 3(c) that corresponds to the average correlation function over the entire set of events. The correlation coefficient is calculated as the maximum correlation between the N-averaged correlations and the reference wavelet. Such a comparison is commonly done in noise-based processing, for example, when traveltime fluctuations are calculated on the coda of correlation functions (Hadziioannou et al. 2011). Figs 7(a)–(c) show the frequency and time-dependent correlation coefficient for N = 25, 50 and 100. The correlation coefficient is generally close to 1, since all correlation functions resemble each other, but the Duzce mainshock clearly manifests itself as a 1–2 per cent drop in the correlation coefficient. The size of the observed effect on the coherency of the correlation is one order of magnitude larger than the traveltime change (∼0.01–0.1 per cent) measured in the coda of the noise correlations around the San Andreas Fault after the Parkfield earthquake (Brenguier et al. 2008).
The frequency dependence of the mainshock signature in Fig. 7 exhibits high sensitivity around 0.8 Hz. Fig. 8 displays the time-evolving correlation coefficients at 0.8 Hz for N = 25, 50 and 100. The recovery time after the mainshock is ∼4 d, as seen also in Figs 6 and 7, which is much shorter than the timescale of the Omori's law. This indicates that the temporal distribution of the events does not obscure the timescale of the physical process that occurs in the fault zone following the Duzce earthquake. We note that there are only 13 aftershocks with ML > 3.75 in the study period (small stars in Fig. 1), with most occurring in the first 2 hr after the 7.1 Duzce mainshock. To confirm that the fluctuations of the correlation coefficient during the recovery time are not dominated by the evolving flux of seismic energy, we examine the evolution of several quantities related to the generated seismic energy in the region.
Using the empirical quadratic scaling relation of Ben-Zion & Zhu (2002), we convert the magnitude of each event to the scalar seismic potency (moment/rigidity) P0, which is proportional to the source energy, and inspect the temporal evolution of P0/r2 from all events, where r is the distance from each hypocentre to the midpoint between stations VO and FP. The results display noisy fluctuations with no clear pattern in the ∼4 d after the Duzce mainshock. We also examine the temporal changes of P0 for all events in radii of 5–15 km around station VO. The results show a rapid decay over ∼1 d that does not correspond to the recovery pattern in Fig. 8. Finally, we compute the observed radiated energy at station VO filtered in the frequency bandwidth of interest (0.5–1.0 Hz) and find again ∼1 d decay time. The results indicate that the observed fluctuations of correlations between JD 316 and 320 in Figs 6–8 are not dominated by a reactivation process associated with the ongoing aftershocks. We also note that the evolving complexity of the waveforms with many more overlapping events in the early aftershock sequence cannot produce the observed pattern, since the performed correlations of average narrow-band seismograms account for different waveform complexities. As discussed in the next section, the fluctuations likely represent motion of fluids within the core damage zone below station VO.
4 DISCUSSION
Seismic waves can be used to monitor the evolution of rock properties, stresses and fluids in fault zones and other important environments, such as volcanoes, oil reservoirs and subsurface storage facilities. The monitoring of changing properties can provide fundamental information on processes leading to and following material failure. In this work, we developed a new monitoring methodology based on evolving coherency of correlation of earthquake waveforms. The technique is applied to data recorded by a pair of stations (VO and FP) ∼400 m apart in the immediate vicinity of the NAF during a time window that includes the 1999 Mw 7.1 Duzce earthquake (Fig. 1).
In contrast to classical processing of ambient noise, the present methodology does not require clipping the waveforms or selecting particular time windows. Since each waveform produces high coherence at different stations, the time step is only limited by the rate of events, which is very high after large earthquakes. To overcome clock problems, the individual correlation functions are realigned on the main peak as was done in previous analysis of the data (e.g. Peng & Ben-Zion 2006). This leads to a loss of traveltime information that precludes measuring traveltimes and monitoring seismic velocities. However, changes in the coherency of successive correlation functions reveals clear co- and post-Duzce changes in the wave propagation properties of the sampled medium between the stations.
Station VO resides within the core damage zone of the NAF above a waveguide structure with width of ∼75 m and ∼50 per cent velocity reduction relative to the surrounding rock (Ben-Zion et al. 2003). Previous analyses of spectral ratios of data recorded at VO and FP and waveforms generated by clusters of repeating earthquakes indicated >20 per cent co-Duzce drop of seismic velocity in the top few hundred metres of the crust underneath station VO (Wu et al. 2009), followed by log(t) recovery that was not complete after 3 months (Peng & Ben-Zion 2006). This study shows a clear co-Duzce drop of coherence in the wavefield between VO and FP at a narrow frequency band centred around 0.8 Hz. Note that no change is observed at higher or lower frequencies. The observed localized changes around 0.8 Hz are likely associated with resonance effects of the fundamental mode of trapped (Love or Rayleigh type) waves in the core damage zone. Estimating the depth section sampled by the fundamental Rayleigh mode requires a detailed knowledge of the seismic velocity profile in the shallow crust, which is not available in the study area. However, assuming very low velocities in the top few hundred metres, based on measurements in California boreholes (e.g. Ely et al. 2010; Taborda & Bielak 2013), suggests that the changes occur in the top 100–300 m as in agreement with previous studies.
The strong reduction of coherency (up to 1–2 per cent of the correlation coefficient) is followed by an overall recovery on a timescale of 4 d (Fig. 8). The much shorter recovery timescale compared to that observed by Peng & Ben-Zion (2006) points to a process that is different from a plain recovery of rock damage. The observations likely reflect post-seismic fluid motion in the highly damaged trapping structure below station VO that was initiated at the time of the Duzce event. While quantitative data on fluid flow near station VO are lacking, the Duzce earthquake generated clear visible fluid effects in the fault zone environment (M. Bouchon, personal communication, 2013). The large coseismic changes of seismic velocity below station VO (Wu et al. 2009) had to be accompanied by increase of permeability by orders of magnitudes (e.g. Wibberley & Shimamoto 2003; Fialko 2004), along with coseismic reduction and increase of stress inside and outside the fault zone, respectively. This would have led to rapid fluid motion into the damaged fault zone rocks just after the Duzce event that produced the observed large drop of coherency, followed by fluid flow away from the damaged fault zone rocks leading to recovery of coherency. The overall recovery timescale of 4 d can be explained in terms of fluid diffusion over a length scale of ∼100 m (width of the highly damaged rocks below VO) and permeability of 0.03 m2 s–1, which is reasonable for highly fractured rocks (e.g. Lockner et al. 2009). The oscillatory recovery process may reflect a complex interplay between healing of the damaged rocks and changes of pore pressure, as seen in recent model simulations that incorporate rapid changes of fluid flow and damage rheology (Shalev et al. 2013). Further interpretation of other fluctuations in the obtained results should be done with care, keeping in mind the evolving averaging timescales (associated with given event numbers).
5 CONCLUSIONS
This work describes a new methodological procedure for monitoring fault zone properties using correlations of earthquake waveforms at a pair of stations. The method benefits from the recording of a large number of seismic events right after the occurrence of mainshocks, which can finely sample the temporal changes of properties between the stations. The monitoring is based on the amplitude of correlations of each seismic waveform recorded by nearby stations close to the fault. At low frequencies, the coherence level is high and the correlation function appears to be very stable and independent of the particular properties of the used seismic events. Each correlation function of data between nearby stations can be interpreted as an autocorrelation that would measure the wavelet energy at the receiver located in the damaged zone. The time evolution of the correlation coefficient of groups of events with a reference wavelet shows a significant co-mainshock change that recovers with some prominent oscillations over a timescale of a few days. The reduction and recovery of coherency are interpreted to reflect the interplay between generation and healing of rock damage and permeability in the fault zone combined with fluid motion. The method can be generalized to combine analysis of both earthquake and noise waveforms. This will be done in a follow-up work.
The paper benefited from useful discussions with Michel Campillo and constructive comments by Toshiro Tanimoto, an anonymous referee and Editor Mike Ritzwoller. The study was supported by grant EAR-1141944 of the National Science Foundation.