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⁻¹) 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 Q_p<100. The application to S phases fails due to scattered P-phase energy interfering with S phases. The Q_p 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.

1 Introduction

The intrinsic and scattering attenuation of the amplitudes of seismic waves is described by the dimensionless factor Q. The mapping of spatiotemporal changes of Q is an important step in seismology, since Q is controlled by temperature, rock porosity, fluid saturation and rock composition (Toksöz and Johnston, 1981). Hence, this factor may help to unravel the possible causes of fluid-induced earthquakes, or thermal anomalies in crustal regions affected by magmatic intrusions. For instance, northwest (NW) Bohemia is regularly affected by earthquake swarms lasting several days or weeks with thousands of recorded events with largest magnitudes up to M_l=4.4 (Fischer et al., 2014). The causes of the repeated earthquake swarms which occur in narrow focal zones remain under debate. Relative earthquake localizations are very precise because of the high waveform quality recorded with a dense local permanent network (Bouchaala et al., 2013). Different tomography studies revealed consistent figures of the 3-D velocity structures (Alexandrakis et al., 2014). The attenuation structure in the source region of the earthquake swarms is scarcely discussed. Some previous studies on whole ray-path Q exist and can be used for verification and benchmarking. However, the main aim of this study is to test whether the method developed here can enhance the resolution of near-source Q and therefore enable more robust conclusions on source dynamics and the role of fluids in the swarm cycle.

Several studies investigated the regional attenuation of northwest Bohemia by integrating along the full path from sources to receivers. Gaebler et al. (2015) estimated intrinsic and scattering attenuation of S waves (Q_s) by means of 14 selected events. Their frequency-dependent results indicate mean ${\stackrel{\mathrm{‾}}{Q}}_{\mathrm{s}}$ of approximately 1000. A study by Michálek and Fischer (2013) investigated source parameters and inferred a station-dependent, regional Q_p from P-phase spectra. They estimated mean Q_p ranging between 100 and 450. They also discuss effects of directivity on Q, concluding that the directivity has little influence due to the position of stations with respect to radiation patterns.

A tomographic study of northwest Bohemia done by Mousavi et al. (2017) indicated a regional average attenuation of approximately Q_p≈100 to Q_p≈300 and a pronounced highly attenuative source region where Q_p<100.

Bachura and Fischer (2016) employed two different methods to resolve the regional coda Q_c from the source volume to receivers. They used 13 selected events of the 2011 swarm and found a variation of Q_c between 100 and 2500 within the exploited frequency range of 1–18 Hz.

Recent work by Wcisło et al. (2018) used a newly developed differential attenuation estimation technique focused on the source region. The authors employed the peak frequency method which relates the half period of P pulses to attenuation. They also used a differential approach to map the inter-event attenuation using a single station (Nový Kostel; NKC) and found Q_p≈120 and Q_s≈80 in the source region.

Most previous Q studies focusing on NW Bohemia were inherently of low spatial resolution. This was, firstly, either because Q was 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. For example, Mousavi et al. (2017) used frequencies between 1 and 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 decreasing spectral slope at high frequencies ω above the corner frequency of the earthquake, often referred to as the ω² 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 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 6 to 13 October 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 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 to overcome this problem, Thomson (1982) proposed the multitaper method, which we employ using mtspec (Prieto et al., 2009; Krischer, 2016).

2 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)

$$\begin{array}{}\text{(1)}& {\displaystyle }{A}_{i,j}\left(\mathit{\omega }\right)=S_j\left(\mathit{\omega }\right)I_i\left(\mathit{\omega }\right)R_i\left(\mathit{\omega }\right)G_{i,j}\cdot e^{\frac{-\mathit{\omega }{t}_{i,j}^{*}}{\mathrm{2}}},\end{array}$$

where ω is the angular frequency. S_j(ω) describes the source spectrum and I_i(ω) the instrument response. R_i(ω) is the receiver 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:

