Statistics of the Seismic Sequence and Rupture Directivity of the M5.5 Earthquake in Orkney, South Africa

Despite our general knowledge of earthquake processes, it is still not fully understood how earthquake ruptures nucleate and propagate and why they stop. Also, the controlling factors of the frequency and of the size of earthquakes are subject of ongoing research. We aim to address these questions with a comprehensive study of seismicity in deep South African gold mines. We find here the unique situation that the seismicity consists of both induced earthquakes and aftershocks triggered by the M5.5 Orkney earthquake which occurred in August 2014. We separate the cataloged seismicity and group the events into 5 three classes: the aftershock sequence, seismicity induced by fluids, and seismicity induced by mining activities. We examine statistical properties of earthquakes in each of the three classes. We conclude that the magnitude statistics of both aftershocks and induced earthquakes are influenced by the finite size and geometry of the rock volume of stress perturbation resulting in an absence of larger magnitude events. The magnitude frequency distributions obey the Lower Bound model of magnitude probability. The statistics of dynamic stress drop of aftershocks and induced earthquakes satisfy log-normal distributions but 10 the value range is different, it means that aftershocks are generally characterized by higher stress drops. Another key aspect in our study is the imaging of the propagating rupture of the M5.5 earthquake. We apply the back projection imaging approach using seismological data from two different local networks and retrieve similar results. We conclude that the rupture of the M5.5 earthquake propagated predominantly unilaterally, nearly from North to South over a distance of about 6km. The hypocenters of the aftershock sequence are situated unilaterally in respect to the hypocenter of the main shock and are aligned to the South 15 which confirms the obtained rupture propagation image and the directivity of the main shock.

