Event couple spectral ratio Q method for earthquake clusters: application to northwest Bohemia

We develop an amplitude spectral ratio method for event couples from clustered earthquakes to estimate seismic wave attenuation (Q−1) in the source volume. The method allows to study attenuation within the source region of earthquake swarms or aftershocks at depth, independent of wave path and attenuation between source region and surface station. We exploit the high-frequency slope of phase spectra using multitaper spectral estimates. The method is tested using simulated full wave-field seismograms affected by recorded noise and finite source rupture. The synthetic tests verify the approach and show that solutions are independent of focal mechanisms but also show that seismic noise may broaden the scatter of results. We apply the event couple spectral ratio method to northwest Bohemia, Czech Republic, a region characterized by the persistent occurrence of earthquake swarms in a confined source region at mid-crustal depth. Our method indicates a strong anomaly of high attenuation in the source region of the swarm with an averaged attenuation factor of Qp < 100. The application to S phases fails due to scattered P-phase energy interfering with S phases. The Qp anomaly supports the common hypothesis of highly fractured and fluid saturated rocks in the source region of the swarms in northwest Bohemia. However, high temperatures in a small volume around the swarms cannot be excluded to explain our observations.

estimated for the integral ray path between sources and stations (except for the work by Wcisło et al. (2018)) or secondly, because they focused on low frequencies, or both. E.g. Mousavi et al. (2017) used frequencies between 1 to 30 Hz and Gaebler et al. (2015) up to 32 Hz.
In this study, we aim to increase the spatial resolution and to resolve Q for waves traveling only within the small source region of the earthquake swarms. The developed event couple spectral ratio method is based on the assumption of an exponentially 25 decreasing spectral slope at high frequencies ω above the corner frequency of the earthquake, often referred to as the ω 2 model (Aki, 1980). From the ratio of the spectral slopes of two events one can estimate the attenuation of P and S phases for the ray path between the two events, given the differential travel time of both events. Matsumoto et al. (2009) exploit amplitude spectral ratios of direct P phases and normalize the spectra with the coda energy to compensate for source effects. Opposed to their approach we focus on the higher frequency content to achieve a higher resolution needed to map the compact source 30 volume.
The spatially compact seismic clusters in NW Bohemia provide us with a favorable case study scenario due to the high similarity of source characteristics (Michálek and Fischer, 2013). We test our method on data recorded from October 6 until October 13, 2008 and a double-difference relocated event catalog of 3841 events with local magnitudes between -0.9 and 3.5 (Fischer and Michálek, 2008). The high density of events during earthquake swarms clustering within a small and confined 35 2 region allows to infer the local attenuation from event couples, by applying the spectral ratio method (Aki, 1980) to their high frequency amplitude spectra as we will explain in the following section. One major issue when calculating spectral content of very few data samples is spectral leakage, as a result of the finiteness of the time window under study. In order overcome this problem Thomson (1982) proposed the multitaper method, which we employ using mtspec (Prieto et al., 2009;Krischer, 2016). 5

Method
A velocity spectrum A(ω) of a direct body wave phase originating from a source j recorded at a station i can be described as (Sanders, 1993): where ω is the angular frequency. S j (ω) describes the source spectrum and I i (ω) the instrument response. R i (ω) is the receiver 10 site effect. G i,j is the frequency independent geometric loss. The exponential term depends on the angular frequency ω and t * , the path integrated attenuation from the source to the receiver: with Q as the dimensionless quality factor, velocity v of the medium and ds a segment along the ray path from the source to the receiver. Attenuation is considered here as a combination of intrinsic and scattering losses. Instead of estimating a 15 total t * describing the full ray path's attenuation we estimate a local t * from velocity spectra of two earthquakes sharing the greater part of their ray paths from the seismogenic zone to a receiver. Site effects as well as the receivers response functions cancel out when two spectra recorded at the same site are analyzed by means of amplitude ratios. Let A j,0 and A j,1 be two velocity amplitude spectra of events E 0 and E 1 (in the following referred to as first and second event of a couple) recorded at a station j ( Figure 1). We assume that their source spectra S 0 (ω) and S 1 (ω) resemble each other to a degree where the effect 20 of random perturbations at high frequencies average out when the proposed method is applied to many couples. Taking the natural logarithm of the spectral ratio of A i,0 and A i,1 yields:

