Identification of seismic doublets occurred on Rudna mine, Poland

Seismic events with similar focal mechanisms and similarity of Green`s functions exhibit common waveform shapes recorded on the same seismic station. The similarity of seismograms can also be observed in the case of man-induced seismicity on the continuous excavated area. In this work, we focused on strong and potential mining-induced seismic doublets, our study was using concepts of signal analysis, using an in-mine network with quite low dynamic range. Since in such case records for strong events are usually affected by clipping, signal analysis were performed with both raw signal cross-correlation (CC) and binary signal cross-correlation (BCC). We found that for events which fulfilled the established criteria, waveform similarities are significantly higher than for other possible doublets. Using BCC and CC approaches, we noticed crucial influences of double couple nodal planes’ as well as P and T axis orientations on waveform similarities although the focal mechanisms were characterized by very high non-double couple components. Finally analysis confirmed that human activity represented by mining is able to produce strong (M > 3) repeating mining tremors.


Introduction
The seismic signal recorded on a seismological network contains a lot of information needed to understand source physics. Of greatest importance are seismic phases used to solve for hypocenter location, focal mechanism or velocity model. Even a simple view of seismograms can provide some advanced analysis. One of the most interesting signal based studies are those dealing with repeated earthquakes with common shapes of recordings and clustered on particular areas. Doublets or more wider multiplets of events are often observed on active tectonic zones (e.g., Kagan and Jackson 1999), low-seismicity and volcanic areas (e.g., Gregersen et al. 2007;Quintanar et al. 2004) or human induced seismicity (e.g., Gibowicz 2006;Wehling-Benatelli et al. 2013) and provide a possibility to improve seismological analysis. Some conditions have to be fulfilled for multiplets in a sense of waveform similarities (Geller and Mueller 1980): two events should have almost the same location, with separation less than one quarter of the shortest wavelength (k/4-criterion) also both source mechanisms and source time functions (STF) should be very similar. Such conditions might have a direct physical interpretation and the most popular is that similar earthquakes occurred on the same part of fault by rerupturing of the same asperity at a different time. This simple concept provides a lot of possibilities to improve seismological studies, from the new events' location procedures up to new conceptual source models (e.g., Poupinet et al. 1984;Got et al. 1994;Waldhauser and Ellsworth 2000;Anooshehpoor and Brune 2001;Schaff and Waldhauser 2005;Baisch et al. 2008;Rudziński and Dębski 2011;Spottiswoode and Milev 1998;Godano et al. 2015).
In contrast to natural earthquakes, the mining induced events never occur without human activity which produces fast changes of stress distribution within rocks close to excavation panels (Gibowicz and Lasocki 2001). Even though dynamic response of the rock to industrial activity can change the velocity model, typical pattern of similar earthquakes might be observed in different induced environments (Evans et al. 2005). Due to clustering of the induced seismicity, waveform similarities might be treated as a good indicator of processes responsible for rock failure. It is well known that mining-induced seismicity tends to cluster in space and time (Gibowicz and Lasocki 2001). Clustering can be also observed when one is dealing with moment tensors (MT) for a number of close-spaced induced tremors in hard coal mining Sen et al. 2013). Analysis of similar tremors (Wehling-Benatelli et al. 2013) also showed that it might be a useful tool to indicate different types of failures caused by coal seam exploitation.
Multiplets and doublets of seismic events are also not exceptional in Polish copper mines. Based on catalogues studies, Gibowicz (2006) developed essential criteria characterizing multiplets in this particular area. Hypocenters of two events are separated by no more than 200 m, the difference in origin time is not longer than 20 days, and the difference in moment magnitude is not larger than 0.15. The study found that in the Rudna mine-one of the analyzed copper mines, 11 doublets, 4 triplets and 1 quadruplet were recorded in 2004 on different mining panels (Gibowicz 2006). In 2010, Rudna mine started excavation on the new mining panel no. XX/1 (Fig. 1). The panel was not affected by previous work. During excavation eight strong induced events with M [ 3 occurred and gave a possibility to provide an analysis of possible doublets among strong mining tremors. Directly following Gibowicz's criteria, we found and studied possible seismic multiplets using waveforms analysis of records from in-mine network. Even though the network is dense enough to locate sources with epicentral errors less than 100 m, the low dynamic range influences seismograms by clipping. This is a factor which must be considered during signal analysis. In this paper, we try to deal with this problem and finally we summarize results obtained in this case study. Our investigation redefines the criteria for multiplets for clustered events occurring in Polish copper mines.

Site and data description
The Rudna mine is a deep copper mine belongs to the Legnica Głogów Copper District located on southwest Poland. Seismicity in the mine is quite strong with more than 1000 events with M [ 1.5 recorded every year. An ore level is located at the depth from 900 to 1200 m and a copper is, at the same time, excavated from more than five mining panels. The XX/1 was excavated between 2010 and 2012. The productive level of panel XX/1 of Rudna mine is located at the depth of about 1100 m and is overlain by some 60-92 m thick dolomite layer followed upwards by rigid anhydrite strata of around 160 m thickness. Below the deposit is a layer of hard Rotliegendes sandstones of about 300 m thickness. Ore is extracted by means of a room-andpillar exploitation system. From 2010 to 2012, the XX/1 mining panel in the Rudna copper mine was struck with several strong seismic events (Orlecka-Sikora et al. 2014). Based on Gibowicz (2006) criteria, we decided to select eight events (Table 1; Fig. 2) for further analysis of possible doublets existence. Our work was facilitated by the mining catalogue obtained from records of the underground seismic network operated in the Rudna mine. The initial events location provided by mining geophysical surveys were corrected using the location procedure worked with probabilistic inverse approaches, the algorithm already tested in the Rudna mine (Rudziński and Dębski 2011). The final hypocentres location with errors for particular coordinates are included in Table 1. While epicentre coordinates are estimate with high precision, accuracy in depth is much lower. Anyway we believe that all events have a similar origin (mining induced) and occurred on the similar depth.
The seismicity in the Rudna mine is continuously monitored by a local seismic network (Fig. 1). The network belongs to the mine and consists of 32 short-period vertical seismometers Willmore MK-III. Even though the number of stations and 500 Hz sampling rate give an opportunity to locate hypocenters within an accuracy of about 50 m in epicentral distance and 150-200 m in depth (Rudziński and Dębski 2011), less than 70 dB dynamic ranges can cause clipping of seismic signals. This condition is especially obstructive when a very strong event is considered.
For selected events from Table 1, two source parameters, namely moment magnitude M w and source radius, have been determined using spectra of P and S waves using the relationship between spectral and Brune's source parameters. All estimations were performed using the formalism of Andrews (1986), Brune (1970Brune ( , 1971) and Wyss and Brune (1968), described by Gibowicz and Kijko (1994) and Niewiadomski (1997). According to those formalisms, moment magnitude is obtained with seismic moment which is directly related to the low frequency level of the far field displacement spectrum: q-density of the source area, c-P or S wave velocity in source, R-source-receiver distance, X 0 -spectral level, R c -free surface correction, S c -site correction, F c -P or S wave coefficient.
Source radius r 0 was estimated using Madariaga's (1976) circular fault model; in such formalism (Brune 1970), the source radius is inversely proportional to the corner frequency f c of P or S wave: where K c is a constant depending on source model and b 0 is S-wave velocity in the source region. The relevant parts of seismograms of P and S waves were selected manually, and then transformed by Fast Fourier Transformation (FFT). The resulting amplitude spectra were corrected for attenuation effects with Q = 400 and 200, respectively, for P and S waves. For further calculations we set velocities of P and S waves in source at V p = 5900 m/s and V s = 3400 m/s. The station number used for analysis differ from 15 to 25 in case of P-wave and from 4 to 10 stations in case of S-wave according to the unclipped recordings. The source radii R of the studied events were between 160 and 208 m (Table 1), which corresponds to the corner frequency from 1.5 up to 3.5 Hz. The source size estimates are not exceptional when compared with the average event source radius of the events from Rudna mine, which was mainly about 200 m (Lizurek and Wiejacz 2011).
In next step for selected events a focal mechanisms were calculated using the full moment tensor (MT) inversion in time domain. Calculations of the full MT from the records of the underground network of the mine were performed with the use of FOCI software (Kwiatek et al. 2016). The input parameters are the amplitude and polarity information on the first P-wave displacement pulses. According to Fitch et al. (1980), De Natale et al. (1987 and Aki and  X and Y are local coordinates in Cartesian system, Z is a source depth in meters, DX, DY, DZ are location errors in meters. M w is a moment magnitude while R is a foci radius Richards (2002), the recorded displacement for the P-wave phase is: where q is the average medium density, r is the sourcereceiver distance, a is the average velocity of the P wave, M is the seismic moment tensor, l is the P wave direction at the receiver, and c is the P-wave direction at the source. The source time function (STF) was based on the Haskell's source model (Haskell 1953): where T is the rupture time. The MT is obtained by solution of a set of N equations of type (3). The deviatoric, pure shear and full MTs (Table 2) were calculated using the L2 norm as a measure of the misfit (Wiejacz 1992;Awad and Kwiatek 2005). The recorded first arrivals may be of different types: direct P waves or refracted waves from the overlying anhydrite layer or the thick strata of sandstone underlying the ore deposit. This is caused by the geological situation. The layers are almost flat, the dip of the layers being about 4°NE. The ore bearing strata are a dolomite layer of 60-90 m thickness. Above the ore bearing strata, there exists the evaporate series mainly composed of the anhydrite of 160 m thickness, and below the ore bearing strata there is a sandstone layer of about 300 m thickness.

Doublets analysis
In Fig. 2, all doublets are shown with its full MT solutions while Figs. 3, 4, 5 and 6 shows the corresponding fault plane solutions, i.e., the double-couple (DC) contribution (see Table 2). Based on the results of the moment tensor inversion, four doublets are proposed (Figs. 2, 3, 4, 5 and 6).
Figures 3, 4, 5 and 6 show the following pairs, Ev2-Ev4, Ev3-Ev5 and Ev7-Ev8 exhibit significant similarities in nodal planes distributions and Ev1-Ev6 with similar MT solutions. The full MTs ( Fig. 2; Table 2) for Ev1-Ev6, is characterized by almost the same decomposition and the smallest DC component among four events pairs. This observation is not an extraordinary feature, since the nondouble-couple (non-DC) events are quite often observed in induced seismicity (Gibowicz 2009, Rudziński et al. 2016. Comparison between Figs. 6 and 2 show that although the MTs of events 1 and 6 are very similar, the corresponding DC solutions are very different. Both the hypocenter location (inter-event distances) and other source parameters (Tables 1, 2 and 3) support our choice of possible doublets. Final inspections suggest the following four doublets: -Ev1 and Ev6, -Ev2 and Ev4, -Ev3 and Ev5, -Ev7 and Ev8.
Since the highest dominant frequency for P wave is 3.5 Hz and based on local average velocity model the P wave velocity is 5.5 km/s, the shortest distance satisfying the k/4 criterion should be less than about 390 m. Note that all inter-event distances fit into these limits (Table 3). The last parameter characterizing seismic doublets in the Rudna mine is interevent time. Since this parameter is very difficult to set up, in the previous studies has been selected arbitrarily as 20 days. In this work, we assumed that it is significantly less important than the distance criterion. Our selected doublets are characterized with following time differences: Ev2-Ev4: 28 days, Ev3-Ev5: 33 days, Ev7-Ev8: 127 days and Ev1-Ev6: 347 days. Even the shortest time distance exceeded 20 days criterion.    Two groups of doublets can be now considered. The first one consists of three pairs: Ev2-Ev4, Ev3-Ev5 and Ev7-Ev8. The second one contains just Ev1-Ev6 doublet. For selected doublets, we propose the following seismic signals analysis. On the first step let us choose one the most adequate seismic station and roughly compare the seismograms obtained for two events from each pair. The selected sensor should be far enough from clusters centroid to avoid some recordings problems especially seismograms clipping in S-wave coda. To show how the waveforms overlap each other, we used a station located about 6 km away from the XX/1 mining panel (black triangle in Fig. 1). The results are gathered in Figs. 3, 4 and 5 (first group) and in Fig. 6 for the last pair.
The first group of events exhibits very interesting features in both focal mechanisms and waveform similarities. It should also be noted that the rough line connecting two events forming doublet, corresponds to directions of dominant discontinuities in the mining panel (Mining Geophysical Survey personal information 2015). This observation suggests that for this group the DC-part rupture on pre-existing discontinuities. The most important conclusion is that waveforms for doublets observed in this group are almost identical, even though tremors have, among pairs, other magnitude ranges.
Some other conclusions follow from Fig. 6. The pair Ev1-Ev6 in comparison with the previous doublets has different fault plane orientation and the waveforms not overlap each other. On the other hand, the complete MT results plotted on a lower hemisphere are identical between Ev1 and Ev6 (Fig. 2). Moreover, based on miners observations, the line between these two events does not correspond to any existing geological discontinuity. These two features are very interesting, when we take into account almost the same MT decomposition and suggest that the focal mechanisms for both events are rather connected with pillar destruction without contribution dealing with preexisting lines of weakens, what is not exceptional at the Rudna mine (Lizurek and Wiejacz 2011). One can make a statement that Ev1-Ev6 should not be treated as a doublet.
It is worth to calculate the parameter g which measures the degree of fault rupture overlapping (Kagan and Jackson 1999), where L represents the rupture length, considered as double source radii R while D is a distance between events hypocenters. The value of g larger than 1.0 suggests that the rupture zones overlap. Table 4 includes g for all analyzed pairs. All doublets, including Ev1-Ev6, have g larger than 1 strongly suggesting overlapping of its rupture zones. For the first group g coefficients, however, are more than four times greater that for pair Ev1-Ev6. Let us also note that if we consider the distances including the higher location errors in epicentral coordinates, only pair Ev1-Ev6 has a g less than 1. The previous part was mostly based on visual inspection of seismic signals. Now, we are going to investigate waveforms similarities using signal analysis. This approach has many advantages, nevertheless very serious problem, namely signal clipping (Fig. 7) caused by network low  dynamic range is a factor influencing our investigation. To minimize this issue, we used not only normalized signal cross-correlation C(s) (CC) (e.g., Schaff and Waldhauser 2005 and references therein): where X and Y are two discrete time series of length N and s denotes the time lag but also binary cross-correlation (BCC) and frequency cross-correlation (FCC) both eliminate saturation effects. Analysis was done with seismic traces for all possible doublets-stations pairs. During analysis with CC method, to avoid signal clipping problems, correlation measurements were made just for 5 s windows lengths, 1 s before and 4 s after first arrivals.
Next maximal values of correlation coefficients were used for further study. The resulting CC coefficients were gathered into two types of figures. The first one includes histograms, which show a number of stations versus coefficient intervals. The second one describes spatial distributions of waveform similarities in the form of the map, where correlation coefficient intervals are sketched for each station. An overview of the normalized cross-correlation results are given in Figs. 8, 9, 10 and 11. Further all possible event pairs (28 cases without coefficients for pairs consist of the same event) were investigated toward average CC values. The values range between 0.3 and 0.8. Only pairs Ev2-Ev4, Ev3-Ev5, among all studied pairs were characterized by average CC equal or exceeding 0.6. In Fig. 12 (left), average CC coefficients are presented for all events.
As in the previous analysis, different doublets are characterized by different CC features. The pair Ev2-Ev4 exhibits not only high signal overlapping in selected channel but has also the highest average CC. The next two pairs, Ev3-Ev5 and Ev7-Ev8, have average CC [ 0.6, even though for doublet Ev7-Ev8 eight records have very small CC-below 0.2. On the other hand some stations for the pairs Ev7-Ev8 and Ev2-Ev4 reached CC values above 0.9 while for Ev3-Ev5 most signals have an CC 0.7. Let us note that northeastern stations exhibit smaller CC, especially in case of Ev1-Ev6; we expect that it can be, at least partially, caused by some random environmental and cultural seismic noise background. In some way, differences in background noises can be explained by a very large mining waste neutralization repository. The facility is located in the eastern part of the Rudna's area (see Fig. 1). The Ev1-Ev6 pair is characterized by the lowest average CC, although the high similarity of the full MT solution with dominance of compensated linear vector dipole (CLVD) term and almost identical T-axis orientation. The main difference which can influence the signal CC is the nodal plane and P-axis orientation whose strike differ by about 90°in this case. These features of the rupture orientations were crucial for seismic wave radiation pattern of those events, which were obviously different according to the results of average CC of the studied signals. Since there is not a sufficient similarity in waveforms, the pair Ev1-Ev6 is not a doublet.
For CC analysis we used just 5 s signals and for pair Ev7-Ev8 there is not direct evidence that obtained CC values denoted a doublet, we decided to support our double hypothesis by eliminating clipping effects. The technique is based on BCC methodology and in comparison with CC method uses the whole available signal data. This part of the data analysis consists of four steps concerning signals: -Mean and trend removal, -Ambient noise elimination by normalizing all amplitudes larger than a fixed value to zero, -Median filtering -Signal binarization using one-bit normalization method. Figure 13 shows a raw seismograms (top), seismogram after signal processing (middle) and binary signal one-bit  Fig. 12 (middle) as an average BCC coefficient. These results strongly support observation obtained with CC methodology. We noticed that the biggest BCC average values are connected with the proposed doublets including Ev7-Ev8 except Ev1-Ev6. Our results suggest that the BCC approach can be valuable analysis for data affected by clipping. It will be very important issue in case of clipped records within seconds after P wave first arrival. Third method used in our studies was the FCC. In this method authors calculate correlation between histograms of frequencies of signals obtained in previous step.
Frequency vector was calculated as 1/(k i Ádt), where dt is sampling rate, and k i is number of following samples without changing a sign. Histogram was calculated with step 0.05 Hz. Results are gathered in Fig. 12 (right). The results obtained with the FCC are more complicated in comparison with the CC and the BCC solutions. The similarities for pairs Ev2-Ev4 and Ev7-Ev8 are clearly visible. On the other hand the pair Ev3-Ev5 has average coefficient even smaller than pair Ev1-Ev6. Moreover, there are few more possible pairs like Ev1-Ev7 or Ev1-Ev8 with quite high coefficients. Nevertheless, when the results from tree methods are gathered together, we clearly see that doublets Ev2-Ev4, Ev3-Ev5 and Ev7-Ev8 are similar in the sense of seismograms similarities.

Discussion and conclusions
Our study shows that an analysis of waveform similarities can be a robust method to distinguish two similar induced seismic events occurring on mining area even in case of data set recorded on network with low dynamic range. Among all possible combinations of events pairs, just three doublets exhibit signals analysis features to be treated as doublets and all of them belong to the first group. Our expectation is that average similarity threshold should be above CC [ 0.6. The same conclusions can be derived if we used BCC/FCC analysis. Some very interesting features observed within an analysis of particular doublets should be concluded. Apart from the most obvious, namely the influence of events separations and general types of focal mechanisms to waveform similarities, we can observe very strong influence coming from nodal planes orientations resulting from the main stress axes orientation. It is especially interesting since the non-DC part of the full MT solutions within the studied pairs was not influencing the result of average BCC/CC/FCC analysis as much as the nodal plane orientations differences. This observation suggests that non-DC parts in the full MT could be contaminated by noise but on  the other hand the thorough study provided by Orlecka-Sikora et al. (2014) states that source mechanisms occurred on the panel XX/1 are better explained by the full MT (tensile source mechanism) and extended source rather than DC part of point source full MT. Let us note that the difference in strike of the nodal plane was crucial. From four pairs taken under investigation, the only pair with low average BCC/CC/FCC was the one with significantly different nodal plane strikes. Knowing that, we calculate a radiation patterns for all pairs caused by slip on fault. We followed with method described by Kwiatek and Ben-Zion (2013) and in Fig. 14 we present the radiation patterns for Ev1 and Ev6. The pair has small CC, especially on northeastern part of the network (Fig. 11). It can be noticed that the radiation pattern can be another feature responsible for that. The inter event time is a less influencing factor. Even though the doublet with the highest CC occurred within just 28 days, next two doublets have time differences up to 127 days and measured CC higher than 0.6. We can conclude that the P and T axis orientation and, therefore, the nodal planes of DC solution of the full MT plays the most significant role in the similarity analysis. The  more similar nodal planes and P, T axis orientations the more similarity in waveforms we observed. Since the most similar pairs have nodal planes in accord with the main discontinuities in the mining panel, they probably followed the same rupture. This may be an important indicator of the rupture properties within the rocks in the particular mining environment.
The last but not least is that the BCC and FCC signal analysis work very well with clipped seismograms although the results obtained with the BCC much better support proposed doublets. The method should be especially valuable when clipping affects the whole record and CC with raw signal is hardly useful for studies. Results presented here could be treated as a first step to provide automated signals analysis to identify similar events on mining environments.