Velocity ratio variations in the source region of earthquake swarms in NW Bohemia obtained from arrival time double-differences

Abstract

1 INTRODUCTION

Based on laboratory experiments and observations of earthquakes and volcanic unrest in Japan (Mogi 1963) discussed the different seismicity patterns and magnitude frequency relations observed for tectonic after- and foreshock sequences and earthquake swarms. He proposed that differences rely mostly on the degree of rock heterogeneity and the increase in crack density with time. While typical tectonic sequences comprise main event and a power-law decrease of aftershock activity, the seismic energy release in earthquake swarms often grows during the period of activity, or occurs in several bursts, and no outstanding principal event can be recognized. This clearly indicates some fundamental difference in the cause of swarms and tectonic earthquake sequences. However, the statistics and characteristics of swarm seismicity varies a lot for different cases and even within a single region (e.g. Hensch et al.2008), and other simple classification and generation models since 1963 are therefore not widely accepted until today.

Source processes of natural earthquake swarms comprise magmatic intrusion at depth (e.g. Matsumura et al.1991; Aoki et al.1999; Dahm and Brandsdottir 1997; Aoyama et al.2002; Hayashi & Morita 2003; Rivalta & Dahm 2006; Takeo et al.2006), fluid-related activity (e.g. Hainzl & Fischer 2002; Bräuer et al.2003; Hainzl et al.2006; Dahm et al.2008), possible gas and CO₂-induced activity (e.g. Hill 1996; Bräuer et al.2003), or possible hydrothermal activity (e.g. Fialko & Rubin 1999; Hensch et al.2008).

Human-related activities may also induce swarm-like seismicity, for example during controlled fluid-injection in gas-fields (enhancement fracturing) (e.g. Rutledge & Phillips 2003; Fischer et al.2009; Dahm et al.2010), hydrothermal stimulation operations (e.g. Evans et al.2005), thermal heating in mines (e.g. Jansen et al.1993; Köhler et al.2009; Becker et al.2010), during collapse of void space in mine areas (e.g. Malovichko et al.2009) or beneath reservoir impoundments because of the related pore pressure diffusion (e.g. Pandey & Chadha 2003).

More novel seismicity models, as for instance from Dieterich (1994), are able to explain well the different occurrence patterns and magnitude relations of earthquake swarms, and have for instance been used to infer stress steps and/or stress-rate changes in the swarm source region (e.g. Dieterich et al.2000; Toda et al.2002; Morita et al.2006). However, these models can only indirectly infer the physical processes causing swarms. So, one remaining challenge today is to better deduce the changes of the elastic parameter in the source region of earthquakes swarms. For instance, the ratio of P- to S-wave velocity, γ = v^P/v^S, may strongly vary for the different possible causes of earthquake swarms.

The classical Wadati diagram approach is one method to resolve v^P/v^S. Wadati (1928) invented his method to better resolve the trade-off between earthquake depth and origin time for deep earthquakes beneath Japan. The Wadati diagram is a linear regression between arrival time differences of P- and S-waves plotted versus the arrival of the P-wave. Theory predicts a linear relation as long as v^P/v^S is constant over the travel path. Even if this assumption is violated, the Wadati diagram is resolving an average v^P/v^S beneath the seismic network (e.g. Kisslinger & Engdahl 1973). Despite its potential strength, the Wadati method is today most often only applied to help the earthquake location problem.

Another approach to resolve v^P/v^S in the source region of earthquake cluster was proposed by Fitch when he implemented the master event relative location method (Fitch 1975). Although first study of this method was used to retrieve v^P/v^S for deep earthquakes, most of the later applications and developments aimed for improving the relative locations. Poupinet & Ellsworth (1984), Ito (1985), Scherbaum & Wendler (1986) and Maurer & Deichmann (1995) pointed out the potential of cross correlation methods to improve the accuracy of time-difference measurements, and Slunga et al. (1995) and Waldhauser & Ellsworth (2000) introduced double-time differences and multiple master events to further improve the relative locations.

Lin & Shearer (2007) developed a method to estimate local v^P/v^S values from precise P and S differential times. Aim is to estimate near source and near fault rock properties, for example by widespread computation of near-source v^P/v^S ratios from large earthquake catalogues and interpreting these together with seismic tomographic models. An application to southern California indicated slightly reduced v^P/v^S ratios in some earthquake source regions and led to the speculation of the presence of water-filled cracks (Lin & Shearer 2009).

Our study combines the precise relative locations of microearthquakes with the robust double-difference Wadati method in order to resolve spatial and temporal variations of v^P/v^S in the source region of earthquake swarms. The methodical approach is very similar to that of Lin & Shearer (2007), expect that we formulate the method in parallel to the standard Wadati method, which enables simultaneous determination of the velocity ratio on two different scales: the scale of the seismic network and the local scale of the focal volume. Besides, we implement a least median of squares regression to reduce the possible effect of outliers and apply the method to fluid-related earthquake swarms. v^P/v^S ratios and their variations shall help to elucidate the source processes for earthquake swarms, and to possibly distinguish between gas- or fluid-induced swarms and thermal or mass-loading induced stress changes.

