1 Introduction

Stress perturbation due to the mining process can potentially cause rock failure in the form of new fractures that disrupt the stability of the geological structure (Dong et al. 2016, 2019; Lu et al. 2016; Zhao et al. 2017). The fracture propagation process at cave areas in underground mines using the block caving method caused by the induction of mining activities needs to be monitored so as not to cause rock damage. As a risk to bear, uncontrolled fractures in the rock will cause tunnel collapse, rockburst, slope instability in open pit mines, and seismic hazards (Bommer et al. 2006; Cheng et al. 2009; Jiang et al. 2016; Kinscher et al. 2017; Melati et al. 2022; Sausse et al. 2010; Wang et al. 2018). The potential hazards in mining above must be monitored to ensure the continuity of mining production, worker and equipment safety, and infrastructure security (Kinscher et al. 2017; Mendecki 1997; Mercier et al. 2015; Pierce et al. 2020; Sainsbury 2010).

With the highly dynamic nature of rocks during the mining production process, microseismic monitoring is required, which can provide input for making decisions about an effective and safe mining process (Ghosh and Sivakumar 2018; Mendecki 1997). Previous studies have shown the application of microseismic for monitoring ground failure, longwall stability, roof fracture process, hazard assessment, rock mass stability analysis, rockburst prevention etc. (Chambers et al. 2015; Cheng et al. 2009; Ghosh and Sivakumar 2018; Gochioco et al. 2012; Jiang et al. 2016; Qi et al. 2015; Shen et al. 2013). Monitoring fluid injection and fracture zones is research that is often studied in microseismic monitoring in the geothermal field (Bommer et al. 2006; Cuenot et al. 2008; Sausse et al. 2010). Meanwhile, in the oil and gas industry, microseismic is used in various ways: exploration activities, hydrofracturing monitoring, reservoir characterization, Carbon Capture Storage, etc. (Chen and Huang 2020; Maxwell et al. 2010; Maxwell and Urbancic 2001; Rafiq et al. 2016; Wu et al. 2017).

Cave propagation is the process of developing fractures that are formed due to the influence of gravity caving and stress caving so that the rock mass shifts due to the loss of some or all of its cohesion value (Brown 2002; Pierce et al. 2020; Sainsbury 2010). Cave propagation is a complex process related to the heterogeneity of rock masses, geological structures, and insitu stress (Campbell et al. 2016.). The caving process can trigger seismic induction in the mining area, resulting in rock stress changes. Stress accumulation in underground mines generally occurs in the abutment and seismogenic zones (Duplancic and Brady 1999.). The high-stress concentrations in the abutment and seismogenic zones result from the uptake of the rock mass. Ideal conditions are expected according to the mine design; cave propagation will propagate in the direction of updip intrusion. If the induced stress is very high, the rock mass can experience collapse, which can affect the direction of cave propagation. Uncontrolled cave propagation can result in an over break of the lithology boundary. This is not expected because it will cause the high-grade ore to mix with the low grade ore (dilution). To find out this possibility, it is essential to monitor cave propagation.

This research was conducted at PT Freeport’s Deep Mill Level Zone (DMLZ) underground mine, the deepest underground mine in Indonesia (Melati et al. 2022; Rumbiak et al. 2022). DMLZ is one of the deepest underground mines in Indonesia located in Papua Province with a depth of ≥ 1780 m from the surface (Fig. 1). The DMLZ is located at an elevation of 2590 masl (extraction level)−2600 masl (undercut level) with an area of 630,000 m2. The DMLZ mine is located 410 m below the Deep Ore Zone (DOZ) mine. The DMLZ is located 3 km to the southeast of the Grasberg open pit. The DMLZ is part of the East Erstberg Skarn System (EESS) deposit. The EESS intrusion is a single copper/gold Diorite intrusion and skarn body in sedimentary rocks. A strong rock mass with a UCS value of up to 156.5 MPa, complex geological structure, deep mining category (> 1500 m), incline intrusion type, and a footprint width of ~ 450 m result in high-stress perturbation making DMLZ need to be very well monitored (Widodo et al. 2019).

