Towards source region tomography with inter-source interferometry: Shear wave velocity from 2018 West Bohemia swarm earthquakes

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.5 km/s and 4.2 km/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 distribution of scatterers on the polarity and location of the maxima in the cross-correlation functions is discussed.

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.5 km/s and 4.2 km/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 distribution of scatterers on the polarity and location of the maxima in the cross-correlation functions is discussed.
Key points: • Travel times between earthquake locations from cross-correlations of coda waves • Clear move-out with shear wave velocity • Spatial variability of seismic velocity at region of 2018 earthquake swarm

Introduction
The trending 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 (Weaver and Lobkis, 2002;Snieder, 2004;Shapiro and Campillo, 2004). In this inter-receiver setting, the Green's function or the crosscorrelation function for arbitrary time-invariant noise sources can be used to monitor the wave velocity between and near the receivers (e.g. Sens-Schönfelder and Wegler, 2006;Wegler and Sens-Schönfelder, 2007;Brenguier et al., 2008;Richter et al., 2014;Hillers et al., 2015;Hobiger et al., 2016;Sens-Schönfelder and Eulenfeld, 2019). 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 branch of studies used coda cross-correlations of nearby earthquakes to determine the inter-event travel time from the decorrelation of the coda wave field of two earthquakes (referred to as method A, Snieder and Vrijlandt, 2005;Robinson et al., 2011Robinson et al., , 2013Zhao et al., 2017;Zhao and Curtis, 2019). 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. Snieder and Vrijlandt (2005) 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. Hong and Menke (2006) used spatial reciprocity (Aki and Richards, 2002, p. 28) 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. Curtis et al. (2009) further demonstrated the principle of inter-source interferometry by comparing real and virtual recordings of the Sichuan earthquake. Different from Hong and Menke (2006) they used the full seismic recording to mainly reconstruct surface waves and therefore can restrict the used stations to the stationary phase zone. Curtis et al. (2009) 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. Tonegawa and Nishida (2010) selected recordings in the stationary phase zone of body waves propagation between deep earthquakes and extended the inter-source interferometry concept to P and S body waves. Inter-source interferometry as in Curtis et al. (2009) and Tonegawa and Nishida (2010) 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. Sun et al. (2016) 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). Figure 1 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 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 degasing, Quaternary volcanoes and earthquake swarms (Fischer et al., 2014). The topography of West Bohemia and Vogtland is shown in figure 2 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 might allow for a spatial resolution that is superior to conventional seismic tomography. Lin and Shearer (2007) developed a double difference arrival method to determine the ratio of P wave velocity to S wave velocity (v P /v S ratio) for a region embraced by a cluster of earthquake. Dahm and Fischer (2013) applied a similar method to the West Bohemia earthquake swarms of the years 1997, 2000 and 2008 and found a time-dependent v P /v S ratio between 1.38 and 1.70. Bachura and Fischer (2016b) showed that the v P /v S ratio 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 differently from previous studies of inter-event interferometry which focus on the relocation of earthquakes.