(2) mining induced seismicity with about 7,000 earthquakes, and (3) fluid induced seismicity with nearly 48,000 earthquakes. Figure 2 shows the hypocenter locations of the three classes of earthquakes.
2 Magnitude statistics of the M5.5 aftershock sequence First, we study the statistics of magnitudes of the aftershock sequence which we educe from the whole catalog. The sequence contains approximately 15,000 events and their hypocenter locations are shown in Figure 3. The area of aftershock activity extends over 5km in NNW-SSE direction, 0.5km perpendicular to the strike direction and 1.5km in vertical direction dipping to the South. The M5.5 main shock is situated at the northern edge of the aftershock cloud. The locations of aftershocks with respect to the main shock location raises the speculation that the rupture of the main shock propagated unilaterally to the South. The magnitude frequency distribution of the aftershock sequence illustrated in Figure 4(a). We estimate the b-value of the Gutenberg-Richter law (Gutenberg and Richter, 1944) using the maximum likelihood method introduced by Utsu (1966) and later modified by Bender (1983). The corresponding graph atop the observed magnitude distribution reveales that the 70 number of aftershocks with a magnitude M > 1.5 deviates from the strict power law of the Gutenberg-Richter relation. Such a phenomenon is often observed for the magnitude frequency distribution of fluid induced seismicity where the seismically active rock volume is restricted to the domain of pore pressure perturbation (Shapiro et al., 2011(Shapiro et al., , 2013Shapiro, 2015). A finite perturbed rock volume influences the magnitude statistics of fluid induced earthquakes. Here, we find similarly for the M5.5 aftershock sequence. The occurrence of aftershocks is restricted in space to the domain of stress perturbation caused 75 by the preceding main shock. This leads to a finite size and to a nearly fixed geometry of the seismically active rock volume which influences correspondingly the magnitude statistics of the aftershock sequence. Therefore, we apply the so-called Lower Bound model of the magnitude probability (Shapiro et al., 2013), which captures the effect of a finite perturbed rock volume on the magnitude statistics, to the Orkney aftershock sequence. We conclude that the lower bound model represents a better approximation to the observed magnitude frequency distribution as shown in Figure 4(b)). Furthermore, with the application of 80 the Lower Bound model we retrieve a geometry-corrected estimate of the b-value of the magnitude frequencies of aftershocks which we report in Table 1. Note, that the difference in b-value is marginal only but it influences seismic hazard assessment. The Lower Bound model yields also an estimate of the maximum possible magnitude of an earthquake occurring within the spatially confined rock volume of a stress perturbation. In the present case of the aftershock sequence of the Orkney earthquake, the maximum possible magnitude of an aftershock is estimated to be M max = 2.4. This agrees with the observation, however, it 85 contradicts Båth law (Båth and Duda, 1964;Båth, 1981) which empirically states that the magnitude of the largest aftershock is around one order lower than the nagnitude of the preceding main shock. Additionally, the lower bound model provides an alternative way to derive the stress drop ∆σ (Shapiro et al., 2013). This parameter is an averaged representation for the aftershock sequence and it is also reported in Table 1. Now, we compare the statistical characteristics of magnitudes of the aftershock sequence to other seismicity observed in the 90 gold mines in the Orkney area. These are two different types of earthquakes caused by human operations, either induced by fluid or induced by mining, respectively. Figure 5 depicts the corresponding magnitude frequency distributions as well as the Gutenberg-Richter law using maximum likelihood b-value estimates. The magnitude frequency distributions are influenced in such way that larger magnitudes are under-represented in comparison to the power-law of the Gutenberg-Richter relation. It reflects again the influence of a finite perturbed rock volume on the magnitude statistics. We find that this influence is more 95 pronounced for mining induced seismicity. We apply the Lower Bound model of magnnitude probabilities to the magnitude frequency distributions of induced earthquakes ( Figure 6) which result in geometry-corrected estimates of the Gutenberg-Richter b-value and the maximum possible magnitude, and in estimates of averaged stress drop ∆σ (Table 1). Fluid induced and mining induced events differ clearly regarding the statistics of magnitudes. We obtain with the Lower Bound model the highest b-value and lowest stress drop estimates for fluid induced seismicity. In contrast, mining induced seismicity is 100 characterized by the lowest b-value resulting therefore in a higher estimate of the maximum possible magnitude.
We also compare the statistics of stress drops for the three different classes of seismicity ( Figure 7). We find it remarkable that stress drops are similar distributed, precisely log-normal distributions, independent of the event type. However, the value ranges are different, it means the stress drop of aftershocks is approximately two orders of magnitude higher than the stress drop of induced earthquakes. We observe the lowest stress drop values for the seismicity that is induced by fluid. Finally, we 105 analyze how the cataloged earthquake stress drops relate to the averaged stress drop estimate resulting from the Lower Bound model. The median values of static and dynamic stress drops are summarized in Table 1 and we find that the Lower Bound derived stress drops are reasonable values and are in good agreement with the median of dynamic stress drops. This finding holds for all three classes of seismicity. (Pty) Ltd (IMS) which is specialized in seismic monitoring of industrial underground activities (Mendecki, 1996). Highquality waveform data with good signal-to-noise ratios were recorded at 4.5Hz 3-component geophones with a sampling rate of 6000Hz. Figure 8 shows the spatial distribution of stations of the three in-mine networks and the cut-out P-waveforms which we use in the MRPI method. To invert for the hypocenter location and to compute the traveltime-grid, we apply a homogeneous 1-layer velocity model with v P = 5.89km/s and v S = 3.69km/s for P-wave and S-wave,respectively, that was provided by Comparing the rupture track with the spatial distribution of aftershocks, we find that the elongation of the aftershock cloud in N-S direction fits well the rupture length. But, the aftershock cloud is shifted southwards by one to two kilometers which we identify by comparing Figures 9 and 10(a). Also, the strike directions of both the cloud and the rupture track are slightly 140 different. One reason for this finding could be the spatial distribution of stations of the in-mine network which is suboptimal because of the limited azimuthal coverage with respect to the earthquake location (Figure 8(a)). However, Manzunzu et al.
(2017) determined focal mechanisms for the M5.5 earthquake and a selected number of aftershocks (Figure 10(b)). They found that the main shock occurred on a steeply dipping strike-slip fault with a strike direction of 182.6 • N. The strike direction of the rupture track coincides in a first approximation with the reported focal mechanism.