Fig. 1
figure 1

Tectonic map of the Papua area, Indonesia. Major inland fault and subduction zone are plotted. The Insert map shows the location of the Deep Mill Level Zone (DMLZ) and its surrounding blocks operated by PT Freeport Indonesia

Full size image

In underground mining operations, cave monitoring is crucial for ensuring the safety and efficiency of the mining process. One promising technology showing success in cave monitoring at underground mining is Time Domain Reflectometry (TDR). Time Domain Reflectometry is a geophysical technique that measures electrical signals reflected from underground structures (Lanciano and Salvini 2020; Shentu et al. 2022). By utilizing Time Domain Reflectometry, mining operators can track changes in the cave propagation and identify potential risks. TDR is an electrical pulse testing method initially created to pinpoint fractures in power transmission cables. In recent years, this approach has been modified for the purpose of monitoring rock collapse resulting from underground mining activities (Liu et al. 2018; O’Connor et al. 2004). Coaxial cables are embedded into a rock mass anticipated to cave, and as the rock shifts over time, it causes deformation in the cable. Subsequently, alterations in cable deformation lead to variations in TDR pulse reflection patterns which can be analytically linked to the extent and position of rock movement (Benson and Bosscher 1999). To monitor cave propagation in the DMLZ area, microseismic and TDR monitoring is conducted.

In this research we use microseismic events occured in the DMLZ area are monitored by 84 stations scattered around the mining facility area within 57 days between August 29 to October 24, 2020. In that period, TDR observed a significant cave propagation especially, toward NW direction (Fig. 2). According to TDR data taken on September 13, 2020, the cave’s outer limit is indicated by the brown line, the blue and yellow color is the outer limit of the cave on 27 September and 11 October 2020, and the green line is the outer limit of the cave on 24 October 2020 (Fig. 2). The travel time from the microseismic source to the stations were observed from time to time. The change in travel time is the used to see changes in rock mass caused by the caving process. Changes in the shape of the cave will have implications for the direction of the ray path, while changes in velocity in the underground mine area reflect changes in the rock mass and stress distribution.

Fig. 2
figure 2

Visualization of the DMLZ cave propagation section based on 57 days of Time Domain Reflectometry (TDR) data. a Vertical cross-section looking to NE at elevation 2620 to 2840 masl, b horizontal cross-section at elevation 2700 masl. The color of the line shows the progress of the cave every two weeks

Full size image

2 Geology

The world’s largest Cu-Au skarn system is located in the Ertsberg–Grasberg region of Papua Province, Indonesia. The DMLZ underground mine is included in the East Erstberg Skarn System (EESS) deposit and controlled by two main fault EES fault and Yellow Valey syncline fault (Fig. 3). The EESS intrusion is a single intrusion of diorite with copper-gold and skarn deposited bodies in sedimentary rock. EESS is a vertically continuous porphyry and skarn-type mineralized zone that develops in the contact zone between Erstberg diorite (Te) rocks and carbonate and siliciclastic wall rocks formed along the northern contact of the Erstberg intrusion complex (Cloos et al. 2005). The Yellow Valley Syncline’s southern part is home to the Ertsberg Skarn System (Fig. 3). Gunung Bijih Timur, Intermediate Ore Zone, Deep Ore Zone, Mill Level Zone, and Deep Mill Level Zone (Mertig et al. 1994; Rumbiak et al. 2022). The layers are oriented toward NW - SE which is influenced by the dominant structure in the form of folds and faults. The syncline structure in the research area has an NW–SE direction and an up-fault structure in the same direction. There is also a paired shear fault in the NW–SE and NE–SW directions. The position of the rocks varies with a slope angle of 40o–75o with the category of moderate to steep layer slope (Fig. 3).

Fig. 3
figure 3

(Modified from Fatimah et al. 2017; Rumbiak et al. 2022)

Geologic map of the Deep Mill Level Zone (DMLZ) at 2600 masl.

Full size image

