Spatial Distribution of Triggered Earthquakes in the Conditions of Mining-Induced Seismicity

Abstract—The spatial distribution of the triggered seismic events in mining conditions in the tectonically loaded rock masses is studied using the example of seismicity in the Khibiny Mountains. It is shown that the distribution of distances from the triggering to triggered events, on average, obeys the power-law with a parameter independent of the magnitude of the triggering event. The model of the maximum distances from a triggering event’s hypocenter to the triggered shocks expected with a given probability is derived. It is shown that the model is consistent with the real data. Based on the error diagram analysis, the guidelines are proposed for the practical use of the model.


INTRODUCTION
This work continues our research into spatiotemporal patterns of seismicity in mining regions. Conducting our investigation into the subject, we validated the productivity law established in our previous works Shebalin et al., 2020) for the conditions of mining-induced seismicity by showing that the number of the triggered events shocks initiated by an earlier event (productivity) obeys exponential distribution . Here, the single parameter of the exponential distribution is independent of the magnitude of the triggering event.
In this work, by analyzing the seismicity of the Khibiny massif, we show that the distances from the triggered earthquakes to their triggering seismic events obey the power-law distribution. This conclusion is consistent with the well-known results previously obtained for the aftershocks of the tectonic earthquakes in California and Japan (Huc and Main, 2003;Felzer and Brodsky, 2006;Richards-Dinger et al., 2010). The last two cited works reflect the discussion on whether the dynamic stress transfer can induce aftershocks. Felzer and Brodsky (2006) concluded that the observed power-law distribution of distances between aftershocks and their main shocks agrees with the fact that the probability of the occurrence of the aftershocks is practically proportional to the amplitude of seismic waves. In the cited work it is also shown that this distribution is poorly consistent with the rate-state models describing movements along a fault with friction that depends on the changes in static stress, velocity, and state (Dieterich, 1994;Scholz, 1998). Considering this and the fact that static stress changes at the more distant aftershocks are negligible, Felzer and Brodsky (2006) hypothesized that aftershocks can be initiated by the dynamic stress transfer from the main shock. Subsequently, Richards-Dinger et al. (2010) used the algorithm of (Felzer and Brodsky, 2006) to identify the main shocks and aftershocks and have shown that the power-law decay of the density of the repeated shocks with distance is also observed for the aftershocks that occurred before the arrival of seismic wave from the main shock, which violates causality in the case of the dynamic stress transfer. Thus, the power-law decay of aftershock density with distance does not signify that dynamic stress transfer causes repeated shocks.
In our opinion, there are also no grounds to believe that the power-law character of spatial distribution of the aftershocks indicates their initiation by the dynamic stress transfer from the main shock because the same distribution is also observed for the distances between earthquake pairs (no matter whether they are main shocks or aftershocks) on the global and regional levels (e.g., (Kagan and Knopoff, 1980;Kagan, 2007 and references therein)) reflecting the fractal geometry of the seismicity. We note that a fractal structure was also obtained for crack distribution observed in the laboratory experiments on fracturing of Oshima granite (Hirata et al., 1987).
This study is relevant for the matter at hand as it confirms that the power-law distribution of distances from main shocks to their aftershocks established for tectonic seismicity with M ≥ 2 is also valid for the weak mining-induced seismicity (0 ≤ M ≤ 3.3, 10 4 ≤ E ≤ 8.7 × 10 9 J). This indicates that the power-law character of the spatial distribution of the repeated shocks appears to be universal. At the same time, in order to accept the power-law distribution as valid on all energy scales, it is necessary to carry out laboratory studies similar to those described in (Hirata et al., 1987;Smirnov et al., 2019;.
Mineral extraction in the tectonically loaded rock masses causes manmade seismicity (e.g., (Adushkin, 2013;2016;Kozyrev et al., 2018;Adushkin et al., 2020)). In this case, rock pressure in the underground workings of the operating mines disrupts the continuity of the rock mass including its near-contour part which manifests itself in the dynamic forms as peelingoff and shooting of rocks, dynamic culling, microimpacts and rock bursts and manmade earthquakes (Kozyrev et al., 2016). Just as the tectonic (natural) seismicity, mining-induced earthquakes can trigger repeated shocks (aftershocks) (Plenkers et al., 2010;Woodward and Wesseloo, 2015;Kozyrev et al., 2018;Baranov et al., 2019a;. After a mining-induced earthquake, a decision should be rapidly made as to suspending the work, evacuating the people, and removing the equipment from the danger zone. Thus, the research addressing spatiotemporal patterns of the postseismic processes in mining regions has a clear practical relevance.
As a practical application of the previously established productivity law for mining-induced seismicity and the power-law spatial distribution of aftershocks revealed in this study, we analytically developed a model that allows for estimating the size of the zone where the triggered events are expected to occur with a given probability. We note that this result is important for safe mining operations.

INITIAL DATA AND IDENTIFICATION OF INDUCED EVENTS
Just as in , in this study we use the catalog of seismic events recorded by the seismic monitoring network of the Kirovsk Branch of AO Apatit (Korchak et al., 2014) from 1996 to August 2020 ( Fig. 1). The network currently incorporates 50 threecomponent seismic sensors installed in the Kirovsk and Rasvumchorr mines which record the input signals at a sampling rate of 1000 Hz. The network allows the locations of hypocenters of seismic events with energy E ≥ 10 4 J to be determined as accurately as up to 25 m in the region of reliable recording. The hypocenter locations of the lower-energy events are estimated less accurately, e.g., the hypocenter location accuracy for the events with E = 10 3 J is up to 100 m in the region of reliable recording and up to 25 m in the high-accuracy region.
Processing of data of the KB AO Apatit seismic monitoring network includes calculation of event energy E (in J). In this paper, energy was converted into magnitude based on the formula of T.G. Rautian (1960): logE(J) = 1.8M + 4.0.
Since 1996, complete recording of seismic events by the network has been provided starting from energy of E c = 10 4 J which corresponds to the magnitude of completeness M c = 0. The catalog used in our study contains 71883 seismic events with 0 ≤ M ≤ 3.3. The completeness of the catalog and the accurate, up to 25 m, hypocenter location allow study of very weak seismicity which fills the gap between the laboratory experiments and the field observations. This is an additional test for the universal character of the regularities revealed by both the laboratory studies and the analysis of global and regional catalogs of tectonic earthquakes.
The triggering and triggered events were identified by the nearest neighbor method (Zaliapin and Ben-Zion, 2016) based on the use of the proximity function in the space-time-magnitude region (Baiesi and Paczuski, 2004). This function depends on the parameters of seismicity: the Gutenberg-Richter (GR) b-value and fractal dimension of earthquake hypocenters d f . The idea of the method is that for each event in the catalog (except for the first event), we find its "ancestor" determined by the minimum of the proximity function calculated over all the previous events. If the minimum of the proximity function is below a certain threshold η 0 , the ancestor is assumed a trigger of the analyzed event. Otherwise, the connection between these events is rejected. Here we only consider the highest hierarchy level where the triggering event and its initiated shocks constitute one series. If an initiated shock is itself a trigger, it forms another series. Events that have no triggers are considered background events irrespective of whether they initiate repeated. The η 0 value can be selected by various methods, e.g., (Bayliss et al., 2019;Baranov and Shebalin, 2019;Shebalin et al., 2020). Here, we use the model-independent method  which is preferable in the conditions of manmade seismicity .
The use of the nearest neighbor method to analyze the seismicity of the Khibiny natural/manmade system (NMS) was discussed in our previous paper  (Baranov et al., 2019a;. In the cited work, we obtained the following parameter estimates: b = 1.25, d f = 1.55, log η 0 = -6.25.

DISTRIBUTION OF DISTANCES FROM THE TRIGGERING EVENTS TO THE TRIGGERED EVENTS
For the triggering seismic events with M m ≥ 1.5, we construct the distribution of the distances between these events and the triggered shocks with magnitude M ≥ M m -1.5 initiated by them. According to (Huc and Main, 2003;Felzer and Brodsky, 2006;Richards-Dinger et al., 2010), distances r between main shocks and their aftershocks starting from certain value r 0 obey the power-law distribution (1) with density (2) Here, n is the distribution parameter characterizing the slope of the graph plotted on a loglog scale.
It was established that in the case of the seismicity of the Khibiny massif, the distances from the triggering events' epicenters to triggered shocks starting from r 0 = 0.13 km also obey the power-law distribution ( Fig. 2) with parameter n = 2.28 in the different ranges of trigger event magnitudes M m . The standard errors σ (for parameter n) and the r 0 values are presented in the caption of Fig. 2, and the characteristics of the series are presented in Table 1. The estimation was carried out by the maximum likelihood method (Clauset et al., 2009). Moreover, just as in (Felzer and Brodsky, 2006), parameter n is independent of the main shock magnitude.
The similar result is also valid for vertical (along the depth) distances between the triggering and triggered events (Fig. 3). In this case, we denote the distance in formulas (1) and (2) starting from which the distribution satisfies the power-law by h 0 . The values of parameter n, standard errors σ, and distance h 0 are indicated in the caption of Fig. 3, and the characteristics of the series are presented in Table 1. In the case of vertical (along depth) distances, the scatter in the values of parameter n for different magnitudes is larger than for epicentral distances. This is due to the fact that errors in depth determination are larger than epicenter location errors. Another probable cause is the presence of stress field inhomogeneities along depth increasing the probability that rock pressure manifest itself in the dynamic form (Kozyrev et al., 2019). In turn, these manifestations may lead to the variations in the decay rate of density of the initiated shocks (parameter n). In any case, the ±3σ intervals for the values of n overlap indicating insignificant deviations in the values of this parameter.

DISTRIBUTION MODEL FOR THE REGION OF TRIGGERED SHOCKS
Given that the power-law distribution parameter n is practically independent of the magnitude of the triggering event, the radius R of the circle encompassing the expected initiated events with magnitude M ≥ M m -ΔM is determined by the number of shocks of a given magnitude initiated by the triggering event (trigger productivity).
A triggering event can trigger several dependent shocks constituting a series. As we consider only one hierarchical level, we may suppose that the events in the series are independent of each other. We assume that for each series, the number of the initiated events with magnitudes M ≥ M m -ΔM obey the Poisson distribution with the mean Λ (Zoller et al., 2013). In this case, the probability that all k initiated shocks occur closer than at distance x from the trigger is F r (x) k where F r (x) is distribution described by formula (1). Using the total probability formula, we obtain the distribution of the maximum epicentral distance R max from the triggering event to the most distant aftershock in the series: (3) According to the earthquake productivity law  which was validated for the seismicity of the Khibiny massif , the number of the triggered events obeys the exponential distribution with density (4) Here, the average number of the triggered events is the estimate of parameter L. N s is the number of series triggered by triggering events with magnitude M m ; N is the number of triggered events in series; r 0 , km, is the distance starting from which the distribution of epicentral distances from triggering events to triggered events obeys the power-law distribution (1); N(r < r 0 ) is the number of the triggered events with epicentral distances to their triggering events shorter than r 0 ; h 0 is the distance starting from which the vertical (along the depth) distances from the triggering events to the triggered events obeys the power-law distribution (1); N(h < h 0 ) is the number of the initiated shocks with vertical (along the depth) distances to their triggering events shorter than h 0 . For finding the distribution of distances from the triggered shocks to their triggering events over the set of the series, we combine (3) and (4) at x ≥ r 0 . Thus, we obtain the distribution function (5) and density (6) where F r is the power-law distribution function (1) with density f r (2).
Given that the vertical (along the depth) distances from the triggered shocks to their triggering events also obey the power-law distribution (Fig. 3), the similar relations are also valid for the maximum distances H max along the depth.
Formulas (5) and (6) define the model of the distribution of maximum distances where the repeated shocks are expected. The correspondence of this model to the seismicity data for the Khibiny massif is shown in Fig. 4 (the value L ≈ 3 is specified according to ; this estimate can also be obtained from Table 1 by calculating the ratio of column N to column N s values).

PRACTICAL ASPECTS
In this section, we consider practical aspects of using the averaged model of the distribution of maximum distances. The region where the triggered events with magnitudes M ≥ M m -ΔM are expected to occur with a given probability can barely be estimated directly from distribution (5) because the power-law decay is only satisfied starting from a certain, albeit small distance from the triggering event, at which about half of all the initiated shocks take place (Table 1). In order to allow for this feature and for the limited spatial extent of the mining activity region, we used the Molchan error diagram (Molchan, 2010) visualizing the dependence of the fraction of the missed target events (the rate of failures-to-predict) ν versus the fraction of space-time alarms τ.
When estimating epicentral distances, we assume the space Ω containing 100% of the triggered events to be a circle with a radius of 2.5 km and a center at the epicenter of the triggering event. This Ω corresponds to the zone controlled by the joint Kirovsk mine. For the Rasvumchorr mine, we assume the same space Ω. We estimate the epicentral alarm zone encompassing the expected triggered events as the circle with the center at the epicenter of the triggering event and with radius R q calculated for probability q from the inverse function for distribution (5): . (Parameters of distribution (5) are n = 2.28, r 0 = 0.134 km (Fig. 2a), L = 3 ). We denote this zone by G q and its area by S q . Then, the fraction of the alarm space τ is defined by the ratio of S q to the area S Ω of space Ω, i.e., τ = S q /S Ω . The rate of missed events ν is the fraction of the triggered events beyond the alarm region G q .
When estimating distances along the depth (vertical distances), we assume the space Ω in the form of a segment of length H Ω = 1 km centered at the main shock hypocenter, which corresponds to the vertical (along the depth) zone controlled by the mines and contains 100% of the aftershocks. Then, as the alarm zone along the depth accommodating the expected triggered events, we use the vertical segment V q with the center at the hypocenter of the triggering event and with length H q calculated for probability q from the inverse function of distribution F a (5) for depths with parameters n = 2.29 and h 0 = 0.06 km (Fig. 3a) and h 0 = 3 . In this case, the fraction of the alarm space is τ = H q /H Ω . The fraction ν of the missed target events is the fraction of the repeated shocks that fall beyond the segment V q .
This shape of the zone corresponds to a cylinder with radius R q , height H q , and center at the epicenter of the triggering event. This shape of the zone allows the radius and height of the cylinder to be determined independently, subject to the importance of the forecast.
The dependence of ν on τ for different q is the error trajectory. Diagonal (0; 1) (1; 0) corresponds to a random forecast. The stronger deviation of the error trajectory from the diagonal means the better forecast performance. Parameter q specifies the size of the alarm zone: the larger q means the larger area of G q or V q . The error diagram constructed for different q values based on the retrospective forecast of the region of the repeated shocks with M ≥ M m -1.5 (Fig. 5) reflects the trade-off between two kinds of errors: the increase in q reduces the probability of a missed event but increases the alarm zone and vice versa. Thus, the scalar parameter q can be characterized as an alarm function (Zechar and Jordan, 2008;Shebalin et al., 2014).
The q value is chosen subject to the forecast objective. In some cases, it is important that the probability of a second-kind error, i.e., missing an event, is low. The situation when a strong aftershock can cause undesirable consequences is the example. In other cases, it may be necessary to minimize the area where the repeated shocks are expected in order to reduce the cost of continuing the warning status. To formalize the choice of parameter q, we proposed a method of three strategies (Baranov and Shebalin, 2017). The idea of the method is to determine the limiting points on the error trajectory that correspond to the neutral, soft and hard strategies.
The point corresponding to the neutral strategy (point 0 in Fig. 5) is determined from the minimum of the loss function γ = ν + τ which is the sum of the errors of two kinds. This strategy is used when the cost of the two kinds of errors is approximately the same or not known. The point corresponding to the soft strategy (point 1 in Fig. 5) is determined by the position of the tangent line to the error trajectory, at which, because of the trajectory's closeness to vertical, even a small change in the size of the alarm zone will result, through the decrease in q, in a strongly increased probability of missing an event. Finally, the hard strategy corresponds to the point (2 in Fig. 5) at which the tangent line to the error trajectory is such that the increase in the alarm region will not reduce the fraction of the missed events because of the closeness of the trajectory to horizontal. The q, ν, τ, R max , and H max values corresponding to the neutral, soft, and hard strategies are presented in Table 2.
Model (5), (6) constructed over the set of the series can be used as a first approximation of the region where the aftershocks triggered by a seismic event with M ≥ 1.5 are expected immediately after its occurrence. The independence of the estimates of the epicentral distance and the estimates along the depth allows us to use different strategies to select the radius and the height of the cylindrical region subject to the location of the main shock. These estimates for a specific series can be improved by taking into account the first aftershocks. The constructed model in this case can be used as a basic one for testing the models that employ information about the first aftershocks. An example of the use of a basic model for estimating the magnitude of the strongest aftershock and the duration of the hazardous period is presented in (Baranov et al., 2019b;Shebalin and Baranov, 2019).

CONCLUSIONS
Based on the seismicity data for the Khibiny massif, it is shown that the distances between the triggering and triggered events, on average, have a power-law distribution with the exponent that practically does not depend on the magnitude of the triggering event. The established regularity is consistent with the previously obtained conclusions for the aftershocks of the tectonic earthquakes (Huc and Main, 2003;Brodsky, 2006;Richards-Dinger et al., 2010). This result has an important theoretical value for statistical seismology as it, firstly, validates power-law distribution for weak seismicity with magnitudes 0 ≤ M ≤ 3.3 and, secondly, gives ground to believe that the regularities revealed for tectonic seismicity are also valid at mining in the tectonically loaded rock masses.
In this study, based on earthquake productivity law, we have constructed the model of the maximum distances of the expected aftershocks which allows the estimates to be obtained starting from the time immediately after the main shock. The consistency of the model with the real data is demonstrated. The guide-  lines for the practical use of this model are substantiated based on the analysis of the error diagram.

ACKNOWLEDGMENTS
We are grateful to the reviewers for their deep and meaningful analysis of the work and for their valuable comments.

FUNDING
The work includes the results obtained in the project no. 19-05-00812 supported by the Russian Foundation for Basic Research and the results obtained in partial fulfillment of state contract from the Ministry of Science and Higher Education of Russian Federation (under project no. 075-01304-20).

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.