$$\begin{array}{}\text{(2)}& {\displaystyle }{t}^{*}=\int Q^{-\mathrm{1}}/vds,\end{array}$$

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 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₀ and E₁ (in the following referred to as the first and second events of a couple) recorded at a station j (Fig. 1). We assume that their source spectra S₀(ω) and S₁(ω) resemble each other to a degree where the effect of random perturbations at high frequencies averages 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

$$\begin{array}{}\text{(3)}& {\displaystyle }\ln \left(A_{i,\mathrm{0}}\right(\mathit{\omega })/A_{i,\mathrm{1}}(\mathit{\omega }\left)\right)=\ln \left({\displaystyle \frac{G_{i,\mathrm{0}}}{G_{i,\mathrm{1}}}}\right)-\mathit{\omega }\mathit{\delta }{t}^{*}/\mathrm{2},\end{array}$$

with

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathit{\delta }{t}^{*}=t_{i,\mathrm{0}}^{*}-t_{i,\mathrm{1}}^{*}.\end{array}$$

This equation describes a linear relation with frequency-independent geometrical losses to the left of the negative sign in Eq. (3). The slope k of a line fitted to Eq. (3) can be used to derive the attenuation time t^* in between the two sources from which Q⁻¹ can easily be inferred using Eq. (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 ω² resulting 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 Bohemian seismic swarms. We use a relation proposed by Michálek and Fischer (2013) to calculate source radii r based on the moment M₀ of an event:

$$\begin{array}{}\text{(5)}& {\displaystyle }r=\mathrm{0.155}\cdot M_{\mathrm{0}}^{\mathrm{0.206}},\end{array}$$

where $M_{\mathrm{0}}=\mathrm{1.38}{M}_L+\mathrm{10.3}$.

Figure 1 Schematic illustration of the geometrical constraints and the employed parameters. The triangle represents a recording station at the surface. Attenuation is estimated for the traversing distance ray-path segment (D_t, green dashed line). Geometrical constraints respect the passing distance (D_p). The grey shaded area illustrates the Fresnel volume of the first event.

Figure 2 Ray segments with respect to stations Luby (LBC) and NKC. Source region seen from the west (a) and from the south (b). The colors indicate 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.

The resulting source radii r are then converted to f_c using

$$\begin{array}{}\text{(6)}& {\displaystyle }{f}_{\mathrm{c}}=k\mathit{\beta }/r,\end{array}$$

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 an additional 5 Hz with respect to f_c to account for uncertainties in M₀ 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 (Fig. 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 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, which is

$$\begin{array}{}\text{(7)}& {\displaystyle }{N}_{\mathrm{tapers}}=\mathrm{int}(\mathrm{bw}\cdot \mathrm{2})-\mathrm{1},\end{array}$$

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 Fig. 1. Ray tracing is done based on a 1-D 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₀ 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₁ onto the ray between E₀ and the station. We defined a minimum traversing distance of D_t≥1500 m to ensure that the signal of the second event is attenuated sufficiently to be detectable in the couple's spectral ratio.

Figure 3 Mapping potential locations of future seismic stations favorable to resolve Q at the swarm using the described method. The size of the blue points for each tested location represents the relative number of arriving shared ray paths from event couples. Red dots show the seismicity of the investigated swarm. Green triangles indicate locations of the local seismic network (WEBNET) stations.

Subsequent to geometrical preselection upper frequency limits of the analyzed bandwidth are potentially corrected to lower values dependent on the second Fresnel volume in between event E₀ and a station. The wavelength from which this frequency limit can be deduced is given as (Matsumoto et al., 2009)

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathit{\lambda }=D_{\mathrm{p}}^{\mathrm{2}}{\displaystyle \frac{D_{\mathrm{t}}+x}{D_{\mathrm{p}}xn}}.\end{array}$$

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₁ to the receiver (Fig. 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₁ is located within the Fresnel volume (grey shaded area in Fig. 1) for the entire analyzed frequency band.

With this approach, we get an estimate for the attenuation along the traversing distance (green dashed line in Fig. 1) and when repeated for a large number of event couples it 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.

3 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 northwest Bohemia. The size of the blue points in Fig. 3 represents the relative number of arriving shared ray paths from event pairs at possible station sites based on ray tracing through a 1-D layered model (Fig. 4). Largest numbers of pairs are expected along a north–south-striking patch which follows the striking direction of the main fault plane. However, in this case study, we use only those stations which provide continuous recordings for the investigated time span. These are stations Kraslice (KRC), Luby (LBC), NKC, Skalná (SKC) and Vackov (VAC). After geometrical filtering, we expect stations NKC and LBC to produce the highest number of couples since they provide continuous recordings for the entire time period and are in a favorable lateral location. Most other stations are located where no or a negligible number of event pairs are expected. Figure 2 shows the rays which penetrate the source volume and fulfill the geometrical requirements 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.

Figure 4 Synthetic velocity and attenuation model used for validation of the method. The seismogenic zone is marked by light blue. The attenuation in that zone is expected to be decreased with respect to the regional attenuation model.

Figure 5 Synthetic tests targeting Q_p=100 (a) and Q_s=50 (b) with noise-free data. No source time functions have been convolved. The correct values for Q_p and Q_s are indicated by the vertical dashed lines. The station color coding as given in the legend is used consistently throughout all following images.

Figure 6 Synthetic tests with the same setup as in Fig. 5 but with additive real recorded noise from analyzed stations. Again, traces have not been convolved with a source time function.

We use synthetic waveforms calculated using the reflectivity method (Wang, 1999) employing the same 1-D velocity model as for ray tracing (Fig. 4). The model simplifies the true conditions and therefore produces comparably clear phase onsets. 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 Sect. 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 $-\mathrm{30}\pm \mathrm{10}$^∘, respectively, uniformly distributed in all three domains. The mean strike, dip and rake values are the predominant source types stated by Fischer et al. (2014) which were retrieved based on polarity analysis of P phases. The seismogenic zone (depth 8500–10 500 m in Fig. 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.

Figure 7 Synthetic tests with the same setup as in Fig. 6. Synthetic traces were convolved with magnitude-dependent synthetic source time functions.

Figure 8 Normalized synthetic amplitude spectra of applied source time functions with different durations as a factor relative to the blue one (factor of 1). The slope of the high-frequency spectral roll-off is identical.

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.

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 frequency limit ($f_{min}=f_{\mathrm{c}}+\mathrm{5.0}$ Hz, where f_c is calculated using Eq. 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 contained data as well. All synthetic traces where convolved with the transfer functions of the local seismic network (WEBNET) stations to generate realistic velocity traces.

Figure 5 shows distributions of retrieved Q⁻¹ estimates from all event couples of the synthetic test using noise-free traces. In this test, as well as in the following test depicted in Fig. 6, traces have not been convolved with a source time function (i.e., they have impulsive source durations). The resulting distributions show some scattering solutions. Peaks in both cases (Q_p and Q_s) resemble the targeted attenuation model (dashed vertical line).

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 show a broadening of the distributions at all stations. While the distribution of station LBC still centers around the model value, the results of station NKC show only weak correlation with the correct model. This is a result of the location of station NKC close to the nodal plane of the dominant rupturing plane where smallest signal amplitudes are expected. S-phase results match the model at NKC but show strong scattering at LBC as a result of the interference with the added noise as well as the P-phase coda.

In a next step (Fig. 7), we convolve synthetic Green's 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 Fig. 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 their full natural complexity. Nevertheless, they prove that the concept is able to estimate attenuation of the anticipated region.

Figure 9 Two P-phase wavelets recorded at station LBC and their spectra from events with local magnitudes M_l=1.2 and M_l=1.4. The first two panels show the (a) first and (b) second events of the analyzed event couple. The grey shaded area in panel (c) indicates the used frequency band. The cross-correlation coefficient is 0.91 and the attenuation was in this case estimated as $Q^{-\mathrm{1}}\approx \mathrm{0.006}$.

Figure 10 WEBNET stations (triangles) and seismicity (blue circles) which occurred in northwest Bohemia during the investigated time period. Red triangle: temporary station installed by the University of Potsdam in October 2017. Stations in NKC and LBC are highlighted. All other stations produce no or an insignificant number of data points.

Figure 11 Attenuation results at stations LBC and NKC for Q_p (a) and Q_s (b). Median values are indicated by the overline Q.

4 Application to northwest Bohemia

Northwest Bohemia is a favorable case for testing our approach. Several focal mechanism studies on earthquake swarms in this region indicate dominant principle faults striking at 169 and 304^∘ (Vavryčuk, 2011), which have been active in different seismic sequences. Events occurring during a swarm tend to rupture on the same fault. This observation in combination with compactness of seismic clusters (Fig. 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.

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 broadband station. All waveforms are sampled at 250 Hz. An example of a P-phase onset recorded at station LBC is shown in Fig. 9 together with the estimated amplitude spectra. A manual 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 led to distortion at high frequencies of phase spectra and significantly increased the number of outliers during the analysis. The catalog of HypoDD (Waldhauser and Ellsworth, 2000) relocated events is comprised of 3841 events and their associated P- and S-phase picks.

Figure 12 Incidence angles plotted against $Q_{\mathrm{p}}^{-\mathrm{1}}$ of rays originating from event couples recorded at station LBC (a) and NKC (b). $Q_{\mathrm{p}}^{-\mathrm{1}}=\mathrm{0}$ is highlighted with a dashed line.

When applied to station LBC, a total of 641 couples were used, which fulfills the requirements in terms of SNR, cross-correlation and geometrical constraints. Results of P phases evaluated at station LBC (Fig. 11, left) have a median ${\stackrel{\mathrm{‾}}{Q}}_{\mathrm{p}}^{-\mathrm{1}}=\mathrm{0.019}$, equivalent to ${\stackrel{\mathrm{‾}}{Q}}_{\mathrm{p}}=\mathrm{53}$. The distribution shows some negative results which do not have a physical meaning and are related to noise in spectral estimates. Results retrieved based on data from station NKC are significantly more unstable, as Fig. 11a indicates. The distribution shows a large number of negative results. The median attenuation is ${\stackrel{\mathrm{‾}}{Q}}_{\mathrm{p}}^{-\mathrm{1}}=\mathrm{0.001}$, equivalent to ${\stackrel{\mathrm{‾}}{Q}}_{\mathrm{p}}=\mathrm{1000}$. Overall, 1404 couples where used in this case.

Attenuation evaluated for S phases shows almost zero centered distributions at both stations (NKC and LBC), which in turn indicates a significant number of negative and therefore unphysical measures. Median attenuation values are ${\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}^{-\mathrm{1}}=\mathrm{0.0023}$ ($\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}=\mathrm{435}$) at station NKC and $Q_{\mathrm{s}}^{-\mathrm{1}}=\mathrm{0.0037}$ ($\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}=\mathrm{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 with S phases and therefore distorting the anticipated high-frequency content.

Figure 13 Spatial distributions of attenuation at station LBC (a, c) and NKC (b, d). The source region is depicted as seen from the west. The origin of the coordinate system corresponds to the southernmost event. Panels (a) and (b) show average color-coded attenuation results of each used event. Panels (c) and (d) show the attenuation for each single event pair. First (deeper) events of each couple are represented by the x axes, whereas second events of each couple are represented by the y axes. Attenuation at station LBC is mostly positive. Subhorizontal seismicity gaps in panels (a) and (b) are caused by the set minimum distance (D_t) between two events.

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⁻¹ against a multitude of parameters for stations NKC and LBC. Figure 12 shows incidence angles of rays of event couples on the y axis and Q⁻¹ on the x axis. Incidence angles are estimated using a 1-D ray-tracing algorithm (Heimann et al., 2017). By definition, the incidence angle is almost identical for both events of a couple. It becomes evident from Fig. 12b that larger incidence angles (>8^∘) show a tendency at station NKC to produce negative Q⁻¹, while results from events with steep incidence angles produce positive Q⁻¹ values. When compared to the same kind of plots for station LBC, no such trend is evident.

Figure 14 P-phase onsets of couples' first events in the northern and the southern part of the swarm area. Panel (a) shows recordings at station LBC and (b) recordings at station NKC.

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 (Fig. 13) show again that the attenuation is mostly positive (red) at station LBC (Fig. 13a and c) in accordance with Fig. 11. The distribution of attenuation at station NKC (Fig. 13b and d) indicates a trend of decreasing Q⁻¹ values from north to south. This supports the hypothesis of a location-dependent issue. Figure 13d shows accumulated positive results related to first events at 1.4 to 2.0 km north (x axis), which coincide with the position of a small subcluster seen in Fig. 13b below a depth of 9.8 km. From 1.2 to 1.4 km north, the event couple distribution becomes more sparse. The separation of positive attenuation in the north and mostly negative attenuation in the south of that gap implies a systematic change in frequency content from two separated segments of the swarm occurring along ray paths from the source region to station NKC.

Figure 14 shows P-phase waveforms of first events of couples recorded at stations NKC and LBC, filtered between 1 and 30 Hz for events northern (Fig. 14a, c) and southern (Fig. 14b, d) 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 shapes of the waveforms following the first onset pulse differ in the northern and the southern sections. Furthermore, there are P phases with flipped polarities which indicate that station NKC was situated closer to the nodal plane than station LBC. This aspect is in accordance with synthetic tests in Sect. 3 which show worse performance of station NKC and then station LBC for P-phase measures.

5 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 suppress leakage caused by windowing effects. Leakage from source effects such as finiteness of the slip and the rupture are not smoothed as they are expected source characteristics and not filter artifacts.

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 able to reproduce average source volume attenuation of Q_p=100 and Q_s=50 given the 2008 earthquake swarm hypocentral locations. In noise-free conditions, 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 results at station NKC are more robust than those 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 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 3-D 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 takeoff angles.

Another influence may be rooted in the different families of focal mechanisms. Vavryčuk et al. (2013) reported three different 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 takeoff 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. If they differ systematically in the lower and upper source regions, 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 ray paths or a combination of both ray-path 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) investigated source characteristics in NW Bohemia and suggested station-dependent whole-path-integrated mean Q_p values ranging between 100 and 450. We find lower values which can be a result of hydration of the seismogenic zone. Haberland and Rietbrock (2001) also report highly increased attenuation (Q<100) within earthquake cluster regions and postulated that this could be related to hydration or partial melting. For instance, melt migration has been postulated from the size and migration pattern of earthquakes of the 2000 earthquake swarm (Dahm et al., 2009). On the other hand, Alexandrakis et al. (2014) interpret their results on velocity variations by dehydration processes. Our results deduced from station LBC for average attenuation are in line with previous findings pointing to high attenuation in the source volume.

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 four borehole seismometers in NW Bohemia. It can be expected that our method will benefit from these measurements. Improved 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 is concentrated in the north of the earthquake swarm (Fig. 3).

In late September 2017, the University of Potsdam installed a short-period seismometer close to the Czech–German border in Oberzwota (red triangle, map in Fig. 10), which is a favorable location. The station recorded 1000 samples per second for 62 days during a period of relative quiescence. Nevertheless, approximately 30 events were recorded in the swarm area with local magnitudes down to M_l=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 would allow to harness spectral information above the corner frequency of the WEBNET stations also for smallest magnitudes, which indicates that this will improve the resolution and robustness of our method once the ICDP borehole installations are operating.

6 Conclusions

Applying the source couple amplitude spectral ratio method to differential phase measures is an alternative to methods which 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 is the case for natural earthquake swarms or seismic nests but also for hydrofracturing experiments.

The application to data from the 2008 northwest Bohemian earthquake swarm indicates source region Q_p<100 based on data recorded at station LBC. The sensitive region measures only approximately $\mathrm{2000}\times \mathrm{500}\times \mathrm{500}$ m in the north, east and west directions (Fig. 2). Results can therefore be considered of high spatial resolution. Nevertheless, the distribution of solutions significantly scatters, 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 the 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 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.