25
This equation describes a linear relation with frequency independent geometrical losses to the left of the negative sign in equation 3. The slope k of a line fitted to equation 3 can be used to derive the attenuation time t * in between the two sources from which Q −1 can easily be inferred using equation 2.
The far field amplitude spectrum A(ω) of P and S phases can be parameterized as follows: a seismic moment dependent low frequency plateau, the corner frequency f c and the high-frequency spectral decay approximately proportional to ω 2 resulting the station at which the displayed ray segment arrived as well as the underlying source of the segments. Grey points show hypocenters which occurred during the investigated time interval but were not attributed to an event couple. from finiteness of particle rise time and the rupture duration (Aki and Richards, 2002). To remove the dependence on seismic moment we investigate the high-frequency spectral decay, only. Furthermore, exploiting spectral content above the corner frequency also reduces source directivity effects on attenuation estimates (Cormier, 1982).
We infer the corner frequency based on previous studies on NW Bohemia seismic swarms. We use a relation proposed by Michálek and Fischer (2013) to calculate source radii r based on the moment M 0 of an event: where M 0 = 1.38M L + 10.3. The resulting source radii r are then converted to f c using with k = 0.32 (Madariaga, 1976), and β = 3.5 km/s, which is a reasonable assumption for the source region (Michálek and Fischer, 2013). We increase the lower frequency limit (f min ) used for our spectral analysis by additional 5 Hz with respect 10 4 to f c to account for uncertainties in M 0 and to ensure linearity of the high frequency decay. This approach allows frequency bands being wide enough for employing a stable linear regression. The upper frequency limit f max was chosen dependent on the Fresnel volume (see below) of a couple's first event ( Figure 1) and the upper corner frequency of the anti-alias filter of the recording equipment which is approximately 85 Hz. We calculate the power spectral density using the multitaper spectral analysis method (MTM) (Thomson, 1982;Park et al., 1987). With this method the time series is multiplied with several 5 orthogonal slepian tapers which are resistant to spectral leakage. The power spectral density is then reconstructed after Fourier transformation of the tapered samples and a weighted summation of the resulting spectra. A more exhaustive explanation can be found in Park et al. (1987). The applied code is a Python wrapper to the Fortran routine MTSpec (Prieto et al., 2009;Krischer, 2016). MTM achieves stable spectral estimates also for very short time windows. A critical parameter of the MTM is the number of slepian tapers as it balances the smoothness and precision of spectral estimates. We use the implemented default, 10 which is where bw is the bandwidth_factor which we set to bw = 4. Lower values prove to increase the number of outliers due to increased spectral leakage. Higher values did not change the results significantly but are more expensive to compute.
We impose strong geometrical constraints to select event couples with respect to a station as sketched in Figure 1. Ray tracing 15 is done based on a 1D velocity model suggested by Alexandrakis et al. (2014) for the seismogenic region combined with a regional crustal model proposed by Málek et al. (2000) (Fig. 4, left panel). The first geometrical constraint is the traversing distance (D t , green dashed line, Fig 1) between an event E 0 with respect to perpendicular projections of other hypocenters onto that path. The second constraint is the passing distance (D p , red dashed line in Fig 1) of that projection of E 1 onto the ray between E 0 and the station. We defined a minimum traversing distance of D t ≥ 1500 m to ensure that the signal of the second 20 event is attenuated sufficiently to be detectable in the couple's spectral ratio.
Subsequent to geometrical preselection upper frequency limits of the analyzed bandwidth are potentially corrected to lower values dependent on the 2 nd Fresnel volume in between event E 0 and a station. The wavelength from which this frequency limit can be deduced is given as (Matsumoto et al., 2009): 25 D p and D t are passing distance and the traversing distance as defined above. x is the distance from the passing point with respect to the second event E 1 to the receiver ( Figure 1). n is the number of the Fresnel volume. The frequency limit can then be deduced from the wavelength and the wave velocity of the underlying medium. This approach provides a physically meaningful limitation to impose on the frequency bandwidth. It ensures that E 1 is located within the Fresnel volume (grey shaded area in Figure 1) for the entire analyzed frequency band.