As a result of the interaction between igneous rocks and sediments, the alteration rocks and rock units that make up the DMLZ underground mine have different mineral compositions. Skarn, hornfels, and low grade ore with significant changes are a few of the modifications (Fatimah et al. 2017). Dolomitic limestones from the Tertiary Waripi Formation are altered, and copper-gold skarn mineralization is found there. Although in some instances substantial endoskarn alteration forms that are difficult to differentiate from nearby sedimentary rocks, clear and easily identifiable rock boundaries are frequently encountered between diorite and skarn in the footwall. In a hanging wall, the transition between stone and marble is typically visible. The condition, original nature, and location of extant main faults are frequently covered by magma intrusion, hydrothermal breccia, and alteration. The brittle-fracture behavior of strong rock mass, combined with complex structures in contact and altered zones, triggered unstability problems. Spalling and rockburst are the common risks related to rock mechanics challenges in this geological setting.

3 Data and methods

The data is the DMLZ underground mine microseismic catalogue for 57 days (from August 29 to October 24, 2020) consist of 14,821 events. The data were obtained from 84 seismometer stations scattered in the main underground mine facility area (Fig. 4). Eighty-two (82) seismometers concentrated at 2600 to 3100 masl, and two (2) stations at an elevation of 3600 masl. Figure 4 shows the position of the seismometer, microseismic event, cave, and haulage shaft in the DMLZ mine facility area. The purpose is to analize the significant changes in cave propagation which was observed by TDR data.

To see the evolution of velocity and stress distribution, we divide the data into 4 (four) data subsets. Subset 1 data was from August 29 to September 12 totalling 4406 events that consist of 38,812 P phases and 23,179 S phases. Subset 2 data was for September 13–26 totalling 2887 events that consist of 28,806 P phases and 18,943 S phases. Subset 3 data was from September 27 to October 10, 2020 totalling 3710 events that consist of 35,230 P phases and 26,672 S phases. Subset 4 data was from October 11 to October 24, 2020 totalling 3818 events that consist of 73,417 P phases and 64,948. The moment magnitude ranged from − 1.6 to 2.0 Mw.

Fig. 4
figure 4

Seismometer distribution map in the DMLZ underground mine area, the black line represents the existing haulage shaft, the brown pyramid represents the cave, the seismometer station distribution is indicated by the yellow box, and the distribution of microseismic events from 29 August to 13 September 2022 of 4406 is shown with white to black color dots the bigger the dot shows the moment magnitude

Full size image

Travel time Double—Difference tomography method is used to see velocity evolution in the vicinity of the cave. Stress distribution is estimated from the distribution of P and S wave velocity. In addition to currently interpretation of changes in P waves for underground monitoring (Ma et al. 2018, 2019; Wang et al. 2018). In this research, we analyzed S-wave changes to quantify whether cave propagation had an impact on changes in Vs and Vp/Vs ratios. We used the TomoDD program to invert the 3-D velocity structure and relocate the hypocenter.

As research area, consists of 5 underground mines, the wave travel time data is selected to minimize microseismic events originating outside the DMLZ. Data selection aims to sort outlier events originating from regional activities or other mines around the DMLZ mine. As well as to data selection is made by making a plot between travel time and distance travelled and based on the velocity of rock information from the results of rock tests in the laboratory (initial model velocity). Plotting the data between travel time and distance travelled obtains a linear gradient. The selection of travel time data used in this research is data that has a gradient velocity of the P wave between 4.8 and 6.3 km/s and for the S wave of 2.7−3.5 km/s. The results of data selection are shown in Figs. 5 and 6.

Fig. 5
figure 5

P wave velocity gradient (4.8−6.3 km/s), The range of value intervals between the green and blue lines is taken, while the values ​​outside the blue and green lines are ignored

Full size image
Fig. 6
figure 6

S-wave velocity gradient (2.7−3.5 km/s), The range of value intervals between the green and blue lines is taken, while the values ​​outside the blue and green lines are ignored

Full size image

3.1 Tomographic inversion

