Passive Seismic Imaging of Stress Evolution with Mining-Induced Seismicity at Hard-Rock Deep Mines

This work aims to examine the stress redistribution with evolving seismicity rates using a passive seismic tomographic tool. We compiled a total of 26,000 events from two underground mines and partitioned them into multiple clusters in a temporal sequence, each of which contains 1000 events. To image stress redistribution associated with seismicity rates, we then run the tomographic studies using each cluster to yield seismic tomograms and computed the corresponding seismicity rate. We found that high velocity anomalies grew with the increase of seismicity rates, and they switched to a shrinking tendency under low seismicity rates. Results of this study imply that seismicity rates increase with increasing stress concentration and decrease with decreasing stress concentration. This study highlights the value of utilizing passive seismic tomography for estimating stress evolution associated with the change of seismicity rates at underground mines. Our findings illuminate the applications of using mining-induced seismicity to assess stress redistribution associated with seismicity rates at hard-rock mines, providing insights into seismic hazards for deep mining.


Introduction
Mining-induced seismicity in hard-rock underground mines occurs in response to significant variation in stress such as stress buildup and release (McGarr and Green 1975;Marsan et al. 1999). Evaluating the dynamic stability of a rock mass with induced seismicity requires knowledge of how macroscopic shear cracks develop in fault zones and cause the failure of rocks. Characterizing and examining the evolution of the stress state, such as stress concentrations and relaxation, is of importance when trying to mitigate geohazards in subsurface processes. Stress fields are correlated with stress indicators, such as seismic velocity. Previous studies have shown a significant non-linear increase in seismic velocities with a small rise in the hydrostatic pressure (Mayr and Burkhardt 2006;Lozovyi and Bauer 2019). Additionally, similar studies showed that the increase in velocity is limited once the change in pressure exceeds a certain threshold (Mavko et al. 2009;Mavko et al. 2009). The assumed cause of this behavior is that all microcracks in the rock have been closed and they are no longer the dominant factor determining the behavior of seismic waves. Instead, macrostructure, which tends to remain constant with increasing pressure (or deviatoric stress change), would govern the propagation of seismic waves. According to acoustic emissions studies by 1 3 loading rock samples in lab tests, it has been found that seismic events surround the newly developed cracks, therefore, a dramatic increase in the frequency or concentration of seismic events could be used as a pre-signal reflecting the failure of rocks (Goebel et al. 2015).
By upscaling the study of seismic events from lab tests to field-scale observations, it is found that mining-induced seismicity is often a precursor to significant rock failures or rock bursts at underground mines. This characteristics of seismic events have been applied in mining to forecast rock bursts and improve mine safety (Young and Maxwell 1992;Mercer and Bawden 2005;Leake et al. 2017). Techniques applied to crustal earthquake studies, including seismic tomographic imaging, seismicity clustering characterization, and statistical modeling, have been widely used in mininginduced seismicity (Urbancic et al. 1992(Urbancic et al. , 1993Shaw 1993;Hudyma and Potvin 2010;Wesseloo 2018). Passive seismic tomography, a seismic imaging method using mininginduced seismicity to infer the velocity of structures through which waves propagate, plays an increasingly important role in stress distribution and mitigating seismic hazards as mines increase their use of seismic monitoring systems (Young and Maxwell 1992;Meglis et al. 2005;Baig et al. 2017). Studies in seismic imaging have validated that a passive seismic tomography is a useful tool for examining stress distribution and perturbation around mining (Luxbacher et al. 2008;Ma et al. 2016;Westman et al. 2017;Vatcher et al. 2018).
The last decade of research has demonstrated that a significant increase in mining-induced seismicity is usually associated with highly-stressed regions at underground mines (Morissette et al. 2017). Seismic monitoring networks assist in the timely capture of seismic sources and facilitates tracking how seismic events are triggered in response to the progression of mining and extraction (Wesseloo et al. 2014;Ma et al. 2018). In the mining field, directly measuring and mapping stress distribution is often undertaken by overcoring or hydraulic fracturing methods, both of which come at a considerable financial cost. Stress change can be monitored by uniaxial and triaxial stress cells, however the measurements can only be made at single points in the rock mass, and can, therefore, be subject to error due to local structure or changes in rock mass properties. Compared with these methods, passive seismic tomography is an efficient and economical tool for estimating stress evolution associated with seismicity at mines. The purpose of this study is to link and investigate stress redistribution with the evolution of mining-induced seismicity using seismic imaging and to determine what relationship exists between stress anomalies and the rate of mining-induced seismicity.

Data and Methods
We investigated data sets of mining-induced seismicity from two hard rock underground mines: Creighton Mine and Kidd Mine. Both mines use microseismic monitoring systems developed by ESG solutions. Previous studies summarized the geological environment and structure of the Creighton Mine (Beck and Brady 2002;Malek et al. 2008;Snelling et al. 2013). The Creighton Mine is located 20 km west of Sudbury of Ontario, Canada and the mine has been in production since 1901. The current extraction rate of the Creighton Mine is nearly 2200 tones per day. It is the deepest nickel mine in the world and the deepest operations are around a depth of 2500 m. Mining methods at the Creighton Mine have evolved over the mine life: blastholes, cut-andfill, and vertical retreat mining have been used in the upper area. The main production areas currently range in depth from 1800 to 2500 m and the current mining method is the large-diameter blasthole method. Because of the great mining depth, the Creighton Mine is seismically active. To mitigate the seismic risk, the Creighton Mine is using a pillarless slot-and-slash, a top-down center-out sequence, and dynamic ground support systems. The Kidd Mine's full name is Kidd Creek Copper and Zinc Mine and it is located 27 km north of Timmins in Ontario, Canada. The Kidd Mine is the world's deepest base-metal mine and its operations are as deep as 2927 m below the earth's surface. The Kidd Mine is operated with three shafts and the mining method used is blasthole stopping with cemented backfill. The Kidd Mine produces more than 7000 tonnes per day. Kidd Creek is based on a rich, steeply dipping volcanogenic sulphide deposit located in the Archaean Abitibi greenstone belt. The ore is hosted in felsic rocks of the Kidd Volcanic Complex and is cut by mafic sills and dykes. The orebody of Kidd Creek consists of three ore types: massive, banded, and bedded (MBB) ores; breccia ores containing fragments of the MBB ores; and stringer ores including irregular chalcopyrite stringers cutting a siliceous volcaniclastic host. These mines were equipped with similarly configured high resolution seismic monitoring networks for seismic surveys, which consisted of numerous triaxial stations and uniaxial stations. To ensure sufficient coverage of the mining region with a modest cost, a limited number of triaxial stations were employed for characterizing properties of seismic events, and many more uniaxial stations were used to record arrival times and locations of seismic events. During this investigation, ten triaxial stations and 52 uniaxial stations were configured in Creighton Mine. The local seismic network at Kidd Mine included eight triaxial stations and 23 uniaxial stations. Figure 1a-c shows the arrangement of the seismic monitoring network and the virtual distribution of raypaths at Creighton Mine. Seismic events that are distant from  station arrays and which are remote from mining are connected with a limited number of seismic stations and the data are, therefore, not well positioned to be analyzed by seismic tomography due to the lack of coverage by the seismic network. In the data processing, events located beyond the seismic network were ruled out to maintain the accuracy of the computation. A region traversed by more raypaths enables a higher accuracy in the computation. To evaluate the raypath coverage in a region, we used an index called derivative weight sum (DWS), which is the sum of total number of rays crossing a region. In addition, only seismic events being recorded by more than ten stations were considered. After filtering seismic events with this threshold, we used our own double-difference tomographic tool to implement velocity inversion by incorporating the locations of seismic events and travel time pairs from events to stations. In a similar way, the configuration of the seismic monitoring network and raypaths distribution of Kidd Mine is shown in Fig. 1d-f. We used velocity distribution of the P-wave in both a temporal and spatial fashion, displayed in seismic tomograms, to infer stress evolution and distribution in the rock masses associated with the underground mines. The speed of waves propagating in rock masses is correlated with the effective stresses (Holcomb 1981;Young and Maxwell 1992;Mavko et al. 2009). Elevated stress is an important factor in raising bulk modulus, which determines the P-wave velocity travelling through the solids. Previous studies validated the quantitative relationship between these elements as where V P is the P-wave velocity, ρ is the density of rock, K is the bulk modulus, and μ is the shear modulus (Mavko et al. 2009).
We formed a double-difference algorithm, built on a series of double-difference algorithms including hypoDD (Waldhauser and Ellsworth 2000) and tomoDD (Zhang and Thurber 2003) and incorporated it into a Matlab library with an implementation of travel time and event location inversions to obtain the optimal velocity structures. The double-difference algorithm generated the double-difference residual by extracting the catalog of arrival time, which is expressed as, is the difference in the observed arrival times between events i and j, and (t i is the difference in the computed arrival times between events i and j. It is noted that this double-difference tomographic tool is capable of performing the velocity inversion along with cal relocating seismic events (Zhang et al. 2019). The locations of seismic events are improved through this iterative process and every iteration of velocity structure adopts these new locations after the relocation iteration. The iterative results of relocation are not shown in this study due to the emphasis of this study on the velocity change in tomograms, but the relocating function contributed to the improved the velocity structure for seismic tomograms. The process of relocating seismic events in Creighton Mine was introduced in a previous study (Ma et al. 2016). Other studies applied collaborative localization methods using analytical and iterative solutions for microseismic sources in underground mining. These studies filtered the abnormal arrivals using analytical solutions combined with the iterative method without premeasured average velocity to improve the location accuracy (Wang et al. 2018;Dong et al. 2019). These studies emphasized constructing the velocity discontinuities, including Conrad and Moho discontinuities. It was reported that the accuracy of the tomographic study on hole-contained structures can be improved by a velocity-free source location method (Dong et al. 2020). The likely influences from the goafs were not included in this study due to a certain limit on the daily production data. By contrast, the fundamental theory of the double-difference tomographic algorithm is to locate two close events and subtract their similar raypaths for eliminating errors along the raypaths. The raypath coverage of a region is evaluated by the DWS (Derivative Weight Sum) value, which is an index of the sum of rays crossing a region. In addition, we set a threshold of the minimum number of stations recording a seismic event at ten, so that only seismic events being recorded by more than ten stations were considered. All other events that failed to meet this threshold were ruled out and this threshold was applied to both Creighton Mine and Kidd Mine. We performed velocity inversion with an estimated initial velocity model that assigned a uniform velocity value to all nodes of a 3D mesh network for Creighton Mine (Fig. 2) and Kidd Mine (Fig. 3), respectively. More details, including a trade-off analysis for finding the optimal damping of velocity inversion and robustness and the method used for applying passive seismic tomography in underground mining was described in a recent study (Ma et al. 2019). In our computation for Creighton Mine and Kidd Mine, an initial velocity model was set up as a 3D grid of points with a uniform P-wave velocity. We delineated 40 layers in-depth and each layer was subdivided into 40 × 40 grids for building the tomographic model (   process for obtaining an initial velocity value for Creighton Mine and Kidd Mine, respectively. Velocity inversion and event relocation were iteratively processed until the optimal solution with a reliable accuracy was reached. A series of damping values was tested to select the optimal one, which was determined and applied in the iteration computation to avoid overfitting. Velocity structures were displayed after interpolating a 3D grid of points on results of P-wave velocities extracted from the last iteration. The process of this tomographic study in mining-induced seismicity is summarized and is shown as a flow chart in Fig. 5. The catalog of mining-induced seismicity needs to be partitioned in time and/or space for tomographic studies. Two approaches are suitable for our temporal tomographic study with quantifying seismicity rates. One is to partition the total data set of seismic events in the temporal sequence into multiple clusters, each of which has the same number of events. The other is to divide the full time span of the catalog containing all seismic events into multiple consecutive time bins of equal lengths of time. Events that fall into the same time bin are grouped as a single cluster. In this study, we used the former method to partition seismic events into clusters with floating time lengths to ensure similar resolutions of tomograms and calculate the rigorous seismicity rates within time periods of different lengths. Using an identical number of events ensure that each seismic tomographic image has a similar resolution for further comparison across them. Considering that seismic rates of both Creighton and Kidd Mines fluctuated dramatically over a set time interval, we allocated the same number (1000) of temporally consecutive events into a single group to ensure that all tomograms could be generated with a similar resolution. Although the number of waves propagating through a rock mass and the accuracy of computation is the main factors that govern the resolution, the identical number of events in our setup would produce a very similar number of raypaths, guaranteeing the similar resolution of tomograms for further comparison across different periods.
The data set of seismicity cataloged by Creighton Mine contained seismic events that were triggered between June and September 2011. It consisted of 251,436 P-arrival times from 14,000 microseismic events. In addition to using a threshold of raypaths, events studied in double-difference tomographic inversion were selected meeting the criteria that the number of stations picked were more than eight. That said, only events being captured by more than eight stations were eligible for the tomographic computation. Similarly, seismic events recorded at Kidd Mine occurred between September 2008 and May 2011. Seismic events filtered using raypaths and station number included 180,568 P-arrival times from 12,000 microseismic events.
Affected by stress perturbations caused by excavations and production blasts, the seismicity rates of Creighton Mine fluctuated significantly. The cumulative number of microseismic events over time in Creighton Mine is shown in Fig. 6. Creighton Mine experienced about 200 events per day on average over the long-term. The days with high seismicity rates, such as July 6th, July 10th, and August 4th, can be identified by relative slopes of segments on the curve. In tomographic studies on Creighton Mine, seismicity data was divided chronologically into several groups and each group included 1000 microseismic events, thereby providing adequate raypath coverage for generating the tomograms with a reliable resolution.
At Kidd Mine, several large blocks, defined by major faults, significantly displaced by blasts. Stress conditions on faults are influenced by deformation of nearby rock mass and remote rock mass through a gradual 'domino' effect. Historic excavations that change confinement and contribute Mine. The slope of the red line, which were fitted using robust linear regressions, reflects the average velocity of these raypaths. The slope of linear regression lines is used as the initial P-wave velocity that is applied to grid nodes of velocity models for the inversion to induce deformation on structures, leading major faults to slip incrementally and other critical structure rupture. Abrupt structural failures are usually associated with high seismic release within stopes. Affected by deformation of rock mass and fault slip, regional stress fields tend to evolve with shear-induced seismic energy release. Imbalanced deformation of rock mass exacerbates stress concentrations in regional stress fields, resulting in the surge of seismicity. However, the energy released with seismic activity only took up a small fraction of the total potential energy of the whole system at Kidd Mine, possibly contributing to the difference in seismicity frequency between Kidd Mine and Creighton Mine. The long-term seismicity rate at the Kidd Mine (Fig. 7) was only 20 events/day in the studied period from September 2008 to May 2011, which was roughly nine times lower than Creighton Mine (Fig. 6). In a similar format of processing data applied to Creighton Mine, we split the data set of seismicity of Kidd Mine into several groups within a temporal sequence, and each group included 1000 microseismic events. Each sub data set was run independently to implement velocity structure inversion and seismic events relocation.

Results
We quantified and examined temporal seismic rates of Creighton Mine and Kidd Mine using historical seismicity data. The seismic rate at Creighton Mine changed significantly, ranging from 5 to 291 events/hour (Table 1). Comparing with Creighton Mine, Kidd Mine experienced a relatively lower seismicity rate: the lowest seismic rate was 6.93 events/day and the highest one was 33.85 events/day at Kidd Mine (Table 2). Although seismicity rates per unit of time varied significantly between Creighton Mine and Kidd Mine, both mines followed the relationship that high stress precedes the occurrence of increased seismicity.
To identify the correlation between seismic rates and stress distribution, we generated seismic images using the doubledifference tomographic tool we developed. Temporal seismic

Creighton Mine
The seismic rate fluctuated steeply over the period from June to August 2011. As shown in Fig. 8i, we examined the crosssection at depth of 2340 m (7700 ft) which included the main production level in 2011. Microseismic events were consecutively partitioned where each cluster contained 1000 events and each cluster, therefore, was formed over different lengths of time. Each tomogram, being plotted from results of velocity inversion and structure reconstruction using the source of 1000 events, visualizes velocity distributions in a series of consecutive periods to illustrate stress evolution. Microseismic events are mainly localized in the area of Easting 1300 m and Northing 1850 m from (Figs. 8,9). In addition to the main swarm of seismicity, a few microseismic events scattered outside the working drifts. Across all tomograms, Fig. 6 The cumulative number of seismicity in long-term (a) and short-term target periods (b) at Creighton Mine. The magnitude distribution (c) shows that the seismicity count decreases with increasing moment magnitude we recognized that microseismic events were more spatially concentrated as clusters rather than evenly distributed. Further, it was found that high stress precedes the occurrence of increased seismicity by imaging stress redistributions with a spatial distribution of seismicity. The main region with high levels of velocity is contoured around mining production drifts and formed a 'Y' shape in the first period (Fig. 8a). From the first period (Fig. 8a) to the second period (Fig. 8b), the rate of seismicity declined by 27% and the previous high velocity in 'Y' shape returned to the background velocity, which was the initial velocity value assigned to the velocity model. Seismic rates reached their peak at 193 events/ hour between 4:43 AM and 9:54 AM on July 6th (Fig. 8c). Along with the boost of seismicity, a highly-stressed region appeared and concentrated on the center of the drifts. The high seismicity rate dropped by about 100 events/hour during the next period from 09:54 to 19:48, July 6th (Fig. 8d) reducing by about 45% from the peak rate at 192 events/ hour (Fig. 8c). Tomograms indicate that the previous highly-stressed area of last period expanded, leading to the Fig. 7 The cumulative number of seismic events in long-term (a) and short-term target periods (b) at Kidd Mine. The magnitude distribution (c) shows that the seismicity count decreases with increasing moment magnitude domination of high stress in most regions (Fig. 8d). Stress evolution on these tomograms implies that the most pronounced stress variations appeared along the direction from northwest to southeast. Interestingly, the location of event groups failed to migrate immediately even though the high stress region shifted, suggesting a stagnation effect of seismicity with stress shift. In contrast, the induced seismicity swarm was located beyond the high stress region, and was more likely to concentrate in the area between the high stress and low stress zones. With the small decline followed by an increase of seismicity rate after 19:49, July 6th (Fig. 8d), the existing high stress slightly decreased and returned to a fairly strong level until the seismicity rate recovered on July, 10th. Several regions showed low velocity plumes, suggesting stress relaxation in these areas (Fig. 8e). Seismicity rate stabilized and was still higher than 30 events/hour during Period C6 (Fig. 8f) and C7 (Fig. 8g). Accordingly, the velocity distribution of Period C6 was nearly identical with Period C7. The seismicity rate continued to decrease to as low as 14 events/hour in Period C8 (Fig. 8h). Associated with the low seismicity rate in Period C8, the region in the center of mining drifts, which were previously dominated by the high velocity volume, decayed to the background velocity. High stress and the occurrence of increased seismicity were well correlated and seismicity distribution shifted following the newly formed high stress.
To evaluate whether the stress concentration would be associated with higher seismicity rates as well as stress reductions with lower seismicity rates, we examined the other main period from July 29th to August 28th, which contained sub-periods with more variable seismic rates. In these periods, the active volumes with high seismicity rate first appeared on August 4th, 2011 in Creighton Mine, and the tomograms coupled with microseismicity distribution around the seismically active periods were analyzed. The seismicity rate increased significantly and reached 291 events/hour on August 4th, 2011. The tomograms indicate that two roughly symmetrical regions were associated with the high velocity, which also enclosed a small area with low velocity, implying that a stress concentration dominated the region when the seismicity rate was at its maximum ( Fig. 9b). A single swarm of seismic events in this period was located on the intermediate velocity area between the two high velocity areas. We found three consecutive periods after August 5th where the seismicity rate was as low as about five events/hour (Fig. 9d-f). In these periods, a few low-stressed zones were identified and microseimic events were dispersed throughout them. In the first of these three seismically stagnant periods, a strong low velocity region emerged to the west of mining opposite a substantial high velocity region to the east of mining (Fig. 9d). From the first period (Fig. 9d) to the last period (Fig. 9f) during these three consecutive periods of low seismicity, low velocity regions gradually increased and high velocity areas correspondingly decayed, leading to the dominance of low velocity in the last period (Fig. 9f).

Kidd Mine
The overall seismicity rate of Kidd Mine was significantly lower than Creighton Mine. A sum of 12,000 seismic events, occurring from September 2008 to May 2011, were analyzed. We used the same approach as used for Creighton Mine and split seismicity into multiple event groups, each of which consisted of 1000 seismic events. Seismicity at Kidd Mine showed signs of high shear displacement and increased deviatoric deformation along with existing structures with a large release of seismic energy. The seismicity rate ranged roughly from 7 to 30 events each day and details of seismicity rates are listed in Table 2. Velocity structure of each group was generated by performing double-difference  Figure 10 shows the subsurface structure and locations of tomograms cut along a vertical plane visualized for Kidd Mine. Stress distribution was indicated by the velocity distribution of these tomograms, which were superimposed with microseismicity to indicate possible correlations between velocity anomalies and seismicity locations. Tomograms calculated from events groups in a temporal sequence were compared to estimate the stress evolution as periods proceeded. Tomograms of Creighton Mine are arranged in a horizontal cross section. In contrast, tomograms of Kidd Mine we studied are laid out in a vertical cross section due to the extensive vertical distribution of seismicity. Tomograms calculated from Events Group K1 to Events Group K12 are shown in Fig. 11. We found that considerable microseismic events were located in the high velocity regions of tomograms, implying that the occurrence of seismicity was potentially capable of indicating stress buildup. From Events Group K1 (Fig. 11a) to Events Group K4 (Fig. 11d), it is evident that high velocity regions followed the migration of the seismicity. Specifically, regions which displayed high stress were followed by the occurrence of increased seismicity in the same areas. Some swarms of seismicity appeared to be in the boundary of high velocity regions rather than in their centers. For example, the swarm of Event Group K4 was located in the zone between two major high velocity areas (Fig. 11d). Swarms of Event Group K5 (Fig. 11e) and K6 (Fig. 11f) emerged in the zones between high velocity regions and low velocity regions. Within the low velocity areas, these two groups contained few events, whereas spatial centers of seismicity swarms appeared near the border of high velocity regions (possibly identifying regions of high deviatoric stress). It should also be noted that the time period in question was abnormal, in that a pair of large damaging rockbursts occurred in January and June 2009, resulting in an extended period of rehabilitation, which affected the mine production sequence for a period of almost 18 months, ending in approximately middle 2010.
According to seismicity rates (Table 2), Event Group K4 (Fig. 11d) and Event Group K8 (Fig. 11h) experienced the two highest seismicity rates. Group K4 immediately preceded the Mn 3.2 major seismic event of June 9, 2009, and Group K8 was associated with the resumption of mining in the satellite Greywacke lens, which was ultimately determined to be a contributing factor in the development of the major seismic events of 2009, and was associated with the subsequent large magnitude events in mid-2011. It is also interesting to note that Group K8 is coincident with mining in the south abutment which is believed to have failed a regional pillar near the 6400 Level. On the other hand, Event Group K10 (Fig. 11J) and K11 (Fig. 11k) correlated with the two lowest seismicity rates. The remaining periods were recognized as having moderate seismicity rates. High velocity in these two groups with the highest seismicity rates displayed a larger amplitude and more concentrated manner, indicating that higher seismicity rates were associated with elevated stress concentrations. In addition, high velocity anomalies were still found in periods with low seismicity rates (Fig. 11j, k). The notable characteristic of these two groups was that low velocity anomalies were identified and these low velocity regions were free of seismicity. An explanation for this phenomenon is that stress drop or relaxation tends to decrease the occurrence of seismicity in the same area, leading to an average low-velocity rate within large spatial ranges. With the reduction of the seismicity rate from K8 (Fig. 11h) to K13 (Fig. 11i), the first major low velocity zone initiated at Easting 65,840 m, Depth 1300 m (Fig. 11i). This low velocity zone continued to develop and expand with time moving forward and stabilized in the next two periods K10 (Fig. 11j) and K11 (Fig. 11k). The second low velocity zone formed in K11 and had a quite different growing pattern compared with the first one. Its previous location is shown in Fig. 11j was governed by a high velocity anomaly and abruptly converted to the second low velocity anomaly (due to the mining of a stope causing relaxation), which initially appeared in K11 (Fig. 11k). The region adjacent to the east of the second low velocity anomaly was occupied by a high velocity zone which expanded in the depth direction, corresponding to the core of a regional pillar between the main ore zone and the satellite Greywacke lens further in the southwestern abutment. A few widespread seismic events were found in the expanded area associated with the growth of this high velocity region (Fig. 11k). In the next period (Fig. 11i), the seismicity rate returned to a moderate level and both low velocity anomalies disappeared, proving that low-stress zones were more likely to exhibit low seismicity rates. The low velocity zones likely disappeared as the result of filling mining voids (stopes) within the south abutment, which provided confinement to the pillar, and allowed the transfer of stress from the rock mass, into the fill.

Discussion and Conclusions
This study investigated the evolution of stress with fluctuation of seismicity rates and showed that stress concentrations increase seismicity and stress relaxation decreases seismicity. The most prominent stress concentration appeared in the period right after the period with the highest seismicity rate, suggesting that seismicity is triggered by preceding stress concentrations. Induced seismicity accordingly can be used to predict stress concentrations temporally. However, regions with stress concentrations are not always the exact locations where the preceding seismicity occurred. As a result, it is not reliable to forecast regions with stress concentrations merely based on spatial information of seismicity. Geological structures including microscopic fractures, weakness planes, and joints inevitably exist and influence the rock mass behaviors with complex structural interactions within the orebody. Small changes or disturbances in confinement and clamping forces on structures are likely to significantly affect mines and lead to seismicity. Passive seismic tomography assists in determining whether stresses return to equilibrium or form new areas of stress concentrations by imaging stress redistribution at a high resolution. Thus, passive seismic tomography modeled with a microseismic system data can facilitate refinement and calibration for other numerical stress or deformation models. According to the comparison across multiple periods, we found that stress distributions keep varying with the evolution of the seismicity rate. Areas of stress concentrations tend to aggregate during periods with high seismicity rates (possibly indicating the coalescence of joints, and the development of regional failure zones within the rock mass). Expanding regions of stress relaxation form in the periods with low seismicity rates. It is inferred that the increase of seismicity rate was associated with elevated stress, implying that high stress promoted the occurrence of these seismic events. In studies of natural earthquakes, it is recognized that small events enhance stress at the epicenters of them (King et al. 1994). Stress changes triggering seismicity are very small: stress change ranging from 0.1 to 0.3 MPa is sufficient to lead to seismicity. On the other hand, stress drop in the same amount could significantly decrease seismicity. Results from this study agree with the studies related to natural earthquakes that small events can increase stress. In different periods with similar seismicity rates, several regions appeared to display very similar stress distributions. For Fig. 11 Tomograms of Kidd Mine cut along a vertical plane, facing north, with superimposed seismicity (purple circles) example, the similarity between K1 and K12 (Fig. 11a, l) at the Kidd Mine implies that a memory effect probably exists in the regional stress. After the release of accumulated strain energy in the form of induced seismicity, the regional stress would recover to the previous stress distribution, unless the mining which occurred over the intervening time period was sufficient to alter the mine geometry and permanently change the regional distribution of stress.
In crustal earthquake studies, it has been proposed that elevated seismicity rates are correlated with regions of peak transient stresses (Svetlizky et al. 2016). In this study, the evolution of stress distribution mainly reflects static stress change, which dominates the triggering of induced seismicity at underground mines and has a larger amplitude than dynamic stress change. It is well known that static stress change governs the spatial distribution of seismicity (King et al. 1994). Thus, analyzing the spatial distribution of seismicity with tomograms provides a viable way to identify and validate correlations between stress distributions and spatial distribution of seismicity, as well as seismicity rates in multiple periods.
We examined whether high stress grew with the increase of seismicity rates, and the correlation between the stress evolution and migration of seismicity. Tomograms from Creighton Mine and Kidd Mine indicate that most stress distribution changes occur prior to the migration of seismic events. For example, by comparing stress evolution and redistribution of seismicity from Fig. 8c to Fig. 8d, we found that in Fig. 8c few seismic events were located in the major highly-stressed region in the center of production drifts, whereas in the next period more seismicity appeared in the center of production drifts where the major high-stress anomaly existed (Fig. 8d). Another example in Kidd Mine agreed with this pattern as well. From K7 (Fig. 11g) to K8 (Fig. 11h), although a highly-stressed region in Northing 1200 m had diminished to the background velocity, several seismic events were still found in the area that was previously occupied by the high velocity anomaly. This evidence suggests that stress redistribution usually precedes the migration of seismicity. That said, mining-induced seismicity is more likely to appear in, or adjacent to, highly-stressed areas and it is triggered after the formation of highly-stressed regions, usually by the removal of significant volumes of rock, which locally changes the deviatoric stress regime and triggers rock mass failure.
It is important to know that this study has limitations. To calculate the seismicity rate and run tomographic studies, the seismicity was arbitrarily divided into multiple volumes and each of them contained the same number of events. In engineering practice, the seismicity rate is affected by the progress of extraction and production blasts in mining. It is important to partition the seismicity to ensure that the temporal bins of seismicity comply with the production schedule of mining. The tomographic study could be more meaningful when stress evolution is linked with the production cycle of the mine, providing a geomechanical perspective to the influence of the extraction of ore reserves.
Future research surrounding passive seismic tomographic study in mining is intended to investigate more factors affecting seismicity rate, in particular: (1) considering energy release in addition to seismicity event rate, and (2) evaluating the effect of timing and size of production blasts.