Toward source region tomography with inter-source interferometry: Shear wave velocity from 2018 West Bohemia swarm earthquakes
Tom Eulenfeld
physics.geo-ph


Towards Backward-Compatible Representation Learning

    
Yantao Shen[1]    
Yuanjun Xiong    
Wei Xia    
Stefano Soatto 


AWS/Amazon AI

ytshen@link.cuhk.edu.hk,   {yuanjx, wxia, soattos}@amazon.com

    December 30, 2023
=====================================================================================================================================================



 b]l
An edited version of this paper was published by AGU. Copyright 2020 American Geophysical Union.

T. Eulenfeld (2020), Journal of Geophysical Research: Solid Earth, 125, e2020JB019931,

doi: http://dx.doi.org/10.1029/2020JB01993110.1029/2020JB019931.




The concept of seismic interferometry embraces the construction of waves traveling between receivers or sources with cross-correlation techniques. In the present study cross-correlations of coda waves are used to measure travel times of shear waves between earthquake locations for five event clusters of the 2018 West Bohemia earthquake swarm.
With the help of a high quality earthquake catalog, I was able to determine the shear wave velocity in the region of the five clusters separately. The shear wave velocities range between 3.5km/s and 4.2km/s. The resolution of this novel method is given by the extent of the clusters and better than for a comparable classical tomography. It is suggested to use the method in a tomographic inversion and map the shear wave velocity in the source region with unprecedented resolution.
Furthermore, the influence of focal mechanisms and the attenuation properties on the polarity and location of the maxima in the cross-correlation functions is discussed. The intra-cluster ratio of P wave to S wave velocity is approximately fixed at 1.68.
Key points:


  * Shear wave travel times between earthquake locations from cross-correlations of coda waves

  * Spatial variability of seismic velocity at region of 2018 earthquake swarm

  * Enhanced spatial resolution compared to classical tomography



Keywords: inter-source interferometry, inter-event interferometry, shear wave velocity, coda waves, cross-correlation, swarm earthquakes, vp/vs ratio




§ INTRODUCTION




The  concept of seismic interferometry is most often applied by cross-correlating a coherent or incoherent seismic signal at different stations. In the limit of a sufficiently random wave field, the cross correlation function converges towards the Green's function between the two receivers <cit.>. 
In this inter-receiver setting, the Green's function or the cross-correlation function for arbitrary time-invariant noise sources can be used to monitor the wave velocity between and near the receivers <cit.>. 
In the reverse, seismic interferometry can be applied to use earthquakes as virtual receivers at depth. In contrast to inter-receiver interferometry this approach cannot be easily used for monitoring, but has the advantage that the medium between and near the event locations is probed.

A number of studies used coda cross-correlations of closely spaced earthquakes to determine the inter-event travel time from the decorrelation of the coda wave field of two earthquakes <cit.>. For this method, the earthquakes need to be separated by less than the dominant wavelength. Additionally, both events need to have the same moment tensor and the scatterer density is assumed to be constant in all directions. <cit.> show that under these conditions the decorrelation of the coda wave field is proportional to the earthquake separation, because both are proportional to the variance of travel time perturbations associated with different wave paths. The proportionality factor is different for different types of sources. 



<cit.> used spatial reciprocity <cit.> to construct virtual pseudo-noise waveforms at each earthquake location by summing coda waves registered at different stations. These virtual recordings at different earthquake positions were then cross-correlated to obtain inter-event shear wave travel times.
<cit.> further demonstrated the principle of inter-source interferometry by comparing real and virtual recordings of the Sichuan earthquake.
Different from <cit.> they used the full seismic recording to mainly reconstruct surface waves and therefore can restrict the used stations to the stationary phase zone.
<cit.> stated that this approach was superior in their application because more energy was used by the cross-correlation and the virtual recordings therefore were better defined.
<cit.> selected recordings in the stationary phase zone of body waves propagation between deep earthquakes and extended this inter-source interferometry concept to P and S body waves.

Inter-source interferometry as in <cit.> and <cit.> using direct body waves or surface waves recorded by stations located in the stationary phase zone of the two earthquakes is in the following referred to as method B.


<cit.> used coda cross-correlograms to determine the S wave travel time between aftershocks of the Lushan Earthquake for an event relocation procedure (referred to as method C). <cit.> used autocorrelograms of coda waves to determine the event depth from the travel time of waves passing though the event location after a reflection at the surface.
Figure <ref> compares the three introduced methods of inter-event seismic interferometry.








In this article I further explore method C and apply the concept of inter-source interferometry by correlation of scattered coda waves of different earthquakes. It is assumed that event locations are separated by more than a wavelength. Therefore only waves scattered in the stationary phase zone behind the two source locations interfere constructively in the cross-correlation function. The wave path are indicated in figure <ref> in the right panel: Waves radiate from one event through the other event location. Waves from both events are scattered for the first time in the stationary phase zone and are recorded in the coda window with a time difference corresponding to the inter-event travel time. Waves from both sources scattered for the first time in a region outside of the stationary phase zone interfere destructively in the cross-correlation function. Therefore the maxima in the cross-correlation correspond to the inter-event travel time. This concept is basically the reciprocity of inter-station interferometry using earthquake coda waves <cit.> and equally applicable for isotropic scattering and the non-isotropic (forward) scattering expected in this study. 
Because more energy is radiated as S waves than as P waves <cit.> the coda will be dominated by waves originating initially from S waves. Therefore the focus is put on the determination of inter-event travel times of shear waves.


Method C is applied to the 2018 earthquake swarm near Nový Kostel in the Czech-German border region West Bohemia/Vogtland. The area  is of great interest as it hosts several signs of geodynamic activity such as CO_2 degassing, Quaternary volcanoes and earthquake swarms <cit.>.
The topography of West Bohemia and Vogtland is shown in figure <ref> together with WEBNET seismic stations and earthquake epicenters between the years 2000 and 2020.