In the first section, we briefly derive the formulas and demonstrate and test our approach for a typical swarm configuration with a synthetic data set. The second and third chapter reports first applications for three midcrustal (8–10 km) earthquake swarms in NW-Bohemia. We observe strong temporal variations of v^P/v^S and discuss the possible causes.

2 METHOD

2.1 Theory

Figure 1.

Open in new tabDownload slide

Setup of the synthetic data example. The numerical grid covers a cube of 20 × 20 × 10 km with a grid spacing of 0.5 km. The simulated earthquake swarm consist of 20 events (red filled circles) within a inner cube of 2 × 2 × 2 km (source region defined by solid lines). 12 stations are placed at the surface in a distance ranging from about 6 to 12 km. The P and S wave was set to α₀ = 5.5 km s⁻¹ and β₀ = 2.9 km s⁻¹, respectively. The ratio γ = α/β is ≈1.9. To simulate an anomalous source region at the cube the S velocity was increased to β₁ = 3.6 km s⁻¹ while the P velocity was kept constant, α₁ = α₀, so that γ ≈ 1.5 in the source region.

How does eq. (4) change if an anomalous source region of radius r₀ with β = β₁ is introduced (alternatively for α). This leads to similar equation as (4), with slope γ unchanged (γ = α₀/β₀) and with a slightly different constant |$\delta _i =(1 -\gamma )t^0_i +\frac{r_0}{\beta _0}( 1 -\frac{\beta _0}{\beta _1} ) + \gamma \langle t^P \rangle$|⁠, which results in a ‘vertical’ shift of ‘Wadati line’. If this shift is estimated and considered for each event, a multiple event method may improve the estimate of the undisturbed γ of the scale of the seismic network.

Next, we introduce the double-difference version of the Wadati method by assuming that a tiny clustered earthquake swarm has occurred and that all stations are far away from the source region of the swarm, so that rays to a station from the individual swarm events can be assumed to be parallel spokes emerging the source region in direction |${\bf \hat{n}_j}$|⁠.

2.2 Inversion technique

The inversion of the double-difference data is done similar. For each γ_1k, the first step is to substract the event-pair dependent offset Δδ_lest, and the best solution γ_1best is either calculated for the L₁ norm or by using the LMS method. For large number of events a reduced number of selected event pair permutations may be sufficient in order to reduce the computational work.

According to (4), the parameter γ_best defines the average velocity ratio across the network, and according to (7) γ_1best represents the velocity ratio across the source volume. The approach has been compared to the method of Lin & Shearer (2007) and we retrieve identical results for noise-free data or data with additive, random noise. It may lead to different results if the sample size is small and outliers are present.

2.3 Synthetic test

We calculate traveltimes for the model in Fig. 1 by using a Finite Difference Eikonal solver (Podvin & Lecompte 1991). The ratio γ = α/β was set to γ₀ ≈ 1.9 for the homogeneous media (case HOM). Another experiment considers a small cubic, inhomogeneous region at the centroid of the earthquake swarm (Fig. 1), where γ is altered to a value of γ₁ ≈ 1.5 (case INH). Most simulated earthquakes occur within the anomalous region.

The inversion using the multi event, single-difference Wadati method (eq. 4) is illustrated graphically in Fig. 2 for the case HOM (left-hand side) and INH (right-hand side). Random noise with a standard deviation of 0.1 s for S and 0.08 s for P arrivals was added to theortical arrival times. The panels show the graphical form of eq. (5) using a velocity ratio of 1.7 for reduction. Due to the layout of the example the data points in Fig. 2 are clustered around two specific distances, that is P arrival times. Although the inhomogeneous source region affects the t^S − t^P time at each station, the best solution (slope) in the modified Wadati diagram is nearly unchanged to the background ratio of 1.9. The effect of the common path of rays in the anomalous source region (case INH) results in a ‘vertical’ shift of ‘Wadati line’ by about +0.1 s when compared to the HOM case. Using eq. (4) and demeaning the shift empirically for each event results in a smaller scatter of the modified Wadati diagram. A careful inspection of Fig. 2 for the inhomogeneous case indicates that the reduced arrival times of the point-cluster at − 0.5 s align along a smaller slope of about 1.5, while the overall trend of the data points resolves the background value of 1.9. The anomalous source region is here already indicated in the standard single-difference Wadati diagram. Similar features can be identified in some real data cases below.

Figure 2.

Open in new tabDownload slide

The modified single-difference Wadati diagram as defined in eq. (5). A reduction ratio of 1.7 was used. (a) The homogeneous (HOM, velocity ratio γ = 1.9) and (b) the inhomogeneous case (INH, background γ = 1.9 and anomalous source region with γ = 1.5, see Fig. 1). The horizontal dashed line represents a theoretical velocity ratio of 1.7. The dashed lines with positive and negative slopes indicate trends expected for a ratio of 1.9 and 1.5, respectively. The resulting velocity ratios are indicated by solid line.

