Mining is typically accompanied by various degrees of seismic activity. This is known as mining-induced seismicity and it is observed both in hard rock and soft rock environments [1–4]. The mining seismicity is typically induced by different operations, such as rock excavation, rock blasting, and others. This is a direct result of a constantly changing underground stress field that leads to the creation of new fracture zones or reactivation of existing fractures/faults [3–5]. Induced seismicity is one of the hazards in mines and it is observed in potash, coal, metal ore, uranium, and diamond mines within a wide range of geological settings [6–8]. Induced seismicity is observed during other anthropogenic operations such as: hydrocarbon production from tight shale reservoirs, enhanced geothermal energy generation, and high-volume underground fluid injection [9–11].

In the context of mining industry, several geohazards are associated with mining operations. These include: water inflow, the occurrence of earthquakes, rockbursts, development of fracture zones, to mention a few [2, 4, 12, 13]. In case of potash mines, the intrusion of undersaturated brine into a mine has the potential to inflict significant damage to a mine internal structure and result in erosion of pillars, weakening of openings, and rupture of the caprock [12–14]. Among those geohazards, the occurrence of earthquakes and /or microseismicity is a critical factor that has to be monitored and studied. The occurrence of mining seismicity is controlled by various aspects of mine operations and settings such as depth, spatial expanse, geological setting, which includes varied rock formations and layering, the presence of fault systems and fractures, as well as the changes in the ambient stress field and pore pressure [15]. In hard rock mines, seismicity is usually associated with an active mining front due to rock excavation [16]. Whereas in soft rock mining, the host rocks may experience plastic creep due to rock extraction. As a result, the seismicity can be observed in the surrounding rocks [17]. Various aspects of the microseismicity has been also extensively studied in rock fracture experiments, for example in Refs. [18–20].

In this study, we focused on the statistical analysis of microseismicity associated with mining operations in a Saskatchewan, Canada, potash mine. Several statistical models have been considered to investigate various aspects of induced seismicity in the mine. Specifically, a comprehensive analysis of the frequency-magnitude statistics of the mining events has been performed. This included the estimation of the magnitude of the completeness of the analyzed catalog by utilizing several catalog-based methods. Also, we modeled the frequency-magnitude statistics of microseismicity by considering several parametric distributions to establish the most appropriate model that was relevant to that mine. In addition, we also considered the standard exponential distribution, the upper-truncated and tapered Pareto distributions, and the mixture of the above distributions. The maximum likelihood method was utilized to estimate the model parameters. In order to classify the mining seismicity in terms of the mode of triggering and possible rheological regimes, the NND clustering method was applied. The influence of the mining operations on the triggering of microseismicity was quantified by using the cross-correlation analysis. The main goal of this work was to establish and quantify the current pattern of seismic activity in the mine in order to monitor and observe any significant deviations in seismic activity during the future mine operation. The studied seismic catalog also presented a unique opportunity to investigate the statistical properties of microseismicity in the soft rock environment.

Potash deposits in Saskatchewan are mined from the Prairie Evaporite formation. This layer contains one of the largest global potash deposits of sylvinite ore and is mainly used for fertilizer production. Saskatchewan potash deposits are unique in the world because they are very extensive, flat, and generally intact, leading to low-cost extraction. Evidence suggests that the age of these strata of sediments is Middle Devonian. The Prairie Evaporite formation, with a thickness of 100–200 m, is located below 400–500 m of Devonian carbonate with thin layers of shales, anhydrite, and salts followed by 400–500 m of Cretaceous shale and sands with glacial tills on top [21–23]. Besides common microseismicity, several moderate earthquakes have been observed in the past in Saskatchewan's potash mines [12, 13]. Those earthquakes are thought to occur in the layers above the mining horizon, mainly in the carbonate Dawson Bay formation. Changes in the local stress regime and ongoing convergence of the underlying strata into the openings are commonly viewed as one of the leading causes of mining-induced seismicity in potash mines in Saskatchewan [24].

To safely operate a mine, it is important to impose reasonable control on the factors that can trigger earthquakes [3, 4, 25, 26]. In this respect, the statistical analysis of past mining seismicity plays a critical role [27, 28]. By analyzing the patterns of seismicity, one can infer more about their possible causes and relate them to specific mining operations or geomechanical processes. As a result, this can help mitigate the risk associated with induced mining seismicity. The standard approach involves the analysis of the frequency-magnitude statistics and spatio-temporal patterns of mining seismicity [27, 29–31]. This is used to improve the re-entry protocols in mines [32, 33], to study aftershocks and clustering [28, 34], to develop the hazard assessment approaches of mining seismicity [2, 31]. In addition, one can also consider more advanced methods to study mining seismicity. These include the analysis of the clustering aspects of seismicity; any possible deviations from Poisson statistics in the occurrence of events; whether earthquake magnitudes are correlated or not. Therefore, the statistical analysis is a crucial step in assessing the associated seismic hazard and is used to mitigate the corresponding risk.