Travel times of seismic waves between earthquake locations are of great interest as their knowledge allows to estimate the seismic velocity in the source region with a resolution determined by the earthquake distribution. Given the high event density in the volume of the swarm earthquakes this allows for a spatial resolution that is superior to conventional seismic tomography.
<cit.> developed a double difference arrival method to determine the ratio of P wave velocity to S wave velocity () for a region embraced by a cluster of earthquake.

<cit.> applied a similar method to the West Bohemia earthquake swarms of the years 1997, 2000 and 2008 and found a time-dependent between 1.38 and 1.70.
<cit.> showed that the inside different clusters of the 2014 earthquake swarm ranged from 1.59 to 1.73.


In this study the S wave velocity in the source region is determined directly from the travel time observations using a high quality earthquake catalog.
This approach is different from previous studies of inter-event interferometry which focus on the relocation of earthquakes.


When dealing with coda waves, the scattering properties of the medium are of interest. Therefore with the help of the Qopen software package <cit.> scattering and intrinsic attenuation properties of shear waves were estimated for a homogeneous scattering medium in the half space.
The transport mean free path, defined by the average length in which the wave forgets its initial direction due to single or multiple scattering, was determined at approximately 130km with the data set used in this study for a frequency range between 10Hz to 20Hz. The corresponding transport mean free time is 36s for a shear wave velocity of 3.6km/s. The coda waves with a travel time up to 50s used in this study are therefore multiply scattered, but not in the diffusion regime. The intrinsic absorption length is estimated at 55km for a frequency range between 10Hz to 20Hz.
Similar values were obtained by <cit.> and <cit.> for this region with different data sets.










§ DATA AND METHOD






The double difference catalog of the 2018 swarm used in this study was compiled by Martin Bachura <cit.>. It consists of approximately 1000 earthquakes with local magnitudes larger than 1.3. The uncertainties in origin locations are 50m relative to each other.
For the purpose of the present study, the catalog is divided into 5 clusters named a-e with earthquakes larger than magnitude 1.8 delimitated and separated by the times , , , , ,  (all in Coordinated Universal Time, UTC). The clusters are plotted in figure <ref> in map view and depth sections. The small cluster south of the region indicated with the box in figure <ref> is excluded from the analysis.

Cluster a started the 2018 swarm activity at a depth of 9km to 10km. In the following days, the activity migrated to the north (clusters a-c). On May 21 earthquakes (cluster d) started rupturing a region south of the first cluster a at lower depth (around 7km). After June 19 the activity faded out with cluster e.
The event clusters are defined in time, but are with few exemptions equally separated in space as can be noticed from figure <ref>.

Focal mechanisms of earthquakes larger than magnitude 3 are strike-slip mechanisms with a normal faulting component aligned with the Mariánské Lázně fault zone <cit.>.


The waveforms of the earthquakes were recorded on the Czech WEBNET stations with 250Hz sampling rate. For this study, data from 9 WEBNET stations are used (compare figure <ref>).

Waveforms are bandpass filtered in the frequency range 10Hz to 40Hz. Coda time windows are selected from 1s after the S pick to 50s after the origin time. A shorter time window is taken if waves from another event in the catalog interfere. In this case, the coda window ends at the P pick of the following event. When the amplitude falls below a threshold based on the noise level the coda window ends at this time. All waveforms are visually inspected and the time window is adapted if necessary due to earthquakes not present in the catalog but interfering with the coda waves.

To compensate for intrinsic attenuation data is normalized by division by the instantaneous amplitude (i.e. envelope). This normalization step, later referred to as time normalization, ensures that data inside the coda window is equally weighted.
The following processing steps are applied to pairs of earthquakes. 
Pairs with an inter-event distance larger than 1km are excluded due to the low fraction of scattered waves which travel the direct path between the events.

Data of two different earthquakes recorded by the same station and channel are aligned relative to origin times. The overlap of the coda windows defines the time window (t_1, t_2) used for the cross-correlation. Station-component combinations with time windows shorter than 10s are discarded. For each station and component the first signal s_1 from the shallower event and second signal s_2 from the deeper event are cross-correlated according to


    C t = ∫_t_1^t_2 s_1τ s_2τ - t dτ/(∫_t_1^t_2 s_1τ^2 dτ)^1/2(∫_t_1^t_2 s_2τ^2 dτ)^1/2 .




Only correlations between the components Z-Z, N-N and E-E are used; cross-component correlations did not show a peak with a similar height at the expected time interval.

For each event pair all cross-correlation functions for different stations and components are stacked together. To enhance the signal a phase-weighted stack of order 2 is used <cit.> and in the process each trace is weighted by the length of the coda window used for the cross-correlation. Figure <ref> displays examples of cross-correlations for two different event pairs. In the left panel <ref>a peaks corresponding to waves traveling from the location of one earthquake to the location of the other earthquake can be identified on each individual trace. Contrary, in the right panel <ref>b the peaks corresponding to the inter-event travel time only emerge after the stacking procedure.

The relationship between the lag time of the maxima resulting from waves traveling between the earthquake locations and the inter-event travel time is


    =  +  1 -  2 + δerr .


 is the sign function and 1 and 2 are the errors of the origin times of the first and second earthquake, i.e. the difference between catalog origin time and real origin time.
If spatial deviations in seismic heterogeneity are present around the two sources, the lag time of the maxima might slightly deviate from the travel time corresponding to the stationary point by the error δerr. Because of the high dominant frequencies used in this study, a relatively homogeneous scatterer distribution is expected in the stationary phase zone and this error is therefore neglected in the following.