30
With this approach we get an estimate for the attenuation along the traversing distance (green dashed line in Figure 1) and when repeated for a large number of event couples can retrieve a median attenuation for the entire source region. The described method is advantageous over other methods which require handcrafted features like onset duration picking as it can be automatized given that an onset catalog is at hand.

Synthetic Study
In order to evaluate the expected number of exploitable couples given the geometrical constraints we calculate the relative number of pairs at discrete surface points covering the region of North-West Bohemia. The size of the blue points in Figure   3 represents the relative number of arriving shared ray paths from event pairs at possible station sites based on ray tracing through a 1D layered model (Fig. 4). Largest numbers of pairs are expected along a North-South striking patch which follows  Figure 2 shows the rays which penetrate the source volume and fulfill the geometrical requirements 10 described above. It shows that for events recorded at the most significant stations NKC and LBC the highest ray density and therefore sensitivity is in the lower half of the seismogenic zone. This bias is more pronounced for recordings at station LBC.
Also, these ray segments sample the volume up to approximately 500 m to the West of the seismic swarm.
We use synthetic waveforms calculated using reflectivity method (Wang, 1999) employing the same 1D velocity model as for ray tracing ( Figure 4). The model simplifies the true conditions and therefore produces comparably clear phase onsets. 15 However, the recorded data are also dominated by relatively little scattering and sharp onsets (Fischer et al., 2010). The effects of scattering on the final results will be discussed in greater detail in sections 5.
Hypocentral locations and origin times are taken from the double difference relocated catalog of Fischer and Michálek (2008).
All synthetic sources are double couple sources with mean strike, dip and rake set to 170 ± 10, 80 ± 10 and −30 ± 10 degrees, respectively, uniformly distributed in all three domains. The mean strike, dip and rake values are the predominant source types  attenuation in that zone is expected to be decreased with respect to the regional attenuation model. m -10500 m in Figure 4) has a Q p of 100 and Q s of 50. It is overlain by a more complex attenuation structure, characterized by higher Q values.
In order to mimic uncertainties in origin times, locations and velocity model travel times are perturbed by 10 ms, uniformly distributed. The uncertainties of the velocity model have an effect only in the source region since both rays of a couple traverse through the same overlaying velocity model. 5 The window length was 0.15 s for P and 0.3 s for S phases. The minimum allowed cross correlation of event couples in synthetic tests and later application to real data was set to 0.75. Signal-to-noise ratio (SNR) of a phase under consideration compared against a noise sample preceding the P phase has to be larger than 5 across the entire selected frequency band (after slight spectral smoothing to reduce effects of spectral notches). These two requirements efficiently reject outliers. The minimum allowed bandwidth is 10 Hz, which excludes all events with magnitudes of less than 0.5, given the magnitude dependent lower 10 frequency limit (f min = f c + 5.0Hz, where f c is calculated using equation 6). The bandwidth threshold stabilizes the linear fit to the spectral ratio as it limits the minimum number of data points. We evaluate Q p from vertical channels and Q s on North-South and East-West components and average results for each couple.
Data availability of the recorded dataset has been accounted for. Synthetic traces were only produced for an event if the recorded dataset contains data as well. All synthetic traces where convolved with the transfer functions of the WEBNET 15 stations to generate realistic velocity traces.  The next test depicted in Fig. 6 includes additive recorded noise. Data windows without seismic events in the field recorded data have been manually extracted and randomly added to synthetic traces to mimic realistic noise conditions. P phase results  In a next step (Fig. 7) we convolve synthetic Greens functions with realistic magnitude dependent source time functions. The applied source time function is half sine shaped where the slope of the high frequency spectral roll off is not dependent on the width of the applied pulse as can be seen in Figure 8 where normalized synthetic source spectra are depicted for different relative pulse widths. The vertical position of the spectral envelope changes with changing duration but the slope remains the   same for all depicted factorized pulse widths. Other than expected, this stabilizes results. This is a result of the pulse broadening which leads to a stabilization of MTM estimates as onsets become less transient.
The performed synthetic tests cannot reproduce waveforms in its full natural complexity. Nevertheless, they prove that the concept is capable to estimate attenuation of the anticipated region.   (Figure 10) explains the high similarity of waveforms observed for each swarm. We therefore also assume that source characteristics including rupture directivity effects are similar throughout each swarm cycle.