The statistical analysis of seismicity offers various approaches to quantify the occurrence of earthquakes in magnitude and spatio/temporal domains. Among those approaches, the clustering aspects of seismicity can be used to identify possible triggering mechanisms that drive the occurrence of earthquakes. Various clustering identification methods have been developed in the context of the statistical data analysis [35, 36] and several have been used to study the spatio-temporal patterns of natural and mining seismicity [37, 38]. The overview of clustering aspects of microseismicity in Australian and Canadian mines was provided by Hudyma [39]. It was suggested that mining seismicity is primarily induced by the local failure of rock mass due to various factors including: mining-induced stresses, geological structures, and mining operations. That analysis was based on the application of the hierarchical agglomerative clustering methods to identify individual seismic sources and rock mass failure mechanisms [39]. The mean minimum distance, hierarchical, and K-means clustering methods were used in the study of seismicity in a coal mine in Poland [40]. A 3D spatial clustering methodology for short-term mining seismicity was developed [41] by introducing modifications to the density-based clustering of applications with noise method [42]. The expectation-maximization algorithm was used to identify probabilistic kernels associated with active fault segments to classify seismic events based on their seismic source mechanism, location, and its uncertainty [43]. A hierarchical clustering analysis incorporating Ward's minimum variance method was applied to seismicity in a coal mine to quantify the seismic hazard [44]. One of the goals of the clustering analysis is to identify the relationship between seismic sources and the influence of inducing operations in order to assess and quantify the seismic hazard. For example, the effect of declustering of earthquake catalogs on the frequency-magnitude analysis of seismicity was investigated to asses the estimation of the b-value [45].

A separate method to analyze the clustering of seismicity is based on the assumption that earthquakes can be linked by a suitably computed nearest-neighbor distance in a multidimensional space spanned by the time, space, and magnitude domains. It is known as the nearest-neighbor distance (NND) method and it was first proposed by Baiesi and Paczuski [46] and considerably expanded by Zaliapin et al. [47], Zaliapin and Ben-Zion [48, 49]. It was also adapted to perform declustering of earthquake catalogs [50]. The NND method was also applied to analyze the seismicity in the Western Canada Sedimentary Basin due to a significant increase in induced seismicity observed in the last decade [51–53].

The paper is organized as follows: Section 2 presents the background information about the mine and seismicity catalog. In Section 3 the statistical methods to analyze mining seismicity are formulated. The results of the statistical analysis are reported in Section 4. Section 5 discusses the results and conclusions are drawn in Section 6.

In order to analyze the induced microseismicity in a potash mine in Saskatchewan, a manually reviewed microseismic catalog was used. The spatial distribution of microseismic events above magnitude m ≥ −2.97 is shown in Figure 1. The catalog spans the time interval between 01/05/2019 and 07/09/2021. The original catalog contained 11,633 events. For the analysis, it was reduced to 9,268 events by manually removing the events associated with the network calibrating surface shots that were performed between 02/12/2021 and 02/19/2021. The magnitudes of events in the catalog are reported with two decimal points.

To perform the statistical analysis of mining seismicity, we analyze several aspects of the occurrence of events. This was accomplished by performing the analysis of event catalog data in magnitude, temporal, and spatial domains. The following specific methods were employed.

The statistical analysis is crucial to quantify patterns of seismic activity and provides insight into possible physical mechanisms that are responsible for triggering of mining microseismicity. It is standard practice to assume that the occurrence of seismic events can be approximated by a stochastic marked point process in time and space. This approximation assumes that seismic events are governed by certain statistical models in temporal, spatial, and magnitude domains. The approach undertaken in this work was to establish and apply these statistical models to mining microseismicity in order to make specific inferences about the physical and statistical properties of the mine host-rock environment.