Synthetic travel times for P and S waves were resolved using a pseudobending technique via an updated 3D velocity model (Zhang and Thurber 2006). Before carrying out the tomographic inversion process, we parameterized the model. Model parameterization includes grid node size, damping, smoothing, threshold, and initial model velocity. The grid nodes used in the TomoDD program are divided into two parts, namely inner and outer nodes. Inner nodes are nodes on the inside that are denser in size, while outer nodes are on the outside and function as coverage with a larger grid size. We vary the size of the grid nodes (inner node and outer node) to see the maximum resolution that can be used. This is done because the underground mining area has a small space (on the order of meters) so a node size that can represent the research area is needed. Based on the results of the grid node analysis that we have done for all data subsets using a grid size of 40 m × 40 m × 40 m.

The number of grid nodes for the directions nx, ny, and nz are 25, 26, and 18 respectively. The size of the inner node is 40 × 40 × 40 m around the caving area, then widens to 100 m and the outermost point is 5 km in the x direction, 3 km in the y direction, and 1.5 km in the z-direction (Fig. 7). TomoDD uses the Conditional Number Damping (CND) value to show the stability of the inversion. Stable inversions are in the range of CND values ​​between 40 and 80 (Waldhauser 2000). CND value analysis was carried out by making an L curve (trade-off) between the variance data and the variant model with input damping that varied from 50 to 500 (Fig. 8). Figure 7 shows the L curve (trade-off) for choosing the optimal damping value for the case of linear interpolation and partial derivatives of travel time with model parameters calculated directly. The optimal damping value based on CND is in the range of 200–300. The optimal damping value that is applied to all subsets of the data that we invert is 300.

Fig. 7
figure 7

Configuration of 3D grid nodes model for tomographic inversion consists of inner nodes and outer nodes. The brown pyramid is the caving area, and the red lines are the haulage shafts of the DMLZ.

Full size image
Fig. 8
figure 8

L curve (trade-off) between variant data vs. variant model to select the optimal damping value, the damping value includes the Condition Number Damping (CND).

Full size image

We use a 3D initial velocity model obtained from geological information based on observations, mapping, and measurements of intact rock coring results (Prahastudhi et al. 2018; Rumbiak et al. 2022). Intact rock samples from each lithology were measured for their velocity values using an ultrasonic test. The lithology in the research area consists of skarn, diorite, and limestone. The velocity of P waves in skarn, diorite, and limestone measured using ultrasonic velocity is 6.23 km/s, 5.59 km/s, and 4.91 km/s. For S-wave velocity, each of Skarn rock, diorite, and limestone is 3.50 km/s, 3.175 km/s, and 2.84 km/s (Fig. 9).

Fig. 9
figure 9

Initial model of P and S wave velocity, the red color represents the velocity value of diorite, the purple color represents skarn, and the blue color represents the velocity of limestone lithology. The crown color in the middle of the image is a cave

Full size image

3.2 Resolution test

Resolution tests was used to estimate the resolution of tomographic inversion results. We conducted a resolution test using the Checkerboard Resolution Test (CRT) and the Derivative Weight Sum (DWS). CRT is a synthetic resolution test performed to determine the accuracy of the parameters used in tomographic inversion. The synthetic wave travel time data is generated from a velocity model like a chessboard. Resolution test was conducted to determine spatially which area can be interpreted. The area with a similar pattern between the input model and the inversion result show that the area can be interpreted properly.

The perturbation velocity model used for the CRT test resolution is ± 10% of the initial velocity model (Fig. 10). The initial CRT velocity model is made by providing alternate low and high velocity anomalies relative to a certain initial seismic velocity model. The initial velocity model used is fixed for a P wave of 5.59 km/s and a Vp/Vs ratio of 1.76 km/s. The travel time obtained from the forward modeling process obtained from the CRT velocity model is calculated using the location of the hypocenter and station from the earthquake catalogue for each data subset.

Synthetic data of travel time is considered as observational data which is then inverted using the initial model of seismic velocity and hypocenter and station location data. The synthetic data obtained is then searched for the event pair. The resulting forward travel time is then inverted again to get the value of Vp, Vs and the ratio of Vp/Vs. We use a pseudobending ray tracing algorithm to trace the ray (ray trace) between the event and the station (Fig. 11). In the CRT and DWS tests, station inputs and subset 2 microseismic events were used. The CRT and DWS results are shown in Fig. 12.