Application to North West Bohemia
By the time of the 2008 swarm the WEBNET stations were equipped with three component short period seismometers, except for station NKC located in the epicentral area which is a broad band station. All waveforms are sampled at 250 Hz. An example of a P phase onset recorded at station LBC is shown in Figure 9 together with the estimated amplitude spectra. A manual 5 revision of all event waveforms has been done to remove those which show indications of event doublets happening shortly after each other but not being indicated as such in the catalog. Spurious signal leftovers of a preceding event not necessarily cause high distortion of the fundamental P phase waveform and may thus not be removed by setting a cross correlation threshold. However, their effect lead to distortion at high frequencies of phase spectra and significantly increase the number of outliers during the analysis. The catalog of HypoDD (Waldhauser and Ellsworth, 2000) relocated events is comprised of 3841 10 events and their associated P and S phase picks. When applied to station LBC, a total of 641 couples were used which fulfill the requirements in terms of SNR, cross-correlation and geometrical constraints. Results of P phases evaluated at station LBC at station NKC and Q −1 s = 0.0037 (Q s = 270) at stations LBC, respectively. Both values are comparably large compared to Q p estimates from station LBC. A bias of these S phase attenuation measures is introduced by the P phase coda energy interfering 5 with S phases and therefore distorting the anticipated high frequency content.
In order to achieve a better understanding of the method's breakdown for P phase recordings of station NKC we disable the cross correlation threshold and scrutinize Q −1 against a multitude of parameters for station NKC and LBC. Figure 12 shows incidence angles of rays of event couples on the y-and Q −1 on the x-axis. Incidence angles are estimated using a 1D raytracing algorithm (Heimann et al., 2017). By definition the incidence angle is almost identical for both events of a couple. It becomes evident from figure 12, b) that larger incidence angles (> 8 degrees) show a tendency at station NKC to produce negative Q −1 while results from events with steep incidence angles produce positive Q −1 values. When compared to the same kind of plots 5 for station LBC no such trend is evident.
Station NKC is located at the northern edge of the swarm's epicentral region. Hence, incidence angle approximately correlates with latitude, indicating a location dependent problem. The depth sections of results from stations LBC and NKC ( Figure 13) show again that the attenuation is mostly positive (red) at station LBC (Fig. 13, a) and c)) in accordance with Figure 11. The distribution of attenuation at station NKC (Fig. 13, b) and d)) indicates a trend of decreasing Q −1 values from North to South.  Figure 14 shows P phase waveforms of first events of couples recorded at station NKC and LBC, filtered between 1 and 30 Hz for events northern (Fig. 14, left panels) and southern (Fig. 14, right panels) of the aforementioned gap. While the used filter frequencies are actually below the exploited frequency band used in the analysis these waveforms demonstrate that the waveform complexity is significantly higher for events recorded at station NKC than for those recorded at station LBC indicating that scattering plays a major role along the ray paths to station NKC. The average shape of the waveforms following which indicates that station NKC was situated closer to the nodal plane than station LBC. This aspect is in accordance with synthetic tests in chapter III which show worse performance of station NKC then station LBC for P phase measures.