Data and method
The double difference catalog of the 2018 swarm used in this study was compiled by Martin Bachura (Fischer et al., 2019). It consists of approximately 1000 earthquakes with local magnitudes larger than 1.3. The uncertainties in origin locations are 50 m. For the purpose of the present study, the catalog is divided into 5 clusters named a-e with earthquakes larger than magnitude 1.8 separated by the times 2018- 05-10, 2018-05-11 03:10, 2018-05-12 06:00, 2018-05-21, 2018-05-25, 2018-06-19 (all   . The ray paths in blue give rise to the peak in the cross-correlation function -dark blue parts are the same for both events, light blue parts are different for both events and determine the travel time difference. Light gray ray paths are not emphasized by the crosscorrelation, because their arrival times do not fall in the selected time window (method A, B, C) or because the corresponding recorded waveforms are not similar (method B, C). In method B, light orange paths give rise to the peak in the cross-correlation function but are usually not considered because the corresponding station is not located in the stationary phase zone. On the bottom, important aspects of each method are listed.   with the box in figure 2 is excluded from the analysis. Cluster a started the 2018 swarm activity at a depth of 9 km to 10 km. 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 7 km). After June 19 the activity faded out with cluster e. 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 (Plenefisch and Barth, 2019).
The waveforms of the earthquakes were registered on the Czech WEBNET stations with 250 Hz sampling rate. For this study, data from 9 WEBNET stations are used (compare figure 2). Waveforms are bandpass filtered in the frequency range 10 Hz to 40 Hz. Coda time windows are selected from 1 s after the S pick to 50 s 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). The following processing steps are applied to pairs of earthquakes. Pairs with an inter-event distance larger than 1 km are excluded due to the low fraction of scattered waves which travel the direct path between the events. Data of two different earthquakes registered 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 10 s are discarded. The first signal s 1 from the shallower event and second signal s 2 from the deeper event are trimmed to (t 1 , t 2 ) and the cross-correlation C(t) is calculated for each station and component with 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 (Schimmel and Paulssen, 1997) and in the process each trace is weighted by the length of the coda window used for the cross-correlation. Figure 4 displays examples of crosscorrelations for two different event pairs. In the left panel 4a 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 4b the peaks corresponding to the inter-event travel time only emerge after the stacking procedure. The relationship between the lag time of the maxima t max resulting from waves traveling between the earthquake locations and the inter-event travel time t travel  A clear move-out with a velocity of around 3.6 km/s is visible for the coda window. For the direct S wave window, the apparent travel time associated with random peaks in the cross-correlation function is lower than expected inter-event travel time. is sgn is the sign function and ∆t o,1 and ∆t o,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 two peaks t max,1 and t max,2 corresponding to waves traveling between event locations in both directions are visible in the cross-correlation function as in figure 4b, the inter-event travel time and difference in origin time errors can be determined directly with 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: t travel ≈ |t max |. In figure 5a 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 moveout with a velocity of around 3.6 km/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 (0 s, 2 s) 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 5b. 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 6 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.6 km/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 6). The observed move-out shows that most of the peaks arise due to the extraction of waves traveling between the earthquake locations. The maxima therefore reflect the inter-event travel times. 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 5a). 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 200 m (frequency 18 Hz, velocity 3.6 km/s) and marked as a horizontal line in figure 6. 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 1, 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 7 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, which are in vertical proximity of each other (different depth, similar epicenter), maxima on The corresponding S wave radiation pattern (arrows, Aki and Richards, 2002, chapter 4.3) has zero points at the P-, T-and N-axis. The largest amplitudes are expected at fault plane and auxiliary plane. The median focal mechanism with its P-and T-axis is also displayed in the two other panels together with an angular histogram of positive versus negative polarity of maxima of cross-correlation functions (b) and with an angular histogram of the number of event pairs (c). Negative polarity occurs preferentially for specific orientations between the two earthquakes and is near the N-axis. Therefore, event pairs with azimuth and inclination inside the region marked by the gray line are later excluded from analysis. All polar plots use a lower-hemisphere stereographic projection with labels of inclination and azimuth values.
In the middle panel only bins with more than 5 event pairs are displayed.
the right side of the cross-correlation function are predominant.
With the help of the Qopen software package Wegler, 2016, 2017), a scattering transport mean free path of around 100 km is estimated for a frequency range between 10 Hz to 20 Hz for the present data set. Similar values were obtained by Bachura and Fischer (2016a) for this region with a different data set. The scattering transport mean free path is defined by the average length in which the wave forgets its initial direction due to single or multiple scattering. The preference for maxima on the right side for low inclination angles seems to indicate that waves going upwards are less likely to be recorded as coda than waves which leave the source in the downward direction and are scattered or reflected in the lower crust. The low ratio between event depth of up to 10 km and the scattering transport mean free path of around 100 km supports this interpretation. For high inclination angles, events are in horizontal proximity 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 6). 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 8). Simulations performed by Sun et al. (2016) confirm this insight. Figure 8a displays the focal mechanism of 13 larger earthquakes determined by Plenefisch and Barth (2019). 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 observed polarity of the maxima in the cross-correlation function and the number of event pairs are displayed in figure 8b and 8c in polar histograms. 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 Aki and Richards (2002) 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-, Tand N-axis, where the shear wave polarity may be completely different for both events. 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 interevent travel time or there may be no peak at all at this lag time (Sun et al., 2016). Therefore, event pairs with azimuths larger than 320°and inclinations between 20°and 50°(region marked in figure 8), all near the N-axis of the median focal mechanism, are excluded from the further analysis.
Finally, figure 9a-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 (Eulenfeld and Wegler, 2016;Huber, 2014)  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 9. The velocities range from 3.5 km/s to 4.2 km/s with MADs between 0.3 km/s and 0.9 km/s. Cluster a has a significantly larger velocity of (4.2 ± 0.5) km/s and cluster d a slightly larger velocity of (3.9 ± 0.9) km/s than the other clusters b, c and e with velocities between 3.5 km/s and 3.7 km/s. The variation around the mean velocities is high, especially for cluster d.

Discussion
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.2 ± 0.5) km/s of the first cluster a is significantly higher than the velocities of the clusters b, c and e of around 3.5 to 3.7 km/s. The measured velocity for cluster d is around 3.85 km/s and higher than the mean velocity of approximately 3.6 km/s for the three clusters b, c and e. The data points of cluster d have a high median absolute deviation of 0.85 km/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 location. 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 2). 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 and inter-event travel times and by inverting for origin times, locations and a seismic velocity grid in a single joint inversion.