When two peaks 1 and 2 corresponding to waves traveling between event locations in both directions are visible in the cross-correlation function as in figure <ref>b, the inter-event travel time and difference in origin time errors can be determined directly with


    =  1/2| 1 -  2| ,
     1 -  2     = 1/2 1 +  2 .


Equation <ref> can be used to check the quality of the earthquake catalog.
In the scope of this study, errors of origin times are ignored and the inter-event travel time is approximated by the lag time of a single maximum: ≈||.

In figure <ref>a the phase weighted stacks of coda cross-correlations of each event pair are linearly stacked in different bins of event distance and plotted versus event distance. A clear move-out with a velocity of around 3.6km/s indicates that the peaks in the cross-correlation functions are due to shear waves traveling between the earthquake locations and that the errors in origin time are much smaller than inter-event travel times.

For comparison, the same procedure is repeated for the direct S wave window ( 0s,  2s) relative to the S pick without time normalization (method B).
When using the direct waves, spurious arrivals in the cross-correlation functions only cancel out if the station distribution is sufficiently dense. This corresponds to the requirement of a sufficiently random wave field in the inter-receiver setting.

The usual procedure to eliminate spurious arrivals in method B is the exclusive selection of stations inside the stationary phase zone, but because of the governing geometry the stationary phase zone is relatively small at the surface and no stations or only a small amount of stations may be located within the zone. Therefore, cross-correlations at all stations and components are simply aggregated with a phase weighted stack of order 2. The linear stacks in bins of event distance are displayed in figure <ref>b.
Not a single peak, but rather a multitude of peaks at absolute lag times smaller than the expected inter-event travel time emerge.
This expected result limits the usability of method B compared to method C for the present application.



§ RESULTS








Maxima with positive or negative polarity are extracted from the phase weighted stack for each event pair. To avoid the definition of a noise time window, the noise level is arbitrarily defined as the absolute value of the 8th largest relative extremum in the stacked cross-correlation function. Figure <ref> displays the lag times of the maxima with positive or negative polarity for all event pairs. Several interesting features can be observed in this figure.

First of all, the data points resemble a move-out of around 3.6km/s, although a considerable amount is off target. Therefore, around 1850 of 2900 data points are selected for this study by enforcing a signal-to-noise ratio larger than 10. Due to this condition a large amount of outlying data points are discarded (gray points in figure <ref>). The observed move-out shows that most of the peaks arise due to the extraction of shear waves traveling between the earthquake locations. Note that, maxima corresponding to inter-event travel times of compressional waves are not observed.
This indicates that initial P waves are less likely be recorded in the S wave coda than initial S waves which is consistent with the high energy ratio of radiated S to P  waves for earthquakes (see section <ref>).
The maxima therefore reflect the inter-event travel times of shear waves.

This is not the case for small event distances for which the maxima in the cross-correlation function are located at smaller lag times than expected from the distance (also compare with figure <ref>a). Then, coda waves from both events approximately travel along the same paths for all radiation directions and the event pairs are suited for a decorrelation analysis (method A).
This is the near field that is also excluded in inter-receiver interferometry.
A typical wavelength for the data set is 200m (frequency 18Hz, velocity 3.6km/s) and marked as a horizontal line in figure <ref>. Event pairs with a shorter distance are excluded from further analysis.





The second observation is, that considerably more data points are located on the right side than on the left side of the cross-correlation function. Because the waveform of the shallower event is used as the first signal in the cross-correlation calculated with equation <ref>,
maxima on the right side (positive side) of the cross-correlation function correspond to waves traveling from the shallower to the deeper event; maxima on the left side (negative side) vice versa.
Figure <ref> shows a stacked histogram of the number of data points on the right and left side of the cross-correlation function as a function of inclination angle between the two events. For events with a low inclination angle  (different depth, similar epicenter) maxima on the right side of the cross-correlation function are predominant indicating that waves going upwards are less likely to be recorded as coda than waves which leave the source in the downward direction.

A higher intrinsic attenuation in the top 10km could explain this observation, because coda waves that are initially going upwards spend more time in the uppermost crust and would suffer stronger attenuation than the initially downgoing waves.

For high inclination angles, events are located at similar depth and maxima on both sides of the cross-correlation are observed with similar frequency.



Thirdly, maxima of cross-correlation functions predominantly have a positive polarity. Maxima with negative polarity can be observed; often these have a larger time lag and arise after positive maxima of event pairs with similar distance (compare figure <ref>).
This observation could indicate that the correct inter-event travel time is associated with the positive maximum.
However, due to its oscillatory character the correlation function shows a down swing following the positive peak which might in some cases be of larger amplitude due to noise. 
To measure the wave velocity it is therefore desirable to search only for maxima with positive polarity.
Before, cases in which a truly negative peak can occur at the correct time have to be identified and excluded.

Different focal mechanisms of the two events in the pair can lead to a negative peak, because the polarity of radiated S waves into the inter-event direction might be different for both events (figure <ref>).
Simulations performed by <cit.> confirm this insight.

Figure <ref>a displays the focal mechanism of 13 larger earthquakes determined by <cit.>. The variability in these 13 mechanisms gives an idea of the variability of the focal mechanisms of all earthquakes in the swarm under the reasonable assumption that most smaller earthquakes have a similar focal mechanism.

The average of the observed polarity of the maxima in the cross-correlation function and the number of event pairs are displayed in figure <ref>b and <ref>c in binned polar plots as a function of inter-event inclination and azimuth. Note, that near the N-axis both positive and negative polarities can be observed. For other directions positive polarity is predominant. The distribution of number of event pairs reflects the shape of the 2018 earthquake swarm.

For the median focal mechanism with median slip, rake and dip the S wave radiation pattern according to chapter 4.3 of <cit.> is plotted alongside. No energy is radiated along the axis of largest compression and dilatation and neutral axis (P-, T-, N-axis). The maximal energy of shear waves is radiated along the fault plane and auxiliary plane. 