Interpretation and Discussion
We present a newly developed method to estimate attenuation from spectral ratios of event couples. The short analyzed time windows are prone to spectral leakage which we mitigate by applying the multitaper approach. However, this method can only The given geometry and available data limit the perceptive field of the applied method mostly to the lower half of the seismic swarm (Fig. 2). Rays arriving at station LBC show a worse penetration of the focal zone compared to those arriving at station NKC (Fig. 2) where the takeoff angle is almost perpendicular to the normal orientation of the main rupturing fault. Therefore, the setting of sources with respect to station NKC is favorable from a geometrical point of view as it can be considered to have a higher sensitivity for attenuation in the source volume.
Synthetic tests show that the method is capable to reproduce average source volume attenuation of Q p = 100 and Q s = 50 5 given the 2008 earthquake swarm hypocentral locations. In noise free condition a precise result can be achieved for both, P and S phases at both significant stations NKC and LBC. With additive recorded noise, the distribution of Q P results broadens but still resembles the true attenuation with high precision at station LBC. Synthetic waveforms at station NKC suffer from weaker SNR at high frequencies indicating that the applied SNR-threshold of 5 is too optimistic. However, reducing this value would also reduce the number of data points and therefore have a negative impact on the statistical significance of results. Q s 10 results at station NKC are more robust then at station LBC. Both measures, Q s and Q p , improve when convolving synthetic waveforms with synthetic source time functions as this stabilizes the multitaper spectral estimates.
The application to recorded P phase onsets shows fewer negative Q p results at station LBC and than at station NKC which results in a distribution with a clearer offset with respect to Q p = 0. Waveforms recorded at station LBC are significantly higher correlated than those at station NKC where changing waveform polarities and high waveform complexity can be observed. We 15 hypothesize that rays arriving at NKC experienced relatively stronger scattering or that a nearby reflector creates multiples which interfere with signals recorded at NKC. In the latter case, the reflector would have to be situated in a location where interactions with signals arriving at LBC are weaker. Mousavi et al. (2017) assume a highly fractured medium in combination with accumulated free gas or fluids which would cause a high attenuation in the source region and could therefore explain our observations on low Q p . A 3D V p /V s tomography by Alexandrakis et al. (2014) identified a low V P /V S ratio body directly overlaying the focal zone. The increased waveform complexity seen at station NKC can be a result of waveform interaction with that body. The effects do not necessarily have an influence on results at station LBC where rays have different paths and take off angles.
Another influence may be rooted in the different families of focal mechanisms. Vavryčuk et al. (2013) reported three different 5 families of focal mechanisms for the 2008 swarm. While the slope of high frequency spectra should not directly be affected by the radiation pattern, there can be higher order effects like rupture propagation and rupture complexity. Dependent on the take off angle these rupture dynamics can affect the high frequency spectral roll off and therefore map into attenuation estimates (Kaneko and Shearer, 2015). P phase polarity changes at station NKC indicate that the station is located close to the nodal plane of the main rupturing fault. This circumstance can increase the effect of the aforementioned effects seen at station NKC.

10
If they differ systematically in the lower and upper source region, this can lead to biases in attenuation analysis due to the heterogeneous sensitivity across the fault plane. Still, we do not see such effects at station LBC and therefore speculate, that the dominating effect is the differing raypaths or a combination of both, raypath scattering and rupture dynamics.
Our findings in terms of P wave attenuation based on data from station LBC show similarly low values compared to results by Wcisło et al. (2018) who obtained Q p ≈ 120 for the source regions. Previous studies by Michálek and Fischer (2013)  Frequency bandwidth is a critical parameter which is limited mostly by the corner frequency of the recording setup and signal to noise ratio at high frequencies. Future plans of the Intercontinental Drilling Project (ICDP) include the installation of up to 4 borehole seismometers in NW Bohemia. It can be expected that our method will benefit from these measurements. Improved 25 signal to noise ratios allow to sample and exploit information at higher frequencies which will stabilize the spectral estimate.
Furthermore, higher sampling rates allow a better temporal (and therefore spectral) resolution of P and S phases. This will, in turn, also allow to use even shorter time windows. For the method discussed here, it would be favorable if at least one of these borehole stations will be situated in a location where a high number of ray path sharing couples can be found. The most sensitive region follows the NNE -SSW striking of the fault and concentrates in the North of the earthquake swarm (Fig. 3).