Comparison between different methods of inter-source interferometry
All of the three methods introduced in section 1 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 6 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 strength in all directions is required (Snieder and Vrijlandt, 2005). 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 7). Method C only assumes 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. This can be guaranteed for similar focal mechanisms by excluding pairs with inter-event directions near the P-, N-or T-axis (compare figure 8). 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.

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 Dahm and Fischer (2013) who observed a temporal variability and Mousavi et al. (2015) who performed a classical travel time tomography with spatial information. Dahm and Fischer (2013) observed P wave to S wave velocity ratios lower than 1.45 at the beginning of the swarms in the years 1997, 2000 and 2008. Such low values of v P /v S ratio have been observed earlier for the fault region of 1997 and 2008 swarms (Vavryčuk, 2011). Dahm and Fischer (2013) use the Gassmann equations (Mavko et al., 2009, chappter 6.3) 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. Dahm and Fischer (2013) 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. Bachura and Fischer (2016b) explained lower velocity ratios down to 1.59 for some clusters of the 2014 earthquake swarm with a similar argument. In appendix A, I calculate v P /v S ratios for the five defined clusters of the 2018 swarm. A reduced v P /v S ratio involving changes in the fluid content in the material could explain a change of velocity on the time scale of days. However, the velocity ratio is approximately fixed at 1.68 (figure 11) which is a typical value for this region (e.g. Málek et al., 2005); a reduced v P /v S ratio is not detected. Therefore the observed higher shear wave velocity of the first cluster a is not  due to a temporal variability, but rather a spatial pattern. With a v P /v S ratio of 1.68 the observed S wave velocity translates into a P wave velocity of 7.0 km/s for cluster a and P wave velocities down to approximately 6.0 km/s for the swarm stages b, c and e. This value is consistent with the CEL09 seismic profile at depth of around 10 km (Hrubcová et al., 2005). In figure 10 the event clusters are plotted into horizontal and vertical slices of the shear wave model of Mousavi et al. (2015). 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 obtained by the cluster approach in this study is already higher than in a classical tomography as in Mousavi et al. (2015). The resolution is expected to be superior inside the swarm regions 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 intra-source tomography and therefore several features visible in the intra-event crosscorrelation functions were explained. The resolution of seismic velocity of this study is already superior to 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. 1.5 1.6 1.7 1.8 1.9 velocity ratio v P /v S h) 1.5 1.6 1.7 1.8 1.9 i) 1.5 1.6 1.7 1.8 1.9 j) Figure 11: a-e) Velocity ratios v P /v S for each earthquake cluster a-e from double difference travel times. Displayed is the orthogonal distance regression with L1 norm (straight line) and velocity ratios of 1.5 and 1.9 (dotted lines). The number of used events and the determined velocity ratio is displayed in each panel. f-j) For each cluster a-e, the mean absolute error (orthogonal distance) is displayed for tested velocity ratios.

A. Source region velocity ratios
To check weather the observed variability of shear wave velocity involve changes in the v P /v S ratio, the method of Lin and Shearer (2007) is applied to the data set of this study. The method can be explained with the following equationŝ δt ij P and δt ij 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 <δt ij P > j resp. <δt ij S > j have to be subtracted from data points of each event pair. The v P /v S ratio can be obtained by a regression between differential S wave travel timesδt ij S and differential P wave travel timesδt ij 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 δt ij P and δt ij S . In the first step, not the mean, but the median is used to calculate preliminary differential travel timesδt ij P, prel ,δt ij S, prel (equation 6). Erroneous data points with δ t ij S, prel − 1.7δt ij P, prel > 0.1 s are located far away from the Wadati line and are removed. The removed data points resulted from wrongly associated picks in the used catalog.δt ij P, prel ,δt ij S, prel are recalculated by removing the median. In a next step, data points with δ t ij P, prel 2 + δ t ij S, prel /R 2 > 0.35 s are additionally removed. This is in accordance with Dahm and Fischer (2013) 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,δt ij P ,δt ij S are calculated by removing the mean for each event pair (equation 6) and a robust orthogonal L1 regression is performed by testing different ratios with brute force (equation 5). For R=1 I obtain v P /v S ratios of 1.73, 1.70, 1.72, 1.71, 1.71 for the earthquake clusters a, b, c, d and e. Lin and Shearer (2007) 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 v P /v S ratio 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 11). These estimates are slightly lower than the estimates for R=1, but in any case the v P /v S ratio is similar for all analyzed earthquake clusters and consistent with previous observations (e.g. Málek et al., 2005).