145
Additionally, we also apply the back projection imaging using wavefields recorded at the surface network. Figure 11 shows the spatial distribution of stations of this network which is operated by the Council for Geoscience South Africa (CGS) and cut-out recordings of P-waveforms that we use for the MRPI. The azimuthal coverage of the station distribution is slightly enhanced in comparison to the in-mine network but still not perfect. The network is equiped mainly with 3-components broadband seismometers, data are recorded with 200Hz sample rate. Since we consider here the surface network, we apply a more 150 advanced but still rather simple velocity model which is listed in Table 2 behavior indicating an initial and a main rupture phase was also reported by Ellsworth (pers. communication). He carried out a slip inversion for the M5.5 earthquake on a vertical plane and he found that the initial rupture origin is located to the South of the main rupture origin which propagated mainly southward.
The combination of three sources of earthquake data − waveform recordings from a surface network, from a deep in-mine 160 network, and a comprehensive catalog with about 70,000 earthquakes − allowed for detailed investigations of seismogenic processes in a spatially confined, seismically active volume. We separated the cataloged seismicity and grouped the events into three different classes: (1) the aftershock sequence of the Orkney M5.5 earthquake (≈ 15,000 EQ), (2) seismicity induced by fluids (≈ 48,000 EQ), and (3) seismicity induced by mining activity (≈ 7,000 EQ). We studied statistical properties of earthquakes in each of the three classes and compared jointly the results. We focused on the magnitude statistics and on the statistics 165 of dynamic stress drop. Another key aspect in our study was the imaging of the propagating rupture of the M5.5 earthquake. We applied successfully the back-projection method using seismological recordings from two different, independently operating networks and retrieved comparable results.
We concluded that the rupture of the M5.5 earthquake propagated predominantly unilaterally, nearly from North to South over a distance of about 6 km. In general, the coherency of waveforms is high for seismic networks located at teleseismic distances 170 to the earthquake source which facilitates application of the back-projection imaging. In our case, however, the waveforms were recorded at spatially suboptimally distributed stations of local monitoring system. Hence the waveforms can be quite variable due to earthquake radiation pattern, higher frequencies and local heterogeneity of seismic velocities and attenuation structure. To avoid this complexity, there exists a different and simpler geometrical method, the P-wave polarization stacking  Two nucleation slip patches, an initial and a main rupture origin, appeared to locate on the rupture surface. They were seperated horizontally by around 2 km in which the initial origin is south of the main origin. An assimiable observation supporting our findings was reported by Moyer et al. (2017). They identified an immediate foreshock precursory to the main shock which was located in about 1.4 km horizontal distance south of the main shock. Integrating both observations supports the conclusion that the interaction of two patches led to the M5.5 Orkney earthquake. Such a phenomenon is well-studied for subduction zone 185 earthquakes. It was discussed and concluded that interactions and synchronization between multiple patches along the subduction interface led to the M w 9.0 Tohoku-oki earthquake in March 2011 (see Aochi and Ide, 2011;Hori and Miyazaki, 2011;Kiser and Ishii, 2012, e.g.). A similar complex triggering model was reported by Kiser and Ishii (2011) also for the M w 8.8 Maule earthquake which occurred in 2010.
The hypocenters of the aftershock sequence of the M5.5 earthquake are situated unilaterally in respect to the hypocenter of the 190 main shock and are aligned to the South. In general, this observation confirms the obtained rupture propagation image and the directivity of the main shock. Furthermore, the extent of the aftershock cloud in N-S direction is in agreement with the length of the imaged rupture track even though the two are offset slightly. We found that the magnitude statistics of both aftershocks and induced earthquakes are affected by the finite size and geometry of the rock volume of stress perturbation which inhibits the occurrence of larger magnitude events. This influence is most pronounced for the magnitudes of aftershocks whereas it with the Lower Bound model, retain the ranking of the three different classes of seismicity. Interestingly, the Lower Bound stress drops are closer to the median of dynamic stress drops than to the median of static stress drops. Overall, our findings and conclusions have basically the potential to contribute to an improved assessment and a possible mitigation of seismic hazard 205 in active mining districts. In particular, we think that the Lower Bound model combined either with the Omori-Utsu model for aftershock productivity or with a Poissonian model for induced earthquakes would be a first step towards more realistic hazard estimations of mining operations.
Code and data availability. Software codes in order to process and analyze the seismological data were written in MATLAB©®and are available upon request from the corresponding author. The seismological data from the CGS network are downloadable from the CGS 210 webpage and are also available from the supplementary material of Moyer et al. (2017). The in-mine networks data used in this research (catalog as well as stations recordings) are not publicly available. The data were provided by AngloGold Ashanti Ltd. and Harmony Gold Ltd. based on a collaboration within the framw of the DSeis drilling project.
Author contributions. CD carried out the data analysis and wrote the manuscript. JF provided the imaging tools used in this study, and JF and JK assisted substantially in the interpretation of the obtained results. SS and HO are PI of the ICDP funded project DSeis − Drilling into