For exactly the same focal mechanism shear waves of the same polarity are emphasized by the cross-correlation function independent of the proximity of the two events and a peak of positive polarity will emerge in the cross-correlation function. For similar focal mechanisms a peak of positive polarity is still expected for all radiation directions except near the P-, T- and N-axis, where the shear wave polarity may be completely different for both events. Additionally, the amplitude of the radiation in inter-event direction is diminutive near the P-, T- and N-axis.

Therefore, depending on the exact orientation and focal mechanisms of the two earthquakes a peak of positive or negative polarity can arise in the cross-correlation function at the inter-event travel time or there may be no peak at all at this lag time <cit.>. 

In order to circumvent these error-prone situations and in order to exclude a possible negative polarity at the correct time, 
event pairs with azimuths larger than 320° and inclinations between 20° and 50° (region marked in figure <ref>), all near the N-axis of the median focal mechanism, are excluded from the further analysis.




Finally, figure <ref>a-e displays the lag times of positive maxima of the cross-correlations functions over event distance for the different event clusters a-e. Beside the previously mentioned constraints both events need to be from the same cluster. The velocity inside each cluster is determined by a robust mean <cit.> of the apparent velocities of each event pair. The weighting function for the iteratively reweighted least squares (IRLS) algorithm is the Hampel's 17A function with tuning constants a=1, b=2 and c=3. Panel f of figure <ref> displays the robustly determined shear wave velocities inside the five clusters together with its median absolute deviation (MAD) as error. The distribution of apparent seismic velocities in each event cluster is also indicated, together with its median value. The determined values for the shear wave velocities together with the corresponding MADs are given in panels a-e of figure <ref>. The velocities range from 3.49km/s to 4.16km/s with MADs between 0.25km/s and 0.88km/s. Cluster a has a significantly larger velocity of 4.16(51)km/s
and cluster d a slightly larger velocity of 3.85(88)km/s than the other clusters b, c and e with velocities between 3.49km/s and 3.72km/s. The variation around the mean velocities is high, especially for cluster d.




§ DISCUSSION




  §.§.§ Comparison between different methods of inter-source interferometry



























All of the three methods introduced in section <ref> aim at extracting the travel times between event locations. Method C was used in this study. Method B and C make use of event pairs which are separated by a larger distance than the dominant wavelength compared to method A for which event pairs must be separated by less than the dominant wavelength.
More event pairs fulfill this requirement for method B and C and therefore more data can be collected than for method A.
Method B which uses the direct wave train and method A do have the advantage that the waveforms between the event pair are expected to be more similar compared to method C. This makes outlying data points as in figure <ref> less likely.
A problem with method B is the constraint that stations have to be present in, but also be restricted to the stationary phase zone. Depending on the data set this can lead to a very low amount of stations or even no stations which can be used in the analysis; this was the case for the present data set.

The assumptions for method A are stronger than for method C. To calculate the event separation from the decorrelation value the same focal mechanism for the two earthquakes and the same scattering and attenuation strength in all directions is required <cit.>. The second assumption is not valid for the present application, otherwise a balanced distribution of maxima between the left and right side of the cross-correlation function would be expected (compare figure <ref>).

Method C in contrary assumes 1) a heterogeneous scatterer distribution in the stationary phase zone and 2)
similar focal mechanisms to the extent that the polarity of waves which are radiated into the inter-event direction is the same for both events. The same polarity can be guaranteed for similar focal mechanisms by excluding pairs with inter-event directions near the P-, N- or T-axis (compare figure <ref>). 
This step and the underlying constraint of similar focal mechanisms could also be omitted in favor of a more dedicated quality check of the calculated cross-correlations and their maxima or by taking into account the focal mechanisms of the earthquakes and their impact on the polarity of maxima in the inter-event cross correlation functions.

Processing for method A is sophisticated compared to the straight-forward processing of method C.
On the other hand, method A has the advantage over method C to not depend on the time lag of the maxima in the cross-correlation function. Therefore a possible error in origin times does not play a role as in method C.



  §.§.§ Shear wave velocities for different earthquake clusters


The inter-event interferometry was for the first time applied to map shear wave velocity in the source region of a fluid-driven earthquake swarm. 
The measured shear wave velocity v S=4.16(51)km/s of the first cluster a is significantly higher than the velocities of the clusters b, c and e of around 3.49 to 3.72km/s. The measured velocity for cluster d is around 3.85km/s and higher than the mean velocity of approximately 3.6km/s for the three clusters b, c and e. The data points of cluster d have a high median absolute deviation of 0.88km/s which could indicate that this cluster contains subclusters with different apparent velocities.
The higher velocities of clusters a and d are robust, insofar they are consistent if some of the applied constraints are relaxed or changed.

The method as applied in this study assumes that the used earthquake catalog is correct in origin times and locations.
An error in origin times shifts the determined inter-event travel time by the difference of the error in origin time of the first and second event (see equation <ref>).
For errors in origin locations the distance between the events is changed.
Both errors affect the velocity estimate and it is expected that the high variability of apparent velocities of different event pairs in one cluster is due to them. However, these errors do not introduce a systematic bias, but lead to a higher median absolute deviation.

For the future I suggest to extent the presented method to a full-fledged tomography by taking into account direct wave onsets, double difference travel times <cit.> and inter-event travel times and by inverting for origin times, locations and a seismic velocity grid in a single joint inversion.




  §.§.§ Interpretation of variance in shear wave velocity




Because there is neither a temporal nor spatial overlap between the defined earthquake clusters it is difficult to tell from the results alone if the differences in shear wave velocity are of temporal or spatial nature. Therefore, I will compare the results to <cit.> who observed a temporal variability and <cit.> who performed a classical travel time tomography with spatial information.