Fig. 10
figure 10

Perturbation ± 10% of the initial 3D seismic velocity model, the blue color shows the value + 10% of the initial velocity model value, while the red color indicates—10% value of the initial velocity model. The brown area represents cave

Full size image
Fig. 11
figure 11

Raytracing from events to subset 2 stations for the CRT and DWS tests. The yellow squares represent stations, the gray lines represent ray tracing from events to stations, and the brown area represents cave

Full size image

The CRT and DWS images were taken on horizontal section at an elevation of 2700 masl. Elevation 2700 is the DMLZ cave propagation area based on TDR data (Fig. 2). The results of the CRT test show that the area around the abutment area and the seismogenic zone, especially areas experiencing cave propagation, can be recovered using both P and S waves. The CRT test is in line with the DWS value, bright (white) to red hues representing high DWS values for P and S data, especially in the NW area. The non-resolvable area is in the SE area, and the elevation is above 2900 masl. This is due to: no seismometer stations in the SE area (off-facilities), and the uneven seismometer distribution above the elevation of 3200 masl (Fig. 11). CRT and DWS test results are used to evaluate the resolution of the tomographic inversion results. The Derivative Weight Sum (DWS) is a measure that estimates the total length of rays affecting a model parameter, thus providing a more suitable evaluation of ray coverage. By calculating the ray density through DWS, we can determine the weighted total length of rays passing through each node in the inversion grid. The well-resolved areas are defined by the area which the checkerboard pattern is retrieved by the inversion and have DWS of more than 50 (black dashed line). The interpretable areas of the DWS are shown in light to dark red color marked with black dashed line (Fig. 12c and d).

Fig. 12
figure 12

Checkerboard Resolution test and Derivative Weight Sum (DWS) horizontal section at elevation 2700 masl. a P-wave velocity CRT results, b S-wave velocity test results, c P-wave velocity DWS test results, d Vs-velocity DWS test results, brown color is cave

Full size image

4 Results

Time-lapse tomography is a method to obtain Vp, Vs, and Vp/Vs ratio values ​​in a certain time interval. In this case we need to assess the resolution of each subset as the distribution of seismic event and the number of event that are not evenly distributed over time. This condition might cause bias because the inversion method used is Iterative Least Square with damping. In this research, CRT and DWS tests were used to assess the distribution of seismic events for each subset. Each subset restored CRT and DWS zones are reliable, and, therefore, changes in the velocity structure results can be used to analyze the changes in the stress distribution due to mining activities. Bright images that are specific for each subset of the data are displayed in places where the resolution test results indicate a high level of certainty.

Time-lapse tomography is shown by looking at changes in velocity evolution data from subset 1 to subset 4. Changes in the mechanical properties of the rock following the function of time are analyzed to describe the cave propagation with the provisions of the progress of the boundary cave analyzed through the distribution of magnitude values ​​and seismic event clusters in the seismogenic zone area. The cave propagation model is validated with Time Domain Reflectometry data. The velocity distribution resulting from tomography is used to describe changes in the condition of the rock mass around the cave validated by a conceptual model (Brown 2002; Duplancic and Brady 1999; Hidayat et al. 2022).

The cave propagation moves for around 20 m towards the southeast. At the same time, the highest movement is seen in the TDR data of subset 3 toward subset 4 in the NW, with the cave’s progression reaching 60 m at an altitude of 2700 masl (Fig. 2). We examine each data subset seismicity distribution to determine the cave inferred from velocity distribution evolution. Low velocity anomaly variations around the cave area periodically characterize the evolution of Vp, Vs, and Vp/Vs in the subsurface structure.

Fig. 13
figure 13

Vp tomogram of horizontal sections at an elevation of 2700 masl overlaid with seismicity (white to black gradation dots) and cave (brown color) for each subset. The moment magnitude is − 1.6 to 0.7 around the cave area. a is a tomogram Vp data subset 1, b is a tomogram tomogram Vp data subset 2, c is a tomogram tomogram Vp data subset 3, d is a tomogram tomogram Vp data subset 4