When working with a microseismic catalog, one encounters a standard problem of the incompleteness of the catalog below a certain magnitude m_c. The underestimation of this threshold magnitude can introduce a significant bias in the statistical results. Therefore, we performed a comprehensive analysis to establish the magnitude of the completeness for the mining catalog we used in this study. Specifically, we employed the method of maximum curvature (MAXC), the method of b-value stability analysis (MBS), and the method based on the goodness-of-fit test (GFT). The details of the methods used are summarized in Section 3.3. We applied these methods by first assuming that the event magnitude distribution can be described by the left-truncated exponential distribution (4) or a more commonly accepted form of G-R scaling (3). The MBS and GFT methods produced a consistent estimate for the magnitude of the completeness m_c. The MAXC method is known to underestimate m_c, typically, by 0.2 magnitude units. Therefore, we decided to use the value of m_c = −0.94 from the MBS method as it is somewhere in between of the two values returned by the GFT method. This gives that 93.7% of the distribution of magnitudes above m_c = −0.94 can be explained by the G-R fit (Figure 3) and this value of m_c was adopted for the subsequent analysis.

When fitting the distribution of event magnitudes above m ≥ −0.94 with the exponential model (4), we obtained a rather large value for β = 5.5 ± 0.35 or corresponding b = 2.39 ± 0.15. This indicates that the mining microseismicity is distributed over a narrow range of magnitudes above the completeness level and displays swarm-like characteristics. An earthquake swarm is typically characterized by the occurrence of events of similar magnitude sizes. This can be related to the fact that the microseismicity is driven primarily by mining operations and is occurring on newly created fractures with a limited size distribution [68]. This also indicates that the frequency-magnitude statistics in soft rock environment differs from what is typically observed in hard rock settings [28, 69]. The observed higher than typical tectonic b-values were also observed in other instances of induced seimsmicity [70, 71]. b-value can also be indicative of the state of stress in the crust [72–75]. In tectonic settings, high b-values are typically associated with an extensional stress regime resulting in predominantly normal faulting for seismic events [76]. One possible explanation for the obtained high b-value result is the occurrence of events in ductile rock environment as was suggested by [72] for rock fracture experiments.

To investigate this further, we performed a cross-correlation analysis between the microseismicity rate and the corresponding mine-cutting operations in one particular subregion of the mine. The computed cross-correlation coefficient clearly shows a pattern of negative and positive correlations between the two time-series (Figure 7). The pattern disappears when the cutting time-series was randomly reshuffled (Supplementary Figure S13). In addition, most of the events occurred at rather shallow depths in soft rock environments without the presence of any major fault structures. The performed cross-correlation analysis does not fully explain what kind of mechanisms are responsible for triggering of the microseismicity. The main limitation is that the cutting rate is only a proxy for the actual mining operations. In addition, it may provide an indirect indication that some delayed triggering mechanism is also present influenced by the ductile deformation of soft rock material.

The left-truncated exponential distribution (4) or its equivalent, G-R scaling (3), is the most commonly used model to describe the frequency-magnitude statistics of earthquakes. However, the actual magnitude distribution of mining seismicity can deviate from G-R scaling and can exhibit a bimodality or other features. Therefore, we explored several other possible statistical models to fit the observed distribution of event magnitudes. Specifically, we considered the tapered Pareto distribution (8) when the magnitude m is replaced by the corresponding scalar seismic moment M = 10^1.5(m+10.73). We also considered the mixture of the tapered Pareto and Pareto distributions (10).

The application of the goodness-of-fit test indicates that the mixture of the tapered Pareto and Pareto distributions, Equation (10), can fit better the frequency-magnitude statics over a wider range of magnitudes lowering the completeness level to m_c = −1.5 (Supplementary Figures S4, S6, S7). However, this still does not fully resolve the problem that the microseismic events are fully recorded in the catalog up to this level. The histogram of the distribution of microseismic event magnitudes, given in Supplementary Figure S1, shows that there is a decrease in the number of events in each magnitude bin below the value m = −1.1. This indicates the deficiency in the numbers below this level and is typically attributed to the incompleteness of a catalog. When fitting all the three models and considering only events above magnitude m_min = −0.94, the exponential distribution (4) and the mixture of the tapered Pareto and Pareto distributions, Equation (10), fit the data best (Figure 4 and Supplementary Figure S5).

