Potvin illustrates the stability graph method based on the calculation of the hydraulic radius HR in the critical surface and the modified stability number N′. These two factors were plotted on the graph which is divided into the stable zone, transition zone, and the caving zone, so as to express the stability of the underground excavation area (Potvin 1988). The hydraulic radius HR in the stability graph can be obtained using Eq. 1: H R = length × width 2 × length + 2 × width (1)

The stability number N′ (Potvin 1988) is calculated by Eq. 2: N ′ = Q ′ × A × B × C (2) where Q ’ can be obtained by using Eq. 3. Q ′ = R Q D J n J r J a (3) where RQD is the rock quality index, J n is the joint number, J r is the joint roughness number, and J a is the joint change number.

Parameter B in Eq. 1 is the joint adjustment factor, and its value is measured by the difference between the dip angle of the mining face and the main joint set, which is shown in Figure 2A. Parameter C is the gravity adjustment factor, which is related to the failure mode of the goaf under gravity, which is shown in Figure 2B below.

A is the rock stress factor, which is calculated by the ratio of the uniaxial compressive strength (UCS) in intact rock to the compressive stress σ _1max is the maximum induced tangential stress acting at the center of the critical face. The value of A is linearly related to UCS/σ _1max, and its range is 0.1–1. It is calculated by the following Eqs. 4 and 5: R = U C S σ 1 max , (4) A = { 0.1 , I f   R < 2 0.9 8 ( R − 2 ) + 0.1 , I f   2 ≤ R ≤ 10 1 , I f   10 < R . (5)

It can be seen from Eq. 5 that the Potvin modified graph method believes that the stope will be in a stable state under low-stress state (when 10 < R). At this time, the obtained A is the maximum, which leads the calculated stability number N’ to the maximum value. By accounting for the adverse relaxation effect in the stope under low stress (10 < R), Mitri et al. (2011) revised the calculation method for A by introducing the maximum stress factor (MSF) at the center of the critical face: M S F = 1 R = σ 1 max U C S , (6) A ′ = { 0.1 , I f   M S F < 0 0.1 + 0.9 e [ ( − e − ( M S F − 0.3 0.09 ) ) − M S F − 0.3 0.09 + 1.0 ] , I f   0 ≤ M S F ≤ 1 0.1 , I f   1 < M S F . (7)

In the widely used (Potvin (1988) and Mitri et al. (2011)) modified graph method, Eqs 4, 5 of stress factor A and A ′ are composed of UCS and the maximum principal stress of stope surface σ _1max calculated from numerical simulation. This calculation method has the following deficiencies:

1) The stress state at any point in the plane is determined by the maximum principal stress σ 1 , the intermediate principal stress σ ₂ and the minimum principal stress σ 3 while the original stress factor A and A ′ only consider the maximum principal stress σ 1 of the stope surface, which ignores the influence of minimum principal stress σ 3 and the intermediate principal stress σ 2 on the stress factor, Moreover, it do not fully consider the failure of the stope caused by tensile stress.

According to the aforementioned two problems, this article uses the modified Hoek–Brown criterion with tensile cut-off (Dai et al., 2018) as the judgment basis for rock mass engineering instability. The Hoek–Brown criterion is demonstrated by Hoek and Carranza-Torres, 2002, who statistically analyze a larger number of rock-triaxial test data and rock field test results, comprehensively considering the influence of the rock-mass structure, rock strength, stress state, and other aspects to obtain the empirical formula of the rock mass failure, which can reflect the inherent nonlinear failure characteristics of rock mass. Therefore, this criterion is widely used in underground chambers, slopes, tunnels, and other engineering projects in the rock engineering sector. Hoek–Brown criterion with tensile cut-off is shown as follows (Dai et al., 2018): F 1 ( σ ) = σ 1 − σ 3 − σ c i ( m b σ 3 σ c i + s ) a = 0 , (8) F 2 ( σ ) = I 1 − σ t 3 = 0 , (9) m b = m i   exp ( G S I − 100 28 − 14 D ) , (10) s = exp ( G S I − 100 9 − 3 D ) , (11) a = 1 2 + 1 6 ( e − G S I 15 − e − 20 3 ) , (12) where σ ₁ is the maximum principal stress of rock mass failure, σ 3 is the minimum principal stress of rock mass failure, σ c i is the uniaxial compressive strength of intact rock, m _i is the constant factor of rock material, which can be obtained by fitting the triaxial test data of rock, or determined by rock type, D is the perturbation factor; it equals to 0 when it is below 100 m, GSI is the geological strength index, and I ₁ is the first stress invariant.

Based on the MSF calculation method proposed by Mitri et al. (2011) and the Hoek–Brown criterion with tensile cut-off (Dai et al., 2018), the MSF _c calculation under shear failure mode is modified according to Eq. 13, and the MSF _t calculation under tensile failure mode is added according to Eq. 14: M S F c = | σ 1 - σ 3 - σ c i [ ( m b σ 3 σ c i + s ) a - s a ] σ c | , (13) M S F t = | I 1 = σ 1 + σ 2 + σ 3 σ t | . (14)

In Eqs 13 and 14, MSF _c is the maximum stress factor upon shear failure in the critical face of the stope, MSF _t is the maximum stress factor upon tensile failure in the critical face of the stope, σ 1 represents the maximum principal stress at the central rock mass in the critical face; σ 2 is the intermediate principal stress at the central rock mass in the critical face, σ 3 is the minimum principal stress at the central rock mass in the critical face, σ c is the uniaxial compressive strength of the rock mass, σ t is the tensile strength of the rock mass, and σ c i is the uniaxial compressive strength of intact rock. In this article, the tensile stress symbol is negative whilst the compressive stress symbol is positive.

When σ 1 ≤ σ 3 + σ c i ( m b σ 3 σ c i + s ) a , the rock stress factor A′ can be calculated using Eq. 15. A ′ ¯ = { 0.1 , M S F t ≥ 1 1.1 − M S F t , 0.1 < M S F t < 1   and   σ 3 < 0 1     M S F t ≤ 0.1   or   σ 3 ≥ 0 . (15)

When σ 1 ≥ σ 3 + σ c i ( m b σ 3 σ c i + s ) a and I 1 = σ 1 + σ 2 + σ 3 < σ t , the rock stress factor A ′ can be calculated using Eq. 16. A ′ ¯ = { 0.1 , M S F c ≥ 1   or   | I 1 | ≥ σ t 0.1 + 0.9 ⁡ exp [ − exp ( M S F c − 0.3 0.09 ) − M S F c − 0.3 0.09 + 1 ] M S F c ≤ 1 , (16) where σ 1 ≥ σ 3 + σ c i ( m b σ 3 σ c i + s ) a and | I 1 | ≥ σ t in the first equation above demonstrate that the rock mass attains shear failure and tensile failure at the same time. However, it is considered as the tensile failure because the tensile strength of rock mass is far smaller than the compressive strength.

It is worth to explain that the tensile strength of the rock mass is much smaller than the compressive strength. For this reason, this article gives priority to judge the tensile failure when solving the stress factor. By distinguishing F 1 ( σ ) < 0 or F 1 ( σ ) > 0 and | I 1 | ≥ σ t in Eq. 8, it enters the calculation process of tensile failure. At this time, the central stress factor of stope critical face is calculated by the first expression of Eq. 16 or Eq. 15. M S F t is determined according to the first stress invariant I ₁ of the central rock mass of stope critical face and the tensile strength σ _t of rock mass. 1) When I 1 ≥ σ t , the tensile cut-off is the failure zone (3), and the stope is under extreme instability; therefore, A ¯ ′ is equal to the minimal value 0.1. 2) When I 1 < σ t , and if M S F t < 0.1 , it indicates that the tensile stress on the stope is far less than the tensile strength. At this stage, A ¯ ′ is equal to the maximum value 1. If 0.1 < M S F t < 1 , define A ¯ ′ = 1.1 − M S F t . It is particularly emphasized that when 0.1 < M S F t < 1 , A ¯ ′ is assumed to be linear, which is similar to the definition A by Potvin (1988) when dealing with high pressure stress.

When F 1 ( σ ) > 0 and I 1 < σ t in the discriminant Eq. 8, the stope face has a significance to experience compressive shear failure, which is the failure zone (1). The maximum stress factor M S F c in the critical face of the stope is calculated by Eq. 13, where M S F c is calculated from the plane point stress state and the compressive strength of the rock mass defined by the Hoek–Brown criterion with tensile cut-off (Dai et al., 2018). The calculation method of A ¯ ′ is based on the M S F calculation method demonstrated by Mitri (Mitri et al., 2011).

In order to obtain the evolution characteristics of stress and displacement of the roof in the northern goaf of Xiaobeigou, combined with the geological data in the northern goaf of Xiaobeigou and the thickness distribution of the roof in this area, a group of four-point bolt dynamometers and two-point multipoint displacement meters are arranged in each area of A and B in the northern goaf for real-time monitoring of stress and displacement. The sensor measurement method is 24-h automatic uninterrupted monitoring. The layout plane position of the measuring points is shown in Figure 3. The green measuring point in the figure represents the layout position of the multipoint displacement meter, whilst the red measuring point represents the layout position of the bolt dynamometer.

The depth of holes A and B for the bolt dynamometer is about 21 m. The first force sensor is about 7.7 and 7.4 m from the vertical distance to the goaf roof, respectively. Meanwhile, each force sensor is arranged every 3.5 m to complete the layout of the two sets of bolt dynamometers. The depth of holes A and B of the multipoint displacement meter is about 21 m, and the first measuring point is arranged at the bottom of the hole. The first measuring point of hole A is about 7.2 m away from the roof of the goaf, and the second measuring point is about 11.2 m away from the roof of it. The first measuring point of hole B is about 3.4 m away from the roof of the goaf, and the second measuring point is about 7.4 m away from the goaf roof.

According to the field anchor dynamometer and preliminary geological survey results, there is no significant tectonic stress in the field, and the buried depth of the goaf is shallow, so its stability is mainly affected by different lithology, lithology interface, and gravity. Therefore, the initial stress field calculated by the model is constructed according to the self-weight stress of rock mass, and the influence of two kinds of the ore-rock interface is considered, which is shown in Figure 7. The constitutive model used in the model is the mesoscopic elastic creep damage constitutive model, which is referenced by Zhu and Tang, (2004). The model was discretized into 140,057 tetrahedral elements. The representative volume element (RVE) is used in this study to improve computational efficiency. The RVE should be small enough (relative to the macroscopic rock mass) in order to characterize the rock heterogeneity with sufficient precision, but it must be large enough so that it contains enough mesostructural information to be regarded as a material point on a macroscopic scale. According to our previous experience with the numerical simulation at this engineering scale based on the Monte Carlo method, a minimum element size of 0.5 m is selected. The bottom surround boundary of the model is constrained by rollers. In order to avoid the boundary condition influence of the calculation model on the goaf deformation in the calculation process, the boundary sizes of the model are 800 and 600 m, respectively. The degree of surface fluctuation and slope is established by 3Dmine-Rhino, which is achieved by inputting the surface contour elevation data. The results are shown in Figure 7. It can be seen that the interface of different lithology passes through the goaf. The shape of the goaf is determined according to the CMS of the goaf scanning detector. The parameters used in the numerical simulation are shown in Table 1.

Due to historical reasons, the goaf is formed by irregular mining and disorderly excavation. It is impossible to understand when the excavation is carried out, which means the excavation time is unknown, it brings difficulties to the numerical simulation analysis of the creep damage process of surrounding rock after excavation, and it is difficult to establish the starting point of time-dependent deformation under the influence of excavation. In order to study the current stability of the goaf area, and predict the failure time of the goaf area in the future according to the current stability, the simulation method is established as follows:

(1) Firstly, the creep damage process of roof position A and position B (as shown in Figure 8) after goaf excavation is analyzed by numerical simulation, which obtained the creep deformation time curve of the surrounding rock mass. The mechanical parameters of surrounding rock mass are obtained according to the field investigation, and they are summarized in Table 1. Therefore, the main purpose of this step is to reproduce the creep deformation rate of surrounding rock mass at the present stage by numerical simulation combined with field monitoring, without paying attention to the starting point of creep deformation of surrounding rock mass.

The comparison results are shown in Figures 9A, B. It can be seen that the creep deformation rate of the surrounding rock obtained by the numerical simulation and field monitoring is highly consistent in the period of about 20–130 days, indicating that the numerical simulation can well simulate the deformation changing of the goaf roof with time in the field at this stage, and it is able to predict the development of the deformation of the goaf roof in the future and reasonably evaluate the stability of it.

It can be seen from Figures 9A, B that the creep deformation rate at monitoring point B of goaf roof is larger than that at monitoring point A. The vertical creep deformation of point A is less than 1 mm within 100 days, whilst the point B is less than 2 mm within the same period. The vertical displacement of the two points generally sink to the goaf. The numerical simulation results are in good agreement with the field monitoring results. This shows that the numerical simulation method in this article can simulate and predict the timeliness deformation law of goaf roof at the same time.

Based on the simulation method determined in the above section, the development of roof deformation in the southern goaf is analyzed, and the deformation in the next five years will be simulated and predicted. The simulation results are shown in Figure 10. The period of region I in the figure can be considered as an unknown time-dependent deformation caused by historical reasons. The period of region II is the verification period based on the monitoring data. Area III period is predicted by numerical simulation based on area II period. It can be seen that, at monitoring point C, due to the influence of the lithologic interface, the deformation decreases in the vertical position with the extension of time, which means the point moves slowly to the surface. The deformation rate of monitoring point E is the largest. Compared with the monitoring period, the creep deformation increment could reach over 13 mm in 5 years, which is the position where the deformation of the goaf roof is large. This is consistent with the evolution trend of damage in the goaf area obtained in Figure 11. It can be seen that, under the influence of excavation disturbance, there are small-scale damaged areas on the roof of the southern goaf, but the damaged areas are not connected and scattered. With the development of creep deformation, the damaged areas of the roof of the southern goaf gradually evolve and develop, and local damaged areas appear on the roof and floor of the goaf. When the creep deformation time reaches 5 years, the scattered damaged area of the southern goaf roof forms a connected area, and the damage mode of caving arch is formed at the goaf roof. This shows that the creep deformation in the future is more serious than that in the southern goaf, and the damaged is easier to penetrate to the surface. However, the damaged area in the north of goaf expanded slowly and remained stable.

The northern goaf is about 65 m long, 48 m wide, and about 6 m high. The southern goaf is about 130 m in length, 40 m in width, and 5 m in height. According to 3DMine modeling software, the roof contour perimeter and area of the north and south goaf in Xiaobeigou can be obtained, and the hydraulic radius HR of the north and south goaf can be calculated, as shown in Table 2.

In this article, the structural plane investigation of +145m platform rock mass is carried out by ShapeMetrix3D digital measurement system, and 15 photogrammetry points are selected. The rock mass structural plane is synthesized by the obtained image according to the pixel matching of the same point and the visual distance judgment technology. The average value of the digital measurement results is taken to obtain the information on rock mass joints, cracks, and geometry. The basic rock mass quality index BQ value of the northern goaf is 512.72, whilst the southern goaf is 471.04. According to Eqs 17 and 18, the BQ value can be converted to Barton’ s rock mass quality Q value, and the corresponding Q’ value can be obtained (Liu et al., 2017). RMR = 15ln Q + 50 , (17) RMR  =   ( B Q - 80 . 79 / 6 . 09 ) . (18)

According to the numerical calculation model, the excavation elasticity of the goaf area is calculated, and σ 1 , σ 2 , and σ 3 of the middle line in the north and south goaf area are determined by the numerical simulation results. Then, the rock stress factor A is calculated according to the correction formula from the above section, and the calculation results of the relevant parameters stability number N are shown in Table 3.

Figure 12 shows Mathews stability graph calculation results of the north–south goaf in Xiaobeigou. It can be seen from the figure that the roof of the northern goaf is located in the transition zone, which is a relatively stable state, while the roof of the southern goaf is in the caving zone and has a certain caving risk, which is consistent with the numerical calculation results in Section 5.3. In the calculation process of Mathews stability graph, the rock stress factor A ¯ ′ needs to be obtained by numerical simulation, which has a high requirement for the accuracy of the numerical model. In this article, the potential risk area is determined by the previous numerical calculation, and then the stress and displacement monitoring points are reasonably arranged. Based on the actual site monitoring data, the numerical model is dynamically fed back to correct. The accurate calculation of the numerical model is fed back to the theoretical calculation again to reach the theoretical value. Field monitoring and numerical calculation are independent and mutually supportive, which are applied to achieve the joint evaluation of the stability of Xiaobeigou goaf and the prediction of the future trend.