<cit.> observed P wave to S wave velocity ratios (s) lower than 1.45 at the beginning of the swarms in the years 1997, 2000 and 2008. Such low values of  have been observed earlier for the fault region of 1997 and 2008 swarms <cit.>.
<cit.> again observed lower velocity ratios down to 1.59 for some clusters of the 2014 earthquake swarm.
<cit.> use the Gassmann equations <cit.> which describe seismic velocities for porous fluid saturated rocks as a function of porosity. For these equations the shear wave velocity does not depend on porosity. <cit.> argue that the reduced velocity ratio is caused by a transition of the pore fluid to the gaseous phase and a resulting decrease in P wave velocity. 

A temporary decrease in P wave velocity could result in a systematic error of the underlying double-difference earthquake catalog at that period of time which in turn might lead to a higher observed shear wave velocity as in this study.
It is therefore important to calculate and examine the s of the five defined clusters of the 2018 swarm. In the appendix <ref> I show that the velocity ratio of the five clusters is approximately fixed at 1.68 (figure <ref>) which is a typical value for this region <cit.>. Please note that no assertion is made about a possible reduction of at the beginning of the individual clusters. The fixed indicates that the observed higher shear wave velocity of the first cluster a is not due to a temporal variability of /, but rather a spatial variation in shear wave velocity.

<cit.> speculate that a high attenuating zone ascertained at the same location as the velocity anomaly of this study, could be an indication of uprising fluids in permeable channels.

With a of 1.68 the observed S wave velocity translates into a P wave velocity of 7.0km/s for cluster a and P wave velocities down to approximately 6.0km/s for the swarm stages b, c and e. This value is consistent with the CEL09 seismic profile at depth of around 10km <cit.>. In figure <ref> the event clusters are plotted into horizontal and vertical slices of the shear wave model of <cit.>. 
As a matter of fact, the patch ruptured by cluster a with highest observed velocity is located in a model box with high shear wave velocity. Cluster d which also ruptured a region with increased velocity is located partly in the same box meaning that the higher shear wave velocity is consistent with previous observations.
The resolution of the model of <cit.> is bounded on the low end by its block size of 4km× 4km× 3km. The resolution in this study is about the extent of the earthquake clusters (≈ 1km) and already better than in a classical tomography of the same region. The resolution in a tomographic study based on inter-event interferometry depends on the inter-event distances (<100m). Therefore, the resolution is expected to be superior inside the swarm regions compared to classical tomography once the results of this study are used for a tomographic inversion.











§ CONCLUSIONS


A method of coda wave inter-source interferometry was applied to determine the shear wave velocity in different parts of the source volume of the 2018 West Bohemia earthquake swarm. Although, similar methods were used before for earthquake relocation procedures, this is the first study that allowed to map the seismic velocity. This is facilitated by a high quality earthquake catalog that has been relocated with the double difference technique. The work should be considered as a preparatory study for source region tomography and therefore several features visible in the inter-event cross-correlation functions were explained.
The resolution of seismic velocity of this study is already better than classical tomography studies of the same region.  I expect the resolution to improve further when using insights of this study within a joint tomographic approach.



  
Acknowledgments

I thank Martin Bachura for providing data from WEBNET seismic stations  and for providing the double difference earthquake catalog. Furthermore I thank Laura Barth and Thomas Plenefisch for the focal mechanisms.
All data used in this study can be downloaded at <https://doi.org/10.5281/zenodo.3741465> <cit.>.
Discussions with Christoph Sens-Schönfelder culminating in a full-blown review of the manuscript are very much appreciated.
I thank Jenny Borns and Ulrich Wegler for proof-reading the manuscript.
Comments from two anonymous reviewers and associate editor Nori Nakata considerably improved the manuscript.
Data processing and plotting was performed with ObsPy <cit.> and the Python packages NumPy, SciPy, matplotlib and statsmodels <cit.>.
This article can be reproduced with the source code provided at <https://github.com/trichter/inter_source_interferometry> <cit.>.






§ SOURCE REGION VELOCITY RATIOS




To check whether the observed variability of shear wave velocity involve changes in the , the method of <cit.> is applied to the data set of this study.
The method can be explained with the following equations






    S/R    = 1/Rv P/v S P    and
     P     =  P -  P  ,  S  =  S -  S  .


P and S are the travel time difference of P resp. S waves for an event pair with index i at station with index j.
Because all event pairs are inverted together, the means over all stations P resp. S have to be subtracted from data points of each event pair. The can be obtained by a regression between differential S wave travel times S and differential P wave travel times P. R is a scaling factor, which is initially set to 1.

For the present data, picks from all available earthquakes inside each event cluster (including earthquakes with a magnitude lower than 1.9) are used to calculate P and S. In the first step, not the mean, but the median is used to calculate preliminary differential travel times P, prel, S, prel (equation <ref>).

Erroneous data points with

    |S, prel - 1.7P, prel| > 0.1s


correspond to unrealistic s and are removed. The removed data points resulted from wrongly associated picks in the used catalog. P, prel, S, prel are recalculated by removing the median. In a next step, data points with


    √(P, prel^2 + S, prel/R^2) > 0.35s


are additionally removed. This is in accordance with <cit.> and removes event pairs whose locations are far away from each other. Here, data points outside a circle (or ellipse for R≠1) are removed to not affect the subsequent regression. Finally, P, S are calculated by removing the mean for each event pair (equation <ref>) and a robust orthogonal L1 regression is performed by testing different ratios with brute force (equation <ref>). For R=1 I obtain s of 1.73, 1.70, 1.72, 1.71, 1.71 for the earthquake clusters a, b, c, d and e.

