Introduction

Two earthquakes exhibiting very similar waveforms (doublet) are representative of the subsequent activation of a source whose mechanism and location do not change over time. In a doublet, while ballistic (i.e., P- and S-waves) arrivals match, the coda can lose coherency, since multiply scattered waves are sensitive to small changes in the medium. The earthquake coda, composed of multiply scattered waves, samples the medium more effectively than the direct (ballistic) arrivals. For identical co-located sources, acting at different times, any observed difference in waveforms is related to a change in the elastic properties of the medium. The coda wave interferometry technique (CWI) uses multiply scattered waves to detect temporal changes in a medium by using the medium as an interferometer, i.e. comparing scattered waves that sample the medium at different times. Following these principles, this app applies the Coda Wave Interferometry technique in order to estimate possible velocity changes between pairs of similar earthquakes.

Method

The app examines a list of seismograms from the same component of the same station. At first, all the independent pairs of seismograms are cross-correlated, in order to identify doublets of repeating earthquakes. A doublet is defined as a pair of seismograms whose correlation is larger than a threshold input parameter. For any identified doublet, the app then applies the CWI method. This consists in evaluating the cross-correlation function (CCF)  at subsequent lapse times T along the coda waves of the two similar earthquakes. For any window position, the time lag τ corresponding to the maximum of the CCF is a measurement of the mean perturbation in traveltimes of one event with respect to the other. Therefore, the behaviour of τ with lapse time along the coda allows identifying the source of the perturbation. In fact, when the perturbation is generated by a (small) displacement of scatterer or source location, τ oscillates around 0 (Snieder, 2006). If, instead, the perturbation is due to a homogeneous velocity change over the medium sampled, then the time separation between the two waveforms will increase (decrease) linearly for later times in the coda. As a first approximation, this trend may be fitted by a linear regression curve, thus the relative velocity change measurement in terms of percentage variation corresponds to the slope of  τ / T changed in sign. 

Internal Flux

This is the description of the code's internal flux. All these steps are transparent to the user; information on how to run the App are in the next Section 4. Terms in bold refer to the input parameters required by the App (see Table 1)

1. Reading seismograms from the same component of the same seismic station;

2. Preconditioning: trend and mean removal, tapering;

3. Band-pass filtering over the [fmin, fmax] frequency band;

4. Calculation of the N(N-1)/2 cross-correlations between all the independent seismogram pairs;

5. Select ion of the NG seismogram pairs for which the correlation coefficient is larger than the threshold parameter;

6. For each pair of similar seismograms:

   • Alignment of the two time series according to the delay time estimated by crosscorrelation;

   • Starting at tstartcwi seconds after largest amplitude, and spanning a total duration of tlencwi seconds, performing of the CWI loop for subsequent lapse times T:

     1. Calculation of the CCF over a window of length wl second;

     2. Interpolation of the CCF by fitting the 3 samples encompassing the CC maximum with a parabola;

     3. Evaluation of the lag time τ from the maximum of the interpolating function, with a precision of 1/10 of the original sampling rate;

     4. Storing the time lag τ and corresponding maximum of the CCF;

     5. Shifting the correlation window wl by 50% of its length;

     6. Goto point (a)

7. The τ vs T relationship is then fitted to a straight line using a robust, L1-norm approach implemented through the IRLS algorithm (Aster et al., 2013).

8. The percent velocity change ∆v = - τ / T is finally evaluated from the slope of the fitting line.

Figure 1 illustrates a pictorial view of the procedure, with specification of some of the relevant parameters.

Figure 1. Main parameters mentioned in the text.

Explaination: Ts-Tp is the delay time between the S and P waves arrival. tstartcwi is the time lag after the maximum amplitude (usually associated with the S-wave arrival) at which starting the CWI analysis. wl is the length of the time window for correlation estimates. tlencwi is the time interval (in s) after tstartcwi over which extending CWI analysis.

Step by Step

Obtaining the SAC files

Selection of the episode of interest

From the list of episodes on the EPISODES Platform, select the one of your interest (in this case, the Song Tranh reservoir). You'll be redirected to the page of the selected episode (Figure 2). Click on the link 'Event Related Waveform' within the data section.

Figure 2. Episode view.

Selection of the recordings

All the data files (seismic recordings) related to this episode will be listed (Figure 3). Select the ones of your interest (in this case, the first 10 files from the Song Tranh episode).

Figure 3. Listing of files related to the episode.

Adding selected files to your workspace

Once files have been selected, click on the 'ADD SELECTED TO WORKSPACE' button, at the bottom left of the window. You'll be prompted for specifying the directory where files should be copied (Figure 4). Select the appropriate directory. (It is suggested to make sure that the directory where the User is putting data is empty.)

Figure 4. Choice of workspace directory to place the files.

Format conversion and data checking

Once all the data files have been copied into the selected destination, you must proceed to the 'SEED to SAC converter' application and click on 'ADD TO WORKSPACE' button to add converter to the User's workspace. You will be asked for a location where creating the folder associated with the application (it has a default name: 'SeedtoSac'). The application starts automatically with the following input form (Figure 5):

Figure 5. The "Seed to SAC converter" application input form.

By clicking on the 'Select Files' button, you'll be asked for selecting data files to be converted (Figure 6). 