The double-difference Wadati method (eq. 7) is illustrated in Fig. 3. Again, the standard deviation of additive noise was 0.1 s for S and 0.08 s for P arrivals. The spread of δt^Pd times is smaller than for the single-difference method, since the spread of the earthquake cluster defines arrival time differences and is much smaller than the distance range of stations. For both, the homogenous and the anomalous case, the double-difference method retrieves well the local velocity ratio of the source region. It is only sensitive to the velocity ratio in the source region and not influenced by the background model.

Figure 3.

Open in new tabDownload slide

The modified Wadati diagram based on double-difference arrival times according to eqs (6) and (7), for the homogeneous (a) and inhomogeneous cases (b). A reduction ratio of 1.7 was used. The gamma value of the best result represents the velocity ratio of the source region (denoted by γ₁ in the text). See Fig. 2 for further details.

Fig. 4 summarizes the results of the inversion under noise for model INH using the single-difference and double-difference Wadati method. The normal distributed noise of S arrivals was incrementally increased from 0 to 0.2 s. The noise of P arrivals was always 80 per cent of the noise of S arrivals. Further, 20 arbitrary chosen outliers (noise with standard deviation of 0.2 s) were added to S-wave arrivals. Different regression models have been used in order to study the robustness against normal distributed noise and outliers. Fig. 4 shows that the single-difference method resolves stably the background ratio α₀/β₀, even under large noise of 0.2 s and independent of the regression norm used. The estimated uncertainty of velocity ratios increases with increasing noise, but is generally relatively small. It is non-zero for zero random noise since outliers were added to the data set. In the case of double-difference inversion, the L₁ and LMS inversions are robust against outliers and retrieve unbiased source region ratios if the noise standard deviation is below 0.15 s. The formal errors of the inversion are generally larger than for the single-difference method. Under largest noise the best results indicate a small systematic bias towards smaller ratios. Trial and error testing indicate that the bias towards smaller values is mainly controlled by noise in P arrivals and not in S arrivals. If the P arrival noise is set half of the amplitude of the S arrival noise, the bias vanishes for the plotted range of values.

Figure 4.

Open in new tabDownload slide

Result of robustness tests against outliers for both the traveltime difference (dt) and double-difference (dd) regression. Squares stand for single-difference and circles for the double-difference Wadati method. The given uncertainty for L₁ is the standard deviation of the variance of the sample mean, calculated from the residuals after the best solution is estimated. Error bars for LMS are of similar size and are omitted to improve the readability. The synthetic data were simulated for the inhomogeneous model (INH) and by using the L₁-norm (red symbols) and LMS (green symbols) inversion approach. We simulated 20 outliers S-wave arrival times, arbitrary selected with a normal distributed additive S-arrival time with standard deviation of 0.2 s. Increasing additive normal distributed noise was added to all theoretical arrival times.

The synthetic tests performed here are extreme. Noise of S and P arrivals in local distances can often be expected smaller. The number of outliers may be smaller and of smaller value. The geometry under testing is representative for a local network and tiny cluster configuration as plotted in Fig. 1. Larger earthquake clusters are expected to have larger noise breakdown points.

3 APPLICATION TO NW BOHEMIA EARTHQUAKE SWARMS

The NW Bohemia/Vogtland earthquake swarms region is located in Czech Republic—German border region (Fig. 5). Since 1991 earthquake swarms clustered to 90 per cent in the region of the small city Novy Kostel (close to station NKC, Fig. 5) where earthquakes typically occur in a depth range between 6 and 12 km (Fig. 6). Hypocentres occur in a crystalline unit consisting of granite body on the west and mica schists on the east (Havir 2003; Chab 2007). No deep borehole or seismic reflection survey data are available to derive further details on the focal zone structure. The KTB deep borehole, 50 km to SW, has encountered uniform amphibolite metamorphic rocks with free water in open fractures down to its maximum depth of 9.4 km (Wagner et al.1997). Single swarms typically consist of several hundreds to thousands of earthquakes with magnitudes mostly below M_L < 3.5. Since 1990, four major earthquake swarms occurred in 1997, 2000, 2008 and 2011 (Fig. 6) and were recorded by the permanent local seismic networks. An overview to the microearthquake activity near NKC is given in for example Fischer & Horalek (2003); Neunhöfer & Hemmann (2005); Fischer et al. (2010). In this study, we use precise arrival time readings of seismic data recorded by the WEBNET network within an epicentral distance of about 15 km, which surrounds the Novy Kostel focal zone. The WEBNET consists of 13 short-period stations (1 or 2 s, three-component seismographs) sampled at 250 hz. The estimated error of the hand picks is ±1 sample (4 ms) for P waves and ± 2 samples (8 ms) for S waves. Several selected stations are situated directly in the areas of earthquake clusters to fix the hypocentre depth. For the investigated swarms, the epicentral distance of the closest station NKC ranged from 0.5 to 3 km. More details about the seismic network could be found in Fischer et al. (2010).

Figure 5.