Full size image

Figure 13 is a tomogram of P-wave velocity combined with seismicity data, CRT and DWS resolution tests of each subset where interpretable areas are shown with light color shades while dark colors have low confidence. The tomogram of subset 1 (Fig. 13a) shows a cave front area with a low velocity anomaly surrounded by a high velocity anomaly (H1). In subset 2 (Fig. 13b), the high velocity anomaly zone H2 expands, contrasting with L2 which becomes more localized. In subset 3, the L2 zone once again expands in a SW to NE direction. The decrease in the low anomaly value at L3 is offset by the presence of a high velocity anomaly H3. The decline in velocity in L3 before collapse suggests extensive fracturing within that zone (Fig. 13c). This reduction in velocity in L3 impacts the increase of velocity observed in the H3 zone. Following the failure/cave growth in subset 3, the L4 anomaly value gradually reverts to its initial condition (Fig. 13d), accompanied by a gradual decrease in velocity within the H4 zone. The observed decrease in velocity is interpreted as a gradual release of stress. Velocity values that exceed the initial velocity model suggest induced stress due to cave propagation. This is consistent with the caving conceptual model that the seismogenic zone adjacent to the fracture area is an area of high stress (Cumming-Potvin et al. 2018). Rocks are stressed and experience increased velocity due to crack closure and elastic deformation. This continues with the growth of the crack, which is indicated by a decrease in velocity (Chen and Huang 2020; Gondim and Haach 2021).

The distribution of seismicity is distributed around the NW–SE area at various elevations with event magnitudes from − 1.6 to 2.0 Mw. Cluster events occur around the cave area, these events arise due to fractures that occur around the cave area. The event range between 0 and 0.7 Mw (white circle) around the cave area is identified as a seismogenic zone which indicates cave propagation. Microseismic clusters ≤ 0.7 trending SE–NW are marked with a white to black gradation dotted with a low velocity anomaly surrounded by high velocity. This low velocity anomaly is suspected as a stress release zone. The number of microseismic events in subset 2 decreased when compared to subset 1, in subset 3 it increased by 28.5% compared to subset 1 (Fig. 13). The increase in the number of events was correlated with cave propagation activity in the research area. The Vs wave tomogram shows the same thing as the Vp tomogram, on the Vs tomogram it is clear that the low anomalous movement in front of the cave continues until there is a change in the shape of the cave boundaries. fractures in rocks result in a decrease in the velocity value. This change is very clearly observed in subset 2 to subset 3. In this subset there are more fractures. Therefore, the velocity value is lower.

Fig. 14
figure 14

Tomogram Vs horizontal incision at an elevation of 2700 masl overlaid with seismicity (white to black gradation dots) and cave (brown color) for each subset. The moment magnitude is − 1.6 to 0.7 around the cave area. a is a tomogram Vs data subset 1, b is a tomogram Vs data subset 2, c is a tomogram Vs data subset 3, d is a tomogram Vs data subset 4

Full size image

Figure 14 represents a horizontal section at an elevation of 2700 masl obtained through inversion tomography Vs subset 1 to subset 4. In Fig. 14a (subset 1), the cave’s entrance is characterized by a low velocity anomaly. The cluster event at location 1 caused a reduction in velocity within that specific area. Figure 14b (subset 2) shows a distinct decrease in velocity (referred to as zone L2), which is accompanied by an increase in microseismic activity. Figure 14c (subset 3), it is intriguing to note that there is a significant decrease in S-wave velocity before the cave propagation process. The interesting part is the observation of an increase in high velocity anomaly ahead of the cave area prior to the collapse (H3). This increase in high velocity occurs due to the closing phase of the crack in the rock caused by pressure.

Fig. 15
figure 15

Tomogram of Vp/Vs horizontal sections at an elevation of 2700 masl overlaid with seismicity (white to black gradation dots) and cave (brown color) for each subset. The moment magnitude is in the range of − 1.6 to 0.7 around the cave area. a is a tomogram of the Vp/Vs ratio of data subset 1, b is a tomogram of the ratio of Vp/Vs of data of subset 2, c is a tomogram of the ratio of Vp/Vs of data of subset 3, d is a tomogram of the ratio of Vp/Vs of data of subset 4

Full size image

Figure 15 presents a tomographic representation of the Vp/Vs ratio for subsets 1 to 4. As the velocity tomogram, the Vp/Vs tomogram illustrates variations in low and high velocity anomalies at certain locations. The inversion results of the Vp/Vs ratio also provide insights into the development of caves and seismogenic zones in underground mines. The seismicity clusters with low velocity anomalies can be observed as indicators of seismogenic zones. Based on the analysis of velocity anomalies in different zones, it is observed that subset 1 of the L1 zone near the cave area exhibits low velocity anomalies, while the H1 zone shows high velocity (Fig. 15a). The high velocity anomaly in the H2 zone is increasing and extending towards the front area of the cave, whereas there are localized low velocity anomalies in subset 2 at the front area of ​​the cave (Fig. 15b). In subset 3, before the rock mass fracturing due to cave propagation, the Vp/Vs ratio value in the L3 zone drops again and tends to expand, which is limited by the high Vp/Vs value of H3 (Fig. 15c). After the cave growth, the Vp/Vs ratio values in zones L4 and H4 gradually return to the initial velocity model (Fig. 15d). The Vp/Vs ratio tomogram indicates changes in low and high velocity anomalies in various locations, similar to the Vp, Vs tomogram. The interesting part about the Vp, Vs and Vp/Vs ratio tomograms is that microseismic event clusters occur between low and high contrast anomaly. The Vp/Vs ratio inversion results can describe cave propagation and seismogenic zones of the underground mines. The seismogenic zone is visible from the seismicity cluster with low velocity anomalies.

5 Discussion

The cave propagation is indicated by low velocity evolutionary movements. The direction of cave propagation can be seen from the evolution of the cave and low velocities surrounded by high velocities. This low velocity value observed at an area of around 60 m from the cave in the NW area (Figs. 13 and 14, and 15) is interpreted as a compensation for decreased velocity due to the low velocity anomalies interpreted as pressure release, rock breaking and cracking caused by caving activities. This anomaly of cave propagation is observed especially at 2700 masl (100 m above the undercut level). Interestingly, there is a high velocity anomaly surrounding the low velocity. We interpret this zone as a consequence of the stress redistribution due to compression from the cave direction. Changes in low velocity values ​​can be seen from subsets 1–3 where in subset 4 the low velocity anomaly disappears, and the velocity returns to its initial state.

In this study, the stress distribution refers to the induced stress perturbation in underground mining based on research conducted by Wang et al. (2018) and Ma et al. (2019). Wang et al. (2018) conducted P-wave velocity tomography imaging in the underground mining area at Yongshaba Deposit, Guizhou Province, China. From the results of this study, there are low velocity anomalies interpreted as void space, pressure release, rock breaking and cracking caused by mining activities. In addition, there is a high velocity anomaly surrounding the low velocity, which is interpreted as a consequence of pressure redistribution due to mining activities. Similar results were also obtained from a study conducted by Ma et al. (2019) in an underground nickel mining area in Ontario.

Based on the tomogram interpretation of the 4 data subsets, it appears that the mechanism of cave propagation can be explained by changes in the distribution of velocity and seismicity. Velocity increases first due to the redistribution of stress due to cave excavation induced. Then there was a decrease in velocity and an increase in the number of seismic events (0–0.7 Mw). Within two weeks (one-time subset) after the velocity decreased, cave propagation was observed from the TDR data. After that, the velocity increases again due to the redistribution of stress. The cycle continues to repeat each caving process. So, it can be concluded from the research results that the passive seismic tomography method is also able to explain the phenomenon of velocity changes due to cave propagation and stress redistribution as well as the active seismic tomography method (Poulter et al. 2022).

The evolution of Vp, Vs and Vp/Vs ratio velocity observed throughout this study serve as indicators of the geomechanical phenomena occurring within the surrounding rock mass of the cave. The region in the NW side of the cave has high velocity changes during the observed times. We interpreted this area is highly affected by caving due to stress redistribution resulting from mining progress. Assuming areas with low velocities as stress release zones or damaged zones and those with high velocities as high stress zones (as ilustrated in Figs. 13, 14 and 15), we can observe how stress evolves within the rock mass prior to caving.

  1. a.

    During the subset 1, stress is redistributed within 150–200 m in front of NW cave boundary (indicated by the dotted white line). This perturbed stress area is evident trough the presence of both low and high velocity zones that are scattered throughout.

  2. b.

    During subset 2, there is clear separation between low and high velocity zone in the vicinity of the NW area of the cave. The concentration of the low velocity zone becomes more pronounced as the cave propagates. In front of the low velocity area, there is accumulation of higher stress levels in the rock mass induced by the damaged area.

  3. c.

    The stress concentration persisted until subset 3, where it was seen that low velocity zone in front of the cave expanded and became more concentrated. The significance of reduction in rock velocity suggests that the area is experiencing fracturing as a result of stress concentration surpassing the rock mass’s yield limit.

  4. d.

    This is supported by the phenomenon of caving in subset 4, and the redistribution of stress around the new cave boundary area occurs one more, as indicated by the high velocity.

Based on the comparison of Figs. 13 and 14, and 15, it can be inferred that the changes in Vp have a greater impact on the stress evolution compared in Vs. This is because P-waves are more sensitive to fracture propagation. However, upon closer examination, it can be observed that areas with low Vs values and low Vp/Vs ratio play a crucial role in determining the boundaries where collapsed planes form between unstable regions and stable rock masses during caving. This research also shows the notable results that Vs change follows the cave propagation cycle. The decrease of Vs occurred because of rock mass fracturing. Approaching cave propagation, the Vp/Vs ratio decreased very sharply because the decrease in Vs became more significant following uncontrolled fracture propagation that occurred before the rock mass was caved/failed. Thus, monitoring the Vp/Vs ratio is highly recommended as an indicator of the cave propagation cycle.

It is worth to note that the information obtained from monitoring results using tomographic analysis or monitoring using Vp, Vs, and Vp/Vs analysis takes a relatively longer time compared to the daily mining production activities. In general, for monitoring the cave growth progress, 4D tomography could provide sufficient information. However, if there is an urgent warning information related to the underground mining activities, there might be a delay due to the time required for data collection and processing in tomography. It is because tomography relies on sufficient data coverage and distribution, which affects its resolution and technical requirements. Recently, there is a development of monitoring velocity variation in a shorter period of time by using ambient noise tomography data through cross correlation at each station point ( Olivier et al. 2015; Olivier and Brenguier 2016). The shorter velocity variations from ambient noise data could be integrated with tomography body wave tomography to better monitor the rock mass changes in underground mining.

6 Conclusion

Our tomographic method has been successful in monitoring cave propagation over time. Underground mining activities are very dynamic and can change the velocity structure in the short time. From the imaging results of Vp, Vs, and Vp/Vs, there is a pattern that always changes from time to time. Cluster events in the seismogenic zone area are indicators of cave propagation. In this research, the microseismic events associated with rock mass fracturing are in the range of values ​​between 0.0 and 0.7 Mw. The event has a positive effect on cave propagation by increasing cave propagation. The low velocity anomaly indicates the presence of rock mass fracturing in the area around the cave. Meanwhile, the high-velocity anomaly is compensation due to increased pressure that is resulted from stress redistribution due to cave growth. Its optimal time range per subset follows the mining rate and uses the ratio Vp/Vs as an indicator of cave growth. Furthermore, the available Time Domain Reflectometry (TDR) data has been used to validate the interpretation of the cave growth based on the velocity tomogram. The results of this research are expected to provide the essential contribution to seismic for monitoring in underground mines, especially the cave mining method. In maximizing a microseismic monitoring, especially in areas with no/minimal ray path, we suggest the implementation of even distribution for seismometer stations off-facilities areas.