<cit.> argue that errors in the differential S wave travel times due to different take-off angles are theoretically larger by a factor R=v P/v S than the corresponding errors for differential P wave travel times. They therefore suggest to determine the iteratively and scale the S wave differential travel times appropriately. With this approach using the same equations as above with R=v P/v S, velocity ratios of 1.68, 1.66, 1.69, 1.68, 1.69 are obtained for the clusters a, b, c, d and e (see figure <ref>). These estimates are slightly lower than the estimates for R=1, but in any case the is similar for all analyzed earthquake clusters and consistent with previous observations <cit.>.

`^̂M=5








49


[Aki and Richards(2002)]Aki2002
Aki K and Richards PG (2002), Quantitative Seismology, University Science
  Books, 2nd edition.

[Bachura and Fischer(2016a)]Bachura2016
Bachura M and Fischer T (2016a), Coda Attenuation Analysis in the
  West Bohemia/Vogtland Earthquake Swarm Area, Pure and Applied
  Geophysics, 173(2):425–437, 10.1007/s00024-015-1137-3.

[Bachura and Fischer(2016b)]Bachura2016a
Bachura M and Fischer T (2016b), Detailed velocity ratio mapping
  during the aftershock sequence as a tool to monitor the fluid activity within
  the fault plane, Earth and Planetary Science Letters, 453:215–222,
  10.1016/j.epsl.2016.08.017.

[Brenguier et al.(2008)Brenguier, Campillo, Hadziioannou, Shapiro,
  Nadeau, and Larose]Brenguier2008
Brenguier F, Campillo M, Hadziioannou C, Shapiro NM, Nadeau RM, and Larose E
  (2008), Postseismic relaxation along the San Andreas fault at Parkfield from
  continuous seismological observations., Science,
  321(5895):1478–1481, 10.1126/science.1160943.

[Curtis et al.(2009)Curtis, Nicolson, Halliday, Trampert, and
  Baptie]Curtis2009
Curtis A, Nicolson H, Halliday D, Trampert J, and Baptie B (2009), Virtual
  seismometers in the subsurface of the Earth from seismic interferometry,
  Nature Geoscience, 2(10):700–704, 10.1038/ngeo615.

[Dahm and Fischer(2013)]Dahm2013
Dahm T and Fischer T (2013), Velocity ratio variations in the source region of
  earthquake swarms in NW bohemia obtained from arrival time
  double-differences, Geophysical Journal International,
  196(2):957–970, 10.1093/gji/ggt410.

[Eulenfeld(2020)]Eulenfeld2020_sourcecode
Eulenfeld T (2020), Inter-source interferometry with cross-correlation of coda
  waves, 10.5281/zenodo.3936439.

[Eulenfeld and Wegler(2016)]Eulenfeld2016
Eulenfeld T and Wegler U (2016), Measurement of intrinsic and scattering
  attenuation of shear waves in two sedimentary basins and comparison to
  crystalline sites in Germany, Geophysical Journal International,
  205(2):744–757, 10.1093/gji/ggw035.

[Eulenfeld and Wegler(2017)]Eulenfeld2017
Eulenfeld T and Wegler U (2017), Crustal intrinsic and scattering attenuation
  of high-frequency shear waves in the contiguous United States, Journal
  of Geophysical Research: Solid Earth, 122(6):4676–4690,
  10.1002/2017JB014038.

[Fischer et al.(2019)Fischer, Bachura, Doubravová, and
  Horálek]Fischer2019AGU
Fischer T, Bachura M, Doubravová J, and Horálek J (2019), The 2018
  seismic activity of West Bohemia–from swarm to mainshock-aftershock
  sequence., in Geophysical Research Abstracts, volume 21.

[Fischer et al.(2014)Fischer, Horálek, Hrubcová,
  Vavryčuk, Bräuer, and Kämpf]Fischer2014
Fischer T, Horálek J, Hrubcová P, Vavryčuk V, Bräuer K,
  and Kämpf H (2014), Intra-continental earthquake swarms in
  West-Bohemia and Vogtland: A review, Tectonophysics, 611:1–27,
  10.1016/j.tecto.2013.11.001.

[Gaebler et al.(2015)Gaebler, Eulenfeld, and Wegler]Gaebler2015
Gaebler PJ, Eulenfeld T, and Wegler U (2015), Seismic scattering and
  absorption parameters in the W-Bohemia/Vogtland region from elastic and
  acoustic radiative transfer theory, Geophysical Journal
  International, 203(3):1471–1481, 10.1093/gji/ggv393.

[Hillers et al.(2015)Hillers, Retailleau, Campillo, Inbal, Ampuero,
  and Nishimura]Hillers2015
Hillers G, Retailleau L, Campillo M, Inbal A, Ampuero JP, and Nishimura T
  (2015), In situ observations of velocity changes in response to tidal
  deformation from analysis of the high-frequency ambient wavefield,
  Journal of Geophysical Research: Solid Earth, 120(1):210–225,
  10.1002/2014JB011318.

[Hobiger et al.(2016)Hobiger, Wegler, Shiomi, and
  Nakahara]Hobiger2016
Hobiger M, Wegler U, Shiomi K, and Nakahara H (2016), Coseismic and
  post-seismic velocity changes detected by passive image interferometry:
  Comparison of one great and five strong earthquakes in Japan,
  Geophysical Journal International, 205(2):1053–1073,
  10.1093/gji/ggw066.

[Hong and Menke(2006)]Hong2006
Hong TK and Menke W (2006), Tomographic investigation of the wear along the
  San Jacinto fault, southern California, Physics of the Earth and
  Planetary Interiors, 155(3-4):236–248, 10.1016/j.pepi.2005.12.005.

[Hrubcová et al.(2005)Hrubcová, Środa,
  Špičák, Guterch, Grad, Keller, Brueckl, and
  Thybo]Hrubcova2005
Hrubcová P, Środa P, Špičák A, Guterch A, Grad M,
  Keller GR, Brueckl E, and Thybo H (2005), Crustal and uppermost mantle
  structure of the Bohemian Massif based on CELEBRATION 2000 data,
  Journal of Geophysical Research: Solid Earth, 110(11):1–21,
  10.1029/2004JB003080.

[Huber(2014)]Huber2014
Huber PJ (2014), Robust Statistics, in M Lovric, editor, International
  Encyclopedia of Statistical Science, pages 1248–1251, Springer Berlin
  Heidelberg, Berlin, Heidelberg, 10.1007/978-3-642-04898-2_594.

[Hunter(2007)]Hunter2007
Hunter JD (2007), Matplotlib: A 2D Graphics Environment, Computing in
  Science & Engineering, 9(3):90–95, 10.1109/MCSE.2007.55.

[Institute of Geophysics, Academy of Sciences of the Czech
  Republic(1991)]WebnetData
Institute of Geophysics, Academy of Sciences of the Czech Republic (1991),
  West Bohemia Local Seismic Network, 10.7914/SN/WB.

[Lin and Shearer(2007)]Lin2007
Lin G and Shearer P (2007), Estimating local Vp/Vs ratios within similar
  earthquake clusters, Bulletin of the Seismological Society of
  America, 97(2):379–388, 10.1785/0120060115.

[Málek et al.(2005)Málek, Horálek, and
  Janský]Malek2005
Málek J, Horálek J, and Janský J (2005), One-dimensional
  qP-wave velocity model of the upper crust for the West Bohemia/Vogtland
  Earthquake swarm region, Studia Geophysica et Geodaetica,
  49(4):501–524, 10.1007/s11200-005-0024-2.

[Mavko et al.(2009)Mavko, Mukerji, and Dvorkin]Mavko2009
Mavko G, Mukerji T, and Dvorkin J (2009), The Rock Physics Handbook,
  Cambridge University Press, Cambridge, second edi edition,
  10.1017/CBO9780511626753.

[Megies et al.(2011)Megies, Beyreuther, Barsch, Krischer, and
  Wassermann]Megies2011
Megies T, Beyreuther M, Barsch R, Krischer L, and Wassermann J (2011), ObsPy -
  what can it do for data centers and observatories?, Annals of
  Geophysics, 54(1):47–58, 10.4401/ag-4838.

[Mousavi et al.(2015)Mousavi, Bauer, Korn, and Hejrani]Mousavi2015
Mousavi S, Bauer K, Korn M, and Hejrani B (2015), Seismic tomography reveals a
  mid-crustal intrusive body, fluid pathways and their relation to the
  earthquake swarms in West Bohemia/Vogtland, Geophysical Journal
  International, 203(2):1113–1127, 10.1093/gji/ggv338.

[Mousavi et al.(2017)Mousavi, Haberland, Bauer, Hejrani, and
  Korn]Mousavi2017
Mousavi S, Haberland C, Bauer K, Hejrani B, and Korn M (2017), Attenuation
  tomography in West Bohemia/Vogtland, Tectonophysics, 695:64–75,
  10.1016/j.tecto.2016.12.010.

[Oliphant(2006)]NumPy
Oliphant TE (2006), A guide to NumPy, USA, Trelgol Publishing.

[Plenefisch and Barth(2019)]Plenefisch2019EGU
Plenefisch T and Barth L (2019), The May 2018 earthquake swarm in
  Vogtland/NW-Bohemia: Spatiotemporal evolution and focal mechanism
  determinations, in Geophysical Research Abstracts, volume 21, pages
  EGU2019–9356.

[Richter et al.(2014)Richter, Sens-Schönfelder, Kind, and
  Asch]Richter2014
Richter T, Sens-Schönfelder C, Kind R, and Asch G (2014), Comprehensive
  observation and modeling of earthquake and temperature-related seismic
  velocity changes in northern Chile with passive image interferometry,
  Journal of Geophysical Research: Solid Earth, 119(6):4747–4765,
  10.1002/2013JB010695.

[Robinson et al.(2011)Robinson, Sambridge, and Snieder]Robinson2011
Robinson DJ, Sambridge M, and Snieder R (2011), A probabilistic approach for
  estimating the separation between a pair of earthquakes directly from their
  coda waves, Journal of Geophysical Research: Solid Earth,
  116(4):1–17, 10.1029/2010JB007745.

[Robinson et al.(2013)Robinson, Sambridge, Snieder, and
  Hauser]Robinson2013
Robinson DJ, Sambridge M, Snieder R, and Hauser J (2013), Relocating a cluster
  of earthquakes using a single seismic station, Bulletin of the
  Seismological Society of America, 103(6):3057–3072,
  10.1785/0120130004.

[Sato et al.(2012)Sato, Fehler, and Maeda]Sato2012
Sato H, Fehler MC, and Maeda T (2012), Seismic Wave Propagation and Scattering
  in the Heterogeneous Earth, Springer, second edition,
  10.1007/978-3-642-23029-5.

[Schimmel and Paulssen(1997)]Schimmel1997
Schimmel M and Paulssen H (1997), Noise reduction and detection of weak,
  coherent signals through phase-weighted stacks, Geophysical Journal
  International, 130(2):497–505, 10.1111/j.1365-246X.1997.tb05664.x.

[Seabold and Perktold(2010)]statsmodels
Seabold S and Perktold J (2010), Statsmodels: Econometric and Statistical
  Modeling with Python, in Proceedings of the 9th Python in Science
  Conference, volume 57, page 61.

[Sens-Schönfelder and Eulenfeld(2019)]Sens-Schonfelder2019
Sens-Schönfelder C and Eulenfeld T (2019), Probing the in situ elastic
  nonlinearity of rocks with earth tides and seismic noise, Physical
  Review Letters, 122(13):138501, 10.1103/PhysRevLett.122.138501.

[Sens-Schönfelder and Wegler(2006)]Sens-Schonfelder2006
Sens-Schönfelder C and Wegler U (2006), Passive image interferemetry and
  seasonal variations of seismic velocities at Merapi Volcano, Indonesia,
  Geophysical Research Letters, 33(21):L21302,
  10.1029/2006GL027797.

[Shapiro and Campillo(2004)]Shapiro2004
Shapiro NM and Campillo M (2004), Emergence of broadband Rayleigh waves from
  correlations of the ambient seismic noise, Geophysical Research
  Letters, 31(7):L07614, 10.1029/2004GL019491.

[Snieder(2004)]Snieder2004
Snieder R (2004), Extracting the Green's function from the correlation of coda
  waves: a derivation based on stationary phase., Physical review. E,
  Statistical, nonlinear, and soft matter physics, 69(4):046610,
  10.1103/PhysRevE.69.046610.

[Snieder and Vrijlandt(2005)]Snieder2005
Snieder R and Vrijlandt M (2005), Constraining the source separation with coda
  wave interferometry: Theory and application to earthquake doublets in the
  Hayward fault, California, Journal of Geophysical Research: Solid
  Earth, 110(4):1–15, 10.1029/2004JB003317.

[Sun et al.(2016)Sun, Zhang, and Wen]Sun2016
Sun L, Zhang M, and Wen L (2016), A new method for high-resolution event
  relocation and application to the aftershocks of Lushan Earthquake, China,
  Journal of Geophysical Research: Solid Earth, 121:2539–2559,
  10.1002/2016JB012840.

[Tonegawa and Nishida(2010)]Tonegawa2010
Tonegawa T and Nishida K (2010), Inter-source body wave propagations derived
  from seismic interferometry, Geophysical Journal International,
  183(2):861–868, 10.1111/j.1365-246X.2010.04753.x.

[Vavryčuk(2011)]Vavrycuk2011
Vavryčuk V (2011), Tensile earthquakes: Theory, modeling, and
  inversion, Journal of Geophysical Research: Solid Earth,
  116(12):1–14, 10.1029/2011JB008770.

[Virtanen et al.(2020)Virtanen, Gommers, Oliphant, Haberland, Reddy,
  Cournapeau, Burovski, Peterson, Weckesser, Bright, van der Walt, Brett,
  Wilson, Millman, Mayorov, Nelson, Jones, Kern, Larson, Carey, Polat, Feng,
  Moore, VanderPlas, Laxalde, Perktold, Cimrman, Henriksen, Quintero, Harris,
  Archibald, Ribeiro, Pedregosa, and van Mulbregt]SciPy
Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D,
  Burovski E, Peterson P, Weckesser W, Bright J, van der Walt SJ, Brett M,
  Wilson J, Millman KJ, Mayorov N, Nelson ARJ, Jones E, Kern R, Larson E, Carey
  CJ, Polat İ, Feng Y, Moore EW, VanderPlas J, Laxalde D, Perktold J, Cimrman
  R, Henriksen I, Quintero EA, Harris CR, Archibald AM, Ribeiro AH, Pedregosa
  F, and van Mulbregt P (2020), SciPy 1.0: fundamental algorithms for
  scientific computing in Python, Nature Methods,
  10.1038/s41592-019-0686-2.

[Weaver and Lobkis(2002)]Weaver2002
Weaver R and Lobkis O (2002), On the emergence of the Green's function in the
  correlations of a diffuse field: Pulse-echo using thermal phonons,
  Ultrasonics, 40(1-8):435–439, 10.1016/S0041-624X(02)00156-7.

[WEBNET group et al.(2020)WEBNET group, Bachura, and Barth]dataset_swarm2018
WEBNET group, Bachura M, and Barth L (2020), Seismological dataset for 2018 West
  Bohemia earthquake swarm, 10.5281/zenodo.3741465.

[Wegler and Sens-Schönfelder(2007)]Wegler2007
Wegler U and Sens-Schönfelder C (2007), Fault zone monitoring with
  passive image interferometry, Geophysical Journal International,
  168(3):1029–1033, 10.1111/j.1365-246X.2006.03284.x.

[Zang et al.(2014)Zang, Oye, Jousset, Deichmann, Gritto, McGarr,
  Majer, Bruhn, McGarr, Majer, and Bruhn]Zang2014
Zang A, Oye V, Jousset P, Deichmann N, Gritto R, McGarr A, Majer E, Bruhn D,
  McGarr A, Majer E, and Bruhn D (2014), Analysis of induced seismicity in
  geothermal reservoirs - An overview, Geothermics, 52:6–21,
  10.1016/j.geothermics.2014.06.005.

[Zhang and Thurber(2003)]Zhang2003
Zhang H and Thurber CH (2003), Double-difference tomography: The method and
  its application to the Hayward Fault, California, Bulletin of the
  Seismological Society of America, 93(5):1875–1889,
  10.1785/0120020190.

[Zhao and Curtis(2019)]Zhao2019
Zhao Y and Curtis A (2019), Relative source location using coda-wave
  interferometry: Method, code package, and application to mining-induced
  earthquakes, Geophysics, 84(3):F73–F84,
  10.1190/geo2018-0601.1.

[Zhao et al.(2017)Zhao, Curtis, and Baptie]Zhao2016
Zhao Y, Curtis A, and Baptie B (2017), Locating microseismic sources with a
  single seismometer channel using coda wave interferometry,
  GEOPHYSICS, 82(3):A19–A24, 10.1190/geo2016-0404.1.