Open in new tabDownload slide

Epicentral map of the NW-Bohemia earthquake swarms. Earthquakes of the year 2000 are plotted by red-filled circles. The used permanent seismic stations (triangles) are within a distance of about 15 km from the swarm area. Solid lines declare faults, the dashed line outlines the Cheb basin. Mofettes and Holocene volcanic structures (Komorni Hurka and Zelezna Hurka) are indicated.

Figure 6.

Open in new tabDownload slide

Hypocentres of swarm earthquakes at Novy Kostel are projected in three orthogonal planes (Master-event locations; figure adapted from Fischer et al.2010). Colours indicate the occurrence time between 1991 and 2011. The activity of specific swarm periods can be recognized in the magnitude over time panel, where they pop-up as energetic peaks (bursts).

The region is also known for its diffuse CO₂ degassing often manifested by moffets, mineral springs and gas vents (Weinlich et al.1998; Heinicke & Koch 2000; Bräuer et al.2003). The anomalous high concentration of mantle-derived ³He indicates magmatic reservoirs at mantle depth (e.g Bräuer et al.2007). Three surface manifestations of young volcanos (Komorni Hurka, Zelezna Hurka, and nearby newly discovered maar structure, all active during Holocene, see, e.g. Mrlina et al.2007) are located in a distance of 15 and 25 km to NKC (Fig. 5). Volcanic rocks at the surface indicate nephelitic magmas of small viscosities with unknown concentration of CO₂ or other volatiles (Mrlina et al.2007).

Monitoring of CO₂ degassing and its possible relation to the swarm activity shows ambiguous results. The small earthquake swarm in 1994 (M_L ≤ 1.8) was followed during the next months by clear anomalies of the gas contents and isotope ratios of ³He/⁴He and δ¹³ C in the mineral springs in Bad Brambach. On the other hand, the much stronger swarm in 2000 (M_L ≤ 3.3) have shown only weak pre-seismic decrease of ³He/⁴He ratios (Bräuer et al.2008) and no gas anomalies were observed in relation to the 1997 swarm (M_L ≤ 2.9). This indicates a high level of complexity in the relation between stress changes accompanying the seismogenic process and crustal fluids that are represented in the seismogenic depths primarily by water and CO₂. The upper mantle origin of carbon dioxide documented by isotope studies suggests CO₂ activity in fluid-permeated fault zone.

We applied several data pre-processing steps to remove outliers and to improve data quality. This is especially important to obtain reasonable uncertainty estimates. First of all, the minimum number of common stations with P and S picks of event pairs is set to 7 (or 9 and 11 for the 2000 and 2008 swarms, respectively). The higher the number of stations involved, the better is the estimate of δ_i and Δδ_i - this helps to reduce the variance of the solution. In order to remove the possible outliers in picking we require that the P and S arrivals align along a Wadati-line with a residual variance after a linear regression of less than 0.15 s; events with larger variance are omitted. We further added quality criteria for the estimation of Δδ_i to identify possible outlier event pairs. Additionally, the data set is limited to event pairs with a demeaned arrival time differences smaller than 0.35 s, in order to limit the spread of the seismicity cloud to remove events which occurred far from the swarm cluster.

3.1 1997 Earthquake swarm

For the 1997 swarm intervals were chosen to ensure more than 20 events (Fig. 7). The inversion of single-time differences reveals the average velocity ratio in the crust of about α₀/β₀ = 1.69 (Fig. 7). This ratio remains constant during the whole year 1997. Most interesting is the ratio retrieved from double-time differences, which represents the velocity ratio in the source zone of the earthquake swarm. α₁/β₁ is small at about 1.43 ± 0.15 in the beginning of the year when the swarm 1997 occurred and increases to nearly normal crustal values in the post-swarm period (Fig. 7). The strong velocity anomaly in the beginning of the activity has been tested by different methods (L₁, L₂ and LMS) and is a stable feature.

Figure 7.

Open in new tabDownload slide

Application of the wadatiDD inversion (L₁ norm) to sequences of the 1997 swarm. Only event-pairs with seven or more common stations were used. Histogram (top panel) shows the weekly number of events (counts) passing the quality criteria for the single-difference inversion. Dashed lines indicate the selected time intervals. Magnitudes of those event are plotted (open circles) which were selected as suited event-pairs for the double-difference inversion. The numbers (bottom panel) gives the number of event pairs in each time interval. The average ratio γ in the crust, α₀/β₀, is indicated by blue filled squares. The ratio associated with the narrow source region of the earthquake swarm, α₁/β₁, is retrieved from double-differences and given by red filled circles.

The scatter plots of single- and double-time differences used for inversion are shown for two time intervals at the beginning (Fig. 8) and the end of the activity (Fig. 9). The average crustal velocity ratio (background ratio) is very well resolved from the narrow band of data points along the line with slope γ₀ ≈ 1.69. The scatter is slightly higher for the double-time difference inversion. In the standard Wadati approach (single-differences) in Fig. 8, the subclustered data points already indicate a trend towards smaller (negative) slopes. A similar feature has been identificed in the synthetic tests and was indicating a source region with anomalous low-velocity ratios.