The event occurrence rate is a direct manifestation of the underlying physical processes responsible for the induction or triggering of microseismic events. The analysis of the mining seismicity shows that it is primarily controlled by the mining operations and lack any significant aftershock sequences. The cumulative number of events given in Figure 2 shows the almost steady linear increase with time for events above m ≥ −0.94 indicating a relatively constant rate of occurrence. However, the interevent time distribution between consecutive events (Figure 5) shows a deviation from the exponential distribution for short time intervals and is better described by the gamma distribution. This signifies that events that are separated by time intervals shorter than 10 min are driven by the time-dependent rate related to the mining operations. This is also observed by studying the distribution of the magnitude difference between consecutive events (Figure 6). This distribution shows a clear correlation between consecutive magnitudes and an indication of deviation from Poisson statistics on shorter time scales. However, it is also affected by the magnitude completeness of the catalog and the effect is more pronounced for lower magnitude cutoffs (Supplementary Figure S12). This effect was also observed in case of rock fracture experiments [20]. A related approach was recently suggested to use the difference between consecutive magnitudes to estimate the b-value that can handle the catalog incompleteness more efficiently [77]. However, that method relies on the assumption that the earthquake magnitudes are independent. In the case of mining seismicity, our analysis indicates a short-time correlation between magnitudes.

Clustering of seismicity is an invaluable feature that is directly related to the mechanisms of triggering and interaction between seismic events. In our work, we used the NND method (Section 3.4) to identify characteristics of the triggering process and the structure of event family trees. The GMM was utilized to investigate any possible multimodality in the distribution of the rescaled nearest-neighbor proximity between events, η (Figure 8). The GMM with two components achieved the lowest BIC = 4029.7 (Supplementary Figure S14). The corresponding intersection of the two components of the GMM was used to identify the threshold, η = 0.538, in the declustering process (Figure 8). Therefore, the microseismic events can be classified as being background events (isolated singleton events) or as clustered events. The obtained results show that the analyzed microseismic catalog has characteristics that are dominated by the occurrence of isolated singleton events.

Within the NND analysis, seismic family trees that form the clustered events can be subdivided into three main types: swarm-like, burst-type, or aftershocks [48, 49]. The swarm-like families are those chains of seismic events without branches. The burst-type families are those that have only one parent event and many children. The aftershock-type families are the combination of both mentioned types and contain first or higher-order triggered events. By applying this classification to the mining microseismicity we were able to identify 52 families (groups of microseismic events) that had more than two nodes. Using the criterion based on the value of B_I (B_I = 1 – swarm-type, B_I > 0.5 – aftershock-type, and B_I ≤ 0.5 – burst-type), from these 52 seismic family trees, 18 families had swarm-type, eight family trees had aftershock-type, and the rest had burst-type characteristics (Supplementary Figures S19–S21). This indicates that among clustered events the “burst-type” of activity dominates the clustering. It should be mentioned from these 52 family trees, 17 families had a foreshock sequence, and in the other 35 families, the root event was the largest event in the family tree. As the cluster size distribution is dominated by small clusters with a predominantly “burst-type” structure (Figure 10), the cluster leaf depth is typically short.

In this work, several statistical methods were applied to investigate the nature of microseismicity in a potash mine in Saskatchewan, Canada. Specifically, the modeling of the frequency-magnitude statistics was performed by fitting several models: the left-truncated exponential distribution (or equivalently G-R scaling), the tapered Pareto distribution, and their mixtures. By analyzing the rate of seismicity we concluded that the occurrence of events deviates from Poisson statistics on short time scales and for the lower magnitude cutoffs. The magnitude of completeness was estimated using several methods and we used m_c = −0.94 for the analysis. The results of the clustering analysis of microseismicity indicate that the majority of events can be treated as independent background events mostly driven by underground mining operations. However, there is some clustering of seismicity and the formation of limited aftershock sequences. This clustering is predominantly of “burst-type” in the terminology adopted in the NND analysis.

From the analysis performed in this work, we can draw several specific conclusions concerning the observed mining microseismicity. The frequency-magnitude distribution of seismicity exhibits a relatively high b-value (b = 2.39 ± 0.15) indicating that the seismic events are distributed in a rather narrow magnitude range above the completeness threshold m ≥ −0.94. The largest observed event during the study period had a magnitude of 0.56. The interevent triggering is also suppressed and does not show any significant cascade-like propagation of seismicity. This suggests that the probability of having even larger events is very low and does not constitute a significant hazard for the safe operation of the mine.

The data analyzed in this study is subject to the following licenses/restrictions. The earthquake catalog supporting this research is not publicly available due to confidentiality restrictions imposed by Nutrien Ltd. A research request to access the catalog data can be made to Nutrien Ltd. by following an appropriate protocol. Requests to access these datasets should be directed at: Nutrien Ltd., https://www.nutrien.com/.

MS performed research, analyzed data, wrote the computer code, wrote the manuscript, and handled the submission. MB provided the data characteristics, participated in discussions, and provided comments and contributed to the writing of the manuscript. RS formulated the problem, supervised the research, wrote the manuscript and computer code, analyzed data, and acquired funding. All authors contributed to the article and approved the submitted version.