Figure 6. Selection of files to be converted.

Select the 'Files from chosen directory' option, and then select the directory where the SEED data files to be converted reside (in this example, it is the DATA directory). Again: this directory must contain ONLY the data files you want to finally process.

Run the app by clicking on the 'RUN' button (Figure 7). You should see a 'Status: Running' message, and the progressive growing of output file names within the 'OUTPUTS' window at the bottom of the page. At the end of the process, all the newly converted data files will be stored in the 'SeedToSac' directory.

Now it's a good moment for checking your data. From the directory tree at the left, go to the 'SeedToSac' directory and click on one of the newly-converted data files. The corresponding seismogram will be displayed (Figure 7). Here you can get information that are necessary for the subsequent selection of input parameters to the the Coda Wave Interferometry, such as sampling rate, seismogram duration, etc.

Figure 7. Display of the selected seismogram.

Running the Coda Wave Interferometry

After adding the 'Coda Wave Interferometry detection of velocity changes' application to the workspace, the following window appears on the screen (Figure 8):

Figure 8. Input window of the application.

Click on the 'SELECT FILES' button. A new menu window appears; select the 'Individual files' option, and then select the SAC files of interest (at least two are required, they should be from the same station and showing the same channel). Alternatively, you may want to choose the 'Application Output' or 'Files from chosen directory' option and select the 'SeedToSac' folder name. This directory must contain ONLY the data files you want to process.

Finally, you have to select the rest of the parameters as explained in the Table reported in the next Section.

Understanding the Input Parameters

Caveats

• CWI window length wl cannot be longer than seismogram length;

• CWI window length wl cannot be longer than the time interval tlencwi selected for the analysis;

• The ending time of the analysis (tlencwi + tstartcwi) cannot exceed the ending time of the seismogram;

• A reasonable choice for the CWI window length wl is twice the low-pass period, i.e. wl = 1 / fmin;

• The low-frequency corner of the filter fmin must satisfy the relationship : fmin > 1/T, where T is the time length of the seismogram (e.g., for a 10-s-long seismogram, fmin must be greater than 0.1 Hz);

• The high-frequency corner of the filter fmax must be lower than the Nyquist frequency : fmax < 0.5 Fs, where Fs is the sampling frequency of the seismograms (i.e., for a seismogram sampled at 100 Hz, fmax must be lower than 50 Hz);

• The input seismograms must be in SAC format, with the KSTNM (station name) and KCMPNM (component name) header field properly filled-in;

• The input seismograms must be from the same component of the same station.

The main output is a text file (cwi.results.txt) containing a row for each pair of seismograms which have been analysed with the following data:

ind1 ind2 name1 name2 tstart1 tstart2 deltaV errdeltaV

ind1= index of the first seismogram, with respect to the original list

ind2= index of the second seismogram, with respect to the original list name1 = name of seismogram 1 name2 = name of seismogram 2

tstart1 = starting time of seismogram 1

tstart2 = starting time of seismogram 2

deltaV = percent velocity change derived from CWI errdeltaV = error in percent velocity change estimate

In addition, two set of figures are produced:

• A record section (png image, named cwi.recordsection.png with all the analysed seismograms prior to correlation thresholding (Figure 9)

• A png picture for each of the seismogram pair which passed the correlation threshold, illustrating the linear fitting of the delay time vs lapse time (named cwi.outfig.???.png); see Figure 10.

Figure 9. The Record Section plot, illustrating all the seismograms selected for the analysis.

Figure 10. Results from the CWI analysis.

Top: the two earthquakes selected for the processing. The dashed area marks the time interval selected for the evaluation of the lag time vs lapse time relationship. Middle plot: Time lags as a function of lapse times. The continuous line is the line bestfitting the data in the L1-norm sense. The two dashed lines are the error bounds. Bottom plot: maximum of the cross correlation function evaluated at subsequent lapse times.

References

Aster, R., B. Borchers and C.H. Thurber (2013). Parameter Estimation and Inverse Problems (Second Edition). Elsevier, 360 ppgg.

Martini F, CJ Bean, G. Saccorotti, F. Viveiros, N. Wallenstein, (2009). Seasonal cycles of seismic velocity variations detected using coda wave interferometry at Fogo volcano, São Miguel, Azores, during 2003–2004. Journal of Volcanology and Geothermal Research 181 (2009) 231–246

Pandolfi D, Bean CJ, Saccorotti G (2006). Coda wave interferometric detection of seismic velocity changes associated with the 1999 M=3.6 event at Mt. Vesuvius. Geophysical Research Letters, vol. 33, ISSN: 0094-8276, doi: 10.1029/2005GL025355

Snieder, R. (2002), Coda wave interferometry and the equilibration of energy in elastic media, Phys. Rev E, 66, 046615, doi:10.1103/PhysRevE.66.046615.

Zaccarelli L., Pandolfi D., Bianco F., Saccorotti G, Bean C. J., Del Pezzo E., (2009). Temporal changes in seismic wave propagation towards the end of the 2002 Mt Etna eruption. Geophysical Journal International, vol. 178, p. 1779-1788, ISSN: 0956-540X, doi: 10.1111/j.1365-246X.2009.04219.x

Back to top