Figure 8.

Open in new tabDownload slide

Arrival time differences considering earthquakes in the time interval 1997 January 1–30 (day 1–30). Single-difference inversion results (left-hand side) according the eq. (3) and double-difference inversion results (right-hand side) according to eq. (9). A reduction ratio of 1.7 was used. The labelled numercial value denotes the best ratio of the inversion (represented by the solid line). The dashed lines with positive and negative slopes indicates theoretical slopes for γ = 2.1, γ = 1.3, respectively.

Figure 9.

Open in new tabDownload slide

Arrival time differences considering earthquakes in the time interval January 30 to 1997 December 31 (day 30–365). See Fig. 8 for further details.

The double-difference Wadati plots in Figs 8 and 9 show two different trends for the selected time intervals. For the first time interval at the beginning of the earthquake swarm activity, we retrieve γ₁ ≈ 1.43 ± 0.15 (Fig. 8, right-hand panel). The time period after 1997 January 30 is represented by a ratio of γ₁ ≈ 1.56 ± 0.04 (Fig. 9, right-hand panel). The hypocentres of the subcluster in the beginning and final phase of the swarm are projected in a vertical fault section (Fig. 10). The spread of the subclusters is very similar and geometrical or spatial occurrence patterns therefore cannot explain the observed strong temporal variatons of γ₁.

Figure 10.

Open in new tabDownload slide

Hypocentre location (grey circles) and spread of subclusters in different time intervals (black circles) of the swarms in 1997 (top panel), 2000 (middle panel) and 2008 (bottom panel). Events were projected in a vertical fault plane section. The time intervals in the left-hand panels correspond to early phases of swarm activity, where small v^P/v^S ratios are retrieved. The right-hand panels show intervals with normal v^P/v^S ratios.

3.2 2000 Earthquake swarm

Fig. 11 shows inversion results for time sequences of the 2000 earthquake swarm. The 2000 swarm was unusual because of its long duration lasting from August to December. It occurred in two major bursts, one from August 27 to September 26 (day 240–270) and the other from October 9 to November 21 (day 283–326). In comparison to 1997, the number of usable stations and data strongly increased. Therefore, we could select a threshold of nine or more stations for the double-difference method to further improve the reliability of results. The average velocity ratio of the crust was again very stable at γ₀ ≈ 1.69. The velocity ratio of the source volume, as derived by the double-difference inversion, was slightly reduced in the beginning of the first activity phase starting on day 238, recovered to normal values in the later phase of the initial activity but decreased again to an anomalous-low value in the beginning of the second activity phase at about day 280 (Fig. 11).

Figure 11.

Open in new tabDownload slide

Application of the wadatiDD inversion to sequences of the NW Bohemia 2000 swarm. Only event-pairs with nine or more common stations were used. The double-difference inversion was not applied to the December interval, because only one event pair passed the quality filter. See Fig. 7 for further explanations.

This is clearly demonstrated by the scatter plots for two complementary time intervals (Fig. 12). During the first activity phase from day 240 to 255 (August 27 to September 11) γ₁ ≈ 1.7 and resembles the regional value of γ₀ ≈ 1.69 (Fig. 12, left-hand side). However, during the onset of the second activity phase from day 283 to 297 (October 9 to 23) the velocity ratio in the source volume is strongly decreased and aligns with a slope of γ₁ ≈ 1.38 (Fig. 12, right-hand side), although the geometry of the network and swarm location was unchanged. The spread of Δt^P and thus the uncertainty of the inversion is similar in both time intervals. The γ₁ anomaly in the beginning of the second activity phase of the 2000 swarm is slightly smaller than for the 1997 swarm, and is much better resolved here. The position and extent of the subclusters in time intervals are similar and partly overlapping (Fig. 10). Changes of the position of the used events are therefore very unlikely to explain the temporal variations of γ₁

Figure 12.

Open in new tabDownload slide

Arrival time double-differences for the time interval 2000 October 9–23 (day 283-297, left-hand panel) and 2000 November 6–21 (day 311–326, right-hand panel) of the 2000 swarm. The solid lines indicate the retrieved slope associated with the labelled value. The dashed subhorizontal lines represent theoretical slopes for γ = 2.1, γ = 1.3, respectively.

3.3 2008 Earthquake swarm