30
In late September 2017 the University of Potsdam installed a short period seismometer close to the Czech-German border in Oberzwota (red triangle, map 10) which is a favorable location. The station recorded 1000 samples per seconds for 62 days during a period of relative quiescence. Nevertheless, approximately 30 events were recorded in the swarm area with local magnitudes down to Ml=0. Despite the installation directly on top of the weathering layer the recordings showed signal to noise ratios larger than 5 at 120 Hz and above for smallest magnitudes. It becomes evident that even a surface mounted station 35 would allow to harness spectral information above the corner frequency of the WEBNET stations also for smallest magnitudes and which indicates that this will improve the resolution and robustness of our method once the ICDP borehole installations are operating.

Conclusions
Applying the source couple amplitude spectral ratio method to differential phase measures is an alternative to methods which 5 commonly exploit the lower frequency ranges. Theoretically, it is therefore able to achieve better resolution. Our synthetic study validates this. The geometrical constraints of this method require a high density of events as it is the case for natural earthquake swarms or seismic nests but also for hydrofracturing experiments.
The application to data from the 2008 North West Bohemia earthquake swarm indicates source region Q p < 100 based on data recorded at station LBC. The sensitive region measures only approximately 2000 * 500 * 500 meters in North, East and West 10 direction (Fig. 2). Results can therefore be considered of high spatial resolution. Nevertheless, the distribution of solutions significantly scatter and we see room for improvement e.g. through high frequency borehole recordings. We are not able to retrieve stable estimates at station NKC but instead see negative attenuation in southern and positive attenuation in the northern section of the swarm. P phase waveforms of the two sections show systematic differences at both significant stations which indicates a North-South structural difference. Furthermore, this effect does not inflict on measures at station LBC. Given the 15 fact that ray segments at NKC and LBC probe two different but directly neighboring media leads us to the conclusion that the fractured medium is highly concentrated along the source patch and that the surrounding medium can be considered much more dense or intact.
Author contributions. MK: Implementation, testing, evaluation and application to synthetic and recorded data, as well as paper writing.

Competing interests. None
Acknowledgements. This work is part of the HISS project which is funded by the DFG ICDP. We didn't remove this sentence as we think that this is a true feature. See reply no. 12. 12. Fig. 13: It would be more suitable to use a kilometer scale on this plot. Besides, the gap at 50.211 latitude appears enigmatic. I believe this plot should be station independent, because it shows coordinates of the events, so it should be visible at both the stations. Even if different events would be used at different stations, at least some indication of the gap should be visible also at LBC, provided there is some overlap between the event sets. Could you please show the vertical section of used relocations in order to identify the origin of the gap? We use a kilometer scale now and modified the figure to show the contribution of deeper and shallower sources as a vertical section in the upper panels. The former latitude gap at station NKC was partially due to a plotting problem. However, some gap is still visible and resulting from the spatial non-homogeneous distribution of the deeper sources (Panel b)).
13. P12: I am not sure if only scattering is responsible for the reported waveform complexity at NKC. It could be caused by different effects. One of these could be the proximity to the nodal line, where the P-onsets are usually of emergent nature and later arrivals are more visible. The data shown however look quite impulsive, which could indicate that the two mechanisms have different nodal lines, which not close to the NKC projection. Another reason could be overlap of waveforms with opposite polarity, which makes the image more complex. We actually consider different hypothesis to explain the complex signals at station NKC including scattering (P 13/13) but also focal mechanisms (P 13/5 ) and nodal lines (P 13/9). (Line numbers refer to the updated manuscript).
14. • P13/6: Here it would be helpful to discuss more the different ray direction and coverage of the focal zone for NKC and LBC as visible in Fig. 2, which could affect different sensitivity of the stations to attenuation. We extended the paragraph discussing the focal zone penetration with respect to Figure 2.
15. P14/1-8: It would be suitable to refer here to the study of Wcislo et al (2018) who obtained similar Qp and Qs using different method on the same data. Mentioning this in Conclusions is too late; BTW Conclusions usually do not contain citations. We now cite the work of Wcislo earlier in the discussion.