Results of the 2008 earthquake swarm are plotted in Fig. 13. The 2008 swarm was recorded by a denser network of stations, and thus more events could be included in the analysis and smaller time intervals were chosen. We consider for 2008 a threshold of at least 11 common stations. The source region velocity ratio was anomalously low at the beginning of the swarm (Fig. 13), similar to the 1997 and 2000 swarm. However, the recovery to regional, normal ratios lasted few days only. The subcluster in the first time interval from October 6 to 8 is more focused than the one from time interval from October 11 to 13 (Fig. 10). The corresponding scatter plots of the modified Wadati diagrams are shown in Fig. 14, where a clear difference in the overall trend of time differences is visible. The least norm fit indicates a ratio of γ₁ ≈ 1.33 of the tiny cluster during the first interval (Fig. 14). The ratio retrieved for the second interval is γ₁ ≈ 1.58 and for subsequent intervals converging more and more to the background value of about 1.69. Fig. 14 (right panel) shows the scatter plot for a later interval, where data align nearly horizontal in the modified Wadati diagram, which corresponds to a ratio of 1.7. Fig. 15 shows seismograms at stations NKC, VAC and KRC at the whole used range of epicentral distances from 1 to 15 km for a weak (M_L 0.3) and medium-sized event (M_L 1.1). Clear sharp P and S phase onsets are visible. Phase identification is here and in general unambiguous and signal-to-noise ratio allows for precise arrival time picking with errors of ±1 sample for P waves ( ± 4 ms) and ±2 samples for S waves ( ± 8 ms) (see also Bouchaala et al.2013). Waveforms are similar, and a systematic bias from picking can therefore be excluded. The seismograms have in general a very high signal-to-noise ratio and quality, and do not show any unexpected peculiarities.

Figure 13.

Open in new tabDownload slide

Application of the wadatiDD inversion to sequences of the 2008 swarm. Only event-pairs with 11 or more common stations were used. Histrograms are calculated for daily number of events. See Fig. 7 for further explanations.

Figure 14.

Open in new tabDownload slide

Arrival time double-differences for the time interval 2008 October 6 to 7 (day 279–280, left-hand side) and 2008 October 11 to 12 (day 286–287, right-hand panel) of the 2008 swarm. The solid lines indicate the retrieved slope associated with the labelled value. The dashed lines represent theoretical slopes for γ = 2.1 and γ = 1.3. Data from only five stations are plotted, VAC, LBC, SKC, NKC, ZHC.

Figure 15.

Open in new tabDownload slide

Seismograms sampled at 250 sps (velocity, vertical = Z, radial = R, transversal = T) at stations NKC (2-km distance), VAC (5-km distance) and KRC (15-km distance) of a M_L 0.3 event occurred during the 2008 swarm on October 31, 18:31 and a M_L 1.1 event occurred on October 6, 03:47. The arrival times of direct P and S waves are indicated in zoomed Z and T traces.

4 DISCUSSION

The double-difference Wadati method was for the first time applied to mid-crustal, fluid-driven earthquake swarms in order to study possible changes and variations of the γ₁ = v^P/v^S velocity ratio in the source region of these swarms. We could find during different short phases of the swarm activity velocity ratios as small as γ₁ ≈ 1.33 during the initial phase of swarm activity.

4.1 Velocity ratio during swarms

Reduced velocity ratios were mostly observed in the beginning of activity bursts. For instance, during the 1997 swarm in NW Bohemia (M_Lmax = 2.9) the velocity ratios of about 1.42 were indicated. This period was associated with intense activity concentrated along a small patch of the fault of the size less than 0.5 × 0.5 km. Afterwards the concentrated activity ceased and occurred in a broader area of the fault zone. During the later phase of the swarm, the velocity ratio in the source region nearly recovered to normal values of γ₁ ≈ 1.69. The much larger 2000 swarm (with several M_L = 3 + events) started with only slightly reduced velocity ratios in the source region during the first period, which quickly recovered to normal values. This first activity period lasted 4 weeks and was interrupted for about 2 weeks. However, coincident with the beginning of a second strong activity phase after 45 days the velocity ratios abruptly dropped to γ₁ ≈ 1.38, and recovered again to normal values during the next 2 weeks. The 2008 swarm, which was the largest one with few M_L = 3.5 + events, indicated source volume velocity ratios of about 1.33 in the very first starting phase, but recovered much quicker to normal values after only 5 or 6 days. The fast recovery may be related to the shorter duration of the 2008 swarm compared to the 2000 swarm (Fischer et al.2010). As illustrated in Fig. 10, the hypocentres occurring within time intervals with different velocity ratios occur in areas of fault zone that partially overlap. This is namely the case of γ₁ ≈ 1.33, which was resolved for the period from October 6 to 8, when the activity occurred in a relatively small area that totally overlaps the activity from 2000 November and from 2008 October 11 to 13, which showed normal velocity ratio of about 1.69. The different time/space occurrence of anomalous ratios is thus probably not a structural feature of the fault zone, but may rather reflect complex swarm evolution and background processes like the presence of gases under high overpressure. However, the first question concerns the reliability and significance of the estimated values.

4.2 Reliability of the velocity ratios

The suggested method is simple and does dot require the hypocentre locations or velocity model to be known. The method thus seems to be robust against structural features and station geometry and therefore possibly suited for even routine analysis. It has been developed independently of the method by Lin & Shearer (2007) and solves the inverse problem under slightly different norm and demeaning approach. However, both methods lead to comparable results in our synthetic data tests under presence of noise and outliers. In our method, we simultaneously determine two kinds of the v^P/v^S ratio using the same inversion technique: the average ratio of the crust across the area of the seismic network and the local ratio associated to the focal volume.

The choice of a robust norm for inversion is of special importance when using the double-difference data. The reason is that arrival time differences require a higher accuracy of picks, since travel-path differences in the source region of an earthquake swarm are usually small. The relative error of the input data may thus be large. For instance, if the extent of the source region is about 1/10 of the average distance from the events to the stations, the traveltime measurements should be 10 times more precise than single phase picks to arrive at the same resolution for the double-difference method in comparison to a standard Wadati diagram. For example, we may select an event pair with an interevent distance of 1 km, which is recorded at one station. The S-wave velocity shall be unchanged, and we want to resolve a γ change from 1.7 to 1.4. This would then correspond to a drop in P-wave velocity from 6   to about 4.9 km s⁻¹, leading to an expected change of the P arrival time difference of about 40 ms (compare also with Fig. 12). This would cover 10 samples, that is about five times the nominal picking accuracy. If the interevent distance would decrease to only 500 m the expected time difference is only 2.5 samples, which may indicate that an anomalous region of a length of 500 m is already close to the limit of our resolution. However, since we employ a multiple event method and use for instance more than 10 stations for each event, the number of independent input data N is often in the order of several hundreds to several thousands. The uncertainty of the estimated γ is then expected to be reduced by factor of |$1/\sqrt{N}$|⁠, which improves the resolution. Fig. 16 shows the average residuals as a function of γ to better visualize the relative uncertainty in the studied time windows. At least for the 2000 swarm the low γ section has a similar deep minimum as the normal γ section. Using synthetic data, we tested the effect of noisy input data for swarm and station geometries similar to the studied ones. Assuming the typical spread of the P-wave arrival time double-differences in the range of Δt^P ≈ ±0.4 s, (similar to, e.g. the Δt^P of the swarm 2008 in Fig. 14) the standard deviations of Gaussian distributed arrival times up to 0.15 s were acceptable. This is well above the picking error of arrival times of the WEBNET data, which is estimated as one sample (4 ms) for P and two samples (8 ms) for S waves (see also Bouchaala et al.2013). However, errors of real data are possibly not Gaussian.

Figure 16.

Open in new tabDownload slide

Station averaged residuals are plotted versus γ₁ for the time intervals in Figs 9 (1997, left-hand side), 12 (2000, middle) and 14 (2008, right-hand side). The full and dotted lines indicate γ₁ for the first and second time intervals shown in the corresponding figures.

Another possible source of uncertainty of the resulting γ is the event-dependent offset Δδ_i, see eq. (9), which corrects for the unknown origin time. As indicated in some scatter plots (e.g. Fig. 12, right-hand panel) this offset is possibly not always correctly determined for each event. However, the influence of such a possible offset on γ₁ is minor due to only exceptional occurrence of outliers in Δδ_i; note that the number of input data ranges in the order of 10³ or more. Additionally, the robustness test shows that the robust least median of squares or L₁ norm regression is capable to suppress these effects. On the other hand, synthetic testing indicated that strong random noise in P arrival times might bias the velocity ratio γ₁ to smaller values. However, the arrival time picking of P onsets is unambiguous and was the most accurate one.

Despite this optimistic result, the resulting very small or even negative values of the Poisson's ratio (γ < 1.4) might be subject to critical view. For instance, the anomaly of the 1997 swarm might be questioned because of the small volume sampled by earthquakes, which is related to tiny arrival time double-differences. Additionally, the minimum number of stations for 1997 was only seven and thus relatively small. However, neither of these are the case of the 2000 or 2008 swarms where γ ≈ 1.38 for the minimum number of stations of 9 and 11 and Δt^P covering a wide range from −0.35 to 0.35 s. Also, the magnitudes of selected event pairs cover a similar range for high γ and low γ intervals and between the different swarms, so that a predominance of small magnitude events can be ruled out to explain the anomalies. This, together with the fact that for similar Δt^P spreads both normal- and low-velocity ratios are found and that station geometry, waveforms and/or activated cluster patches are all unchanged during the evolution of the swarms, indicates that the temporary reduced, low-velocity ratios are not artefacts.

4.3 Physical interpretation of the low-velocity ratios

What are the possible physical causes for a strongly reduced v^P/v^S ratio?

Figure 17.

Open in new tabDownload slide

Predicted relationship between γ and porosity Φ for brine-filled and gas-filled sandstone. Eq. (11) was used with a shear and bulk modulus of 30 and 50 GPa, respectively, and a Biot-Willies parameter of α = 1. The compressibility ratio of brine and gas is indicated in the legend.

The physical effect is that the S-wave velocity is carried by the skeleton of the porous rock and remains unchanged, while the P-wave velocity is affected by the fluid-replacement in the pore space and strongly reduced.

Since the Novy Kostel swarms occur at the intersection of two or more fault zones at a depth of 7–10 km, we assume that the rock there is fractured and can be represented by a porous model with small porosity. We postulate that the reduced value of γ₁ during some phases of the swarm activity is caused by the spontaneous transition of fluid to gaseous phases. This may be triggered by a sudden magma-water contact or by earthquake-related ground shaking under the presence of oversaturated gas.

Dahm et al. (2008) studied the pattern of the 2000 NW Bohemia swarm in terms of crack and intrusion models. Interesting was the activity burst occurring during 15. October 16–20 h UTC (day 289 in 2000), which could not be explained by a magmatic intrusion. Dahm et al. (2008) postulated that the short and strong subswarm, which was vertically elongated (about 1 km) and migrated clearly upward, was buoyancy driven. They estimated a minimal density difference between the rock and fluid of >1000 kg m⁻³, which would clearly point to a gaseous (possibly CO₂) phase. The low-velocity ratio anomaly of γ ≈ 1.38 falls in the time period of the activity burst from day 289, which gives an independent indication that the source volume was enriched with a gaseous phase at this time.

However, Dahm et al. (2008) also found that the bulk of the 2000 intrusion must have been a dense fluid and cannot be represented by CO₂. This is in accord with the observation that the α₁/β₁ was almost normal during the remaining periods of the 2000 swarm.

Reinhardt (2007) studied in detail the variability of moment tensors and the local stress changes for the 1997 swarm. The analysis indicated that the swarm earthquakes re-activate cross cutting and subparallel en-echelon fault segments, and that the stress field is locally close to fault zones very heterogeneous. This would support the idea that the rocks in the source region are highly fractured and that fluids possibly use these regions and pre-existing fault zones to migrate up- or downward, thereby reducing the strength of these faults.

Few studies tried to answer the question concerning the overpressure in the fault zone. Based on the analysis of moment tensors and fault plane traction Vavrycuk (2002) estimated the fluid pressure to be less than the lithostatic stress only by about 5 MPa for the 1997 NW-Bohemia earthquake swarm. Hainzl et al. (2012) analysed the stress triggering and hypocentre spreading of the 2000 and 2008 swarms to estimate the fluid overpressure responsible for driving the swarm activity. Using two different methods, they arrived at a similar estimate of the overpressure up to 30 MPa. Despite the expected low resolution of these estimates, the two methods agree in the order of fluid overpressure of tens of MPa.

The intrusion model of the 2000 earthquake swarm by Dahm et al. (2008) postulated that a vertical magma dike at slightly deeper levels than the earthquake swarm exists and is cooling and decompressing and thus experiencing slow degassing and bubble formation. The degassing may trigger the onset of the earthquake swarms by creating a pathway to the nearby Pocatky-Plesna fault zone. Bubble formation of magma at a depth of 9–10 km is known to exist (e.g. Caracausi et al.2003). The model would predict that the overpressure is large at the beginning of the intrusion and earthquake activity and becoming smaller when the new fracture volume is formed. This model is strengthened by the observed strong velocity ratio anomaly at the beginning of activity phases in all three studied swarms. Water-magma contact may play an additional role for the abrupt compression of parts of the crust and the triggering of patch-like intrusion growth and velocity ratio anomalies. The free water may use pre-existing zones of weakness to migrate upwards. The violent and vertical migrating shoot-out phase from 2000 October 15 could be a candidate for such a degassing event. The proposed simple velocity ratio monitoring may in future help to possibly detect such types of events.

It should be noted that the low v^P/v^S velocity ratios for the NW-Bohemia swarms were already reported by Vavrycuk (2011) who analysed the moment tensors of tensile earthquakes occurred during the 1997 and 2008 swarms. He employed the linear relation between the ISO and CLVD moment tensor components for tensile earthquakes, where the slope of the linear relation depends on v^P/v^S. Both swarms showed small v^P/v^S in the focal volume; γ₁ = 1.45 for the 1997 swarm and 1.35 for the 2008 swarm, which is quite consistent with our results.

5 CONCLUSION

We used a simple and robust method to estimate the average velocity ratio in the crust beneath a seismic network and in the small source volume of an earthquake swarm. The method only requires measured arrival times of P and S phases, and no earthquake hypocentres and/or origin times. It is suited for automatic applications and monitoring purposes, for instance in regions of earthquake swarms or at volcanos.

The application to three mid-crustal earthquake swarms in NW-Bohemia in 1997, 2000 and 2008 reveals interesting and unexpected results. The P- to S-wave velocity ratio in the source region is highly reduced and variable during the swarms. It partly reaches values as low as about 1.33. We tested the robustness of the anomalies and are convinced that the source volume of the earthquake swarms experiences strong temporal variations, leading to strongly decreased, apparent (effective) v^P/v^S ratios. Typically, the velocity ratio is reduced in the beginning of seismic activity phases. This may be explained by overpressurized gas in the beginning phase of a swarm or subswarm.

Data used in this study were provided by the WEBNET and Sachsen seismic network. Plots were created by GMT plotting tools. We thank Peter Shearer, Joachim Wassermann and Horst Langer for constructive remarks that helped to improve the manuscript. The work of TF was supported by the Czech research projects MSM 0021620855 and P210/12/2451.

REFERENCES