Abstract

Although the impact of Karst Collapse Pillars (KCPs) on water inrush has been widely recognized and studied, few have investigated the fluid-solid interaction, the particles migration inside KCPs, and the evolution feature of water inrush channels. Moreover, an effective approach to reliably predict the water inrush time has yet to be developed. In this work, a suite of fully coupled governing equations considering the processes of water flow, fracture erosion, and the change of rock permeability due to erosion were presented. The inverse velocity theory was then introduced to predict the water inrush time under different geological and flow conditions. The impact of four different controlling factors on the fracture geometry change, water flow, and inrush time was discussed in detail. The results showed that the inverse velocity theory was capable of predicting the occurrences of water inrush under different conditions, and the time of water inrush had a power relationship with the rock heterogeneity, water pressure, and initial particle concentration and an exponential relationship with the initial fracture apertures. The general approach developed in this work can be extended to other engineering applications such as the tunneling and tailing dam erosion.

1. Introduction

Groundwater inrush in coal mines has caused thousands of fatalities across many countries, such as the USA, Russia, Poland, Canada, Australia, Germany, Great Britain, India, and especially China [1]. As mining activities progress further deeper underground, there have been an increasing number of water outburst incidents occurring annually largely due to the impact of some key geological structures such as Karst Collapse Pillars (KCPs), as illustrated in Figure 1 [2, 3]. As a special geological structure, KCPs commonly exist at more than twenty coal fields in the northern China, such as Shanxi province (see Figure 2) [4].

The existence of the geological structure usually functions as an underground water flow path, thus posing a great threat to the mining production safety [3, 5]. According to the statistics, the most serious water disasters in Chinese coal mines were all directly caused by water inrushes associated with the KCPs under the Permo-Carboniferous coal seams located at northern China [6, 7]. Take Fangezhuang coal mine of Kuiluan Group, for instance, during the accident in 1984, the maximum water inflow associated with KCPs reached 2,053 m3/min, resulting in the mine and other three neighboring mines submerged in a very short time, and caused a direct capital loss of around US90 million [8]. Thus, better understanding of the mechanism of water inrush has both economic and safety benefits.

Significant number of investigations have been performed to explore water inrush mechanism of KCPs using single or combined methods of theoretical analysis and numerical simulation as well as experimental studies. For instance, based on the elastic thick plate theory, Tang et al. [4] developed a mechanical model for water inrush of KCPs. Bai et al. [9] established a mechanical model-plug model, which was used to describe the behavior of water seepage flow in coal-seam-floor containing KCPs. Furthermore, the variable mass dynamics and nonlinear dynamics were introduced, and the seepage properties of KCPs associated with particles migration were investigated using numerical simulation [3]. Ma and Bai [10] numerically studied the impacts of mining-induced damage on KCPs and the surrounding rocks and on the formation of the fracture zone and analyzed mining-induced KCP groundwater inrush risk. This work was later extended to study the effect of coal mining operations on KCPs related groundwater inrush using , and the shear stress, damage zone, and their effects on seepage field development were also discussed [11]. Yao [12] experimentally studied the evolution of the crushed rock mass seepage properties under different particle sizes and stresses and analyzed particle migration feature and the KCPs water inrush mechanism. Moreover, there are also some experimental studies focused on the permeability change of the KCPs to investigate water inrush mechanisms [1315].

Recent studies on the water inrush mechanism of KCPs have consistently revealed that complex interactions exist between the evolution of fractures and porosity in solids, water transport, and rock solid particles migration, which accompanies the erosion phenomenon [3, 12]. Significant progress on the understanding of this complex process has been made. For example, the reaction of aqueous solutions mineral components and its impact on water flow were investigated, and a model for porosity change was established [16]. Vardoulakis et al. [17] studied the piping and surface erosion effects based on mass balance, particles migration, and Darcy’s law. The dissolution of fracture surfaces due to water-rock interactions with variable fracture apertures over time was also studied, and a depth-averaged model of fracture flow and reactive transport that explicitly calculated local dissolution-induced alterations in fracture apertures were presented [18]. Habib et al. [19] numerically investigated the erosion rate correlations of a pipe protruding from an abrupt pipe contraction problem.

Despite the significant progress in understanding KCPs water inrush mechanism and water prevention technology over the past several decades [2022], most investigations were carried out with respect to structure failure, and few have investigated this issue by integrating the solid-fluid interaction and particles migration due to erosion. Water inflow rate, as probably the only direct measuring data in the field, has yet to be directly used in the literature to predict the precise time of water inrush. An effective method to predict water inrush does not appear to be available neither. In this work, a set of fully coupled governing equations for KCPs water inrush will be developed by incorporating the fractures and porosity evolution for KCPs, particles migration, and water seepage process. The governing equations will then be implemented into COMSOL Multiphysics software. The heterogeneity of rocks in KCPs will be determined using the Weibull distribution, and the parameters including porosity, seepage, particle concentration, and water inflow as well as seepage channels evolution law will be obtained accordingly. This numerical model will be then used to investigate the water flow characteristics with erosion. Based on the inverse velocity theory, the corresponding water inrush times for different flow conditions will also be discussed in detail.

2. Development of a Fully Coupled Theoretical Model

2.1. Selection of Microstructure Model of KCPs

A KCP typically consists of three parts: solid rock matrix, fluids, and infilling particles. The intact rock blocks form a rock mass, which is filled with fluids, for example, water, and infilling materials such as solid particles (Figure 3). As the permeability of rock blocks is generally very low, the interrock fractures provide the main passages for water flow. This structure matches the characteristics of the dual porosity model, where the rocks are considered as consisting of matrixes and fractures. A numbers of different dual porosity models have been developed including the Warren-Root, Kazemi et al., and de Swaan models [2426], but the Warren-Root model, as shown in Figure 4, was widely accepted as a proper model to represent such rock structure. As such, in this work, the Warren-Root model is selected to examine the processes of particle migration and the evolution of the fractures aperture and porosity.

2.2. Governing Equations
2.2.1. Assumptions

In order to establish a fully coupled theoretical model that integrates particles migration, seepage, and fracture and porosity evolution in KCPs under erosion effects, the following assumptions have been made:(1)The fluid and particles in KCPs are incompressible.(2)The suspended particles share the velocity field with the fluid.(3)As the internal structure of KCPs is generally very loose and the rock mass surrounding the KCPs is generally competent, the impact of the change in the effective stress on KCP’s permeability is insignificant.(4)The impact of erosion on the permeability change is in proportion to the change of particle concentration in the fluid.

2.2.2. Definition

Based on the aforementioned assumptions, a representative element for a typical KCP microstructure is shown in Figure 5. The length (m) and volume (m3) of the element are expressed aswhere (m) is the length of the matrix and (m) is the aperture of the fracture.

The total voidage (%) of an element is written as [28]where (%) is the matrix porosity. It should be noted that .

2.2.3. Mass Conservation Equations for Particles

In order to study the migration characteristics of the infilling particles, the dynamic mass conservation equations for particles were calculated based on a three-dimensional representative element, as illustrated in Figure 6.

Infilling particles migration inside a representative element under erosion effects is affected by both convection and processes. The average seepage velocity of an element can be defined as (m/s). For the convection and diffusion process, in the direction of (), the mass losses of particles flowing out of the element in a unit time can be calculated by (3) and (4), respectively:where (%) is the particle concentration (i.e., particle saturation), (kg/m3) denotes the density of solid particles, and (m2/s) illustrates the particles diffusion coefficient.

In the direction of , the mass flowing into the element due to the convection-diffusion should be expressed as

For the element, the total mass flowing can be expressed as the summation of the three directions:

Under the effects of erosion, the mass loss of the element in a unit time can be summarized aswhere (kg/m3/s) is the particle mass that migrates into the fluid from a unit element in a unit time and is written as

According to the mass conservation law, for the unit element, (6) and (7) are equivalent, giving the mass conservation equation for particles as

Substituting (2) into (9) yields

As particle sizes are normally several orders of magnitude smaller than the fracture aperture, the diffusion effect can be neglected. Equation (10) can therefore be simplified into

2.2.4. Water Mass Conservation Equations

For the water flow, the mass flowing out of the element in the direction of () can be expressed aswhere (kg/m3) denotes the density of the fluid. The mass loss of the fluid in the element in a unit time can be calculated as

By combining (2), (3), (12), and (13), the mass conservation equation for water flow can be derived as

2.2.5. Evolution of the Fracture Aperture and Porosity

Rock porosity plays an important role in determining rock permeability, so full understanding of its evolution is essential to study the seepage characteristics of KCPs. Sakthivadivel and Irmay [29] investigated the erosion problem for porous media by using both experimental and theoretical methods. Vardoulakis et al. [17] summarized the studies on porous media erosion, and according to that study, the evolution of fracture aperture and permeability was affected by the porosity and particles concentration as well as the seepage velocity. The following equations have been developed to determine the evolution of fracture aperture and porosity [17]:where and are constant, is the absolute value of the seepage velocity (m/s), and and illustrate the maximum value that fracture aperture and porosity can reach under fluid erosion.

Equations (15) indicate that the porosity evolution is proportional to the particle concentration, as well as the seepage velocity.

2.2.6. Water Seepage

Fluid transport in porous media is normally described by Darcy’s law, which is derived (a) from balance of momentum for the fluid phase and (b) from a constitutive equation for the fluid-solid interaction force (i.e., the seepage force) [30]. In general, the balance of linear momentum for the fluid phase has the form [31]where denotes the pore pressure (Pa), is the dynamic viscosity of the fluid (Pa·s), and illustrates the permeability of the volumetric element (m2). For this study, the flow velocity is relatively low; thereby, the acceleration term can be neglected, givingwhere is the Darcy velocity (m/s) and is the unit vector in the direction of gravity. For fractured porous media, is the sum of fracture permeability and porosity permeability (m2), which can be expressed as [32]in which and represent the porosity permeability and fracture permeability, respectively.

Equations (2), (11), (14), (15), and (17) together with (18) compose the coupled processes of water transport of KCPs under erosion effects, as illustrated in Figure 7. The above governing equations will be implemented into COMSOL Multiphysics next to understand water inrush mechanism in KCPs, to evaluate water inrush risks, and to predict water inrush time under different rock and flow conditions. There are six unknowns with equations, so the set of governing equations can provide a unique solution.

In petroleum engineering, generally there is a critical velocity associated with sand production in oil and gas reservoirs. However, for the particle movements in KCPs, the critical velocity was not considered because the sizes of the infilling particles are significantly smaller than the fracture apertures (μm versus mm); the movement of the eroded materials does not require a critical pressure gradient (or flow velocity) for most KCP infilling materials are soluble to water forming diluted solution, and the fluid flow inside KCPs is generally much greater than the critical velocity of around 0–7.7 mm/s for the particle size of 1–10μm [33, 34].

2.3. Inverse Velocity Theory

The inverse velocity theory was originally proposed by Fukuzono [35] based on experimental studies. This theory was initially used to predict the failure time of a solid material that experiences slow continuous deformation (i.e., creep). When the reciprocal of the deformation velocity (i.e., the inverse velocity) is plotted as a function of time, its value approaches zero as the velocity increases asymptotically towards failure. The trend line of the inverse velocity intersects with the abscissa which represents time, and the intersectional point can predict the time of failure. There are three types of trend line, namely, concave, linear, and convex, as shown in Figure 8.

The following equation is used to define the envelope:where is the inverse velocity (m/s), and are constants, is the time of failure, and is the time of cutoff or the most recent time of the monitoring data. The inverse velocity curve is concave when , and it is linear when , while it is convex when . Fukuzono [35] indicated that a linear fit could provide a reasonable prediction of failure time.

The inverse velocity theory has been successfully used in soils, rocks, and other materials and has been proved to be reliable in estimating the failure time using the linear fitting curve [36, 37]. Therefore, in this work, this theory (i.e., inverse fluid flow velocity) is adopted as the first attempt to predict water inrush in KCPs.

In order to verify the applicability of theory to water inrush, the field monitoring data from Yao et al. [38] and Bai [39] was used, and the comparison of inrush time between real data and predicted time using inverse velocity theory was conducted. The results were plotted in Figure 9. It can be seen that the predicted inrush time matches the actual water inrush time very well, which proves that this theory can be used to predict water inrush for KCPs.

3. Numerical Model Implementation

3.1. Background of the Coal Mine

In this section, a case study was conducted based on the hydrogeological conditions of the Tuanbai Coal mine, which is affiliated with Huozhou Coal Electricity Group of China. The mine is located in Tuanbai County, Linfen City, Shanxi Province, as marked in blue in Figure 2. The mining operation is current extracting number 10 coal seam at approximately 300 m depth. To date, 97 KCPs have been found during the mining operations at this seam. These KCPs included rock blocks of different sizes, and in most cases, mud infillings and sizes of KCPs varied from tens to hundreds of meters. In addition, some had low level of compaction and cementation and relatively high infiltration capabilities.

Since the distance between number 10 coal seam and the Ordovician limestone was only about 35 m (as illustrated in Figure 13(a)), the mining operation is being carried out under pressure from the Ordovician limestone water of the floor. The water pressure was found to be about 2 MPa. This mine has suffered a number of water inrush accidents associated with the KCPs. For instance, one accident occurred on the belt roadway of the first longwall panel in 2007, with a water flow rate of nearly 470 m3/h. In addition, there were a number of water inrush accidents on longwall panels such as 10-106 as shown in Figure 10, 10-112 and 10-114 at a rate of range from 40 m3/h to 150 m3/h in 2010, resulting in the submergence of the longwall panels. Later investigation of these incidents revealed that all were directly caused by the KCPs, which acted as water channels, with the confined water serving as the source.

3.2. Characterisation of the Infilling Materials

In KCPs, the composition of the infillings and their particle size distribution (PSD) have been found to be highly relevant to the migration of the particles. For this reason, two samples (i.e., sample 1 and sample 2) were collected from the Tuanbai coal mine, and the X-ray diffraction phase analysis (composition analysis) was carried out on the infilling materials of the KCPs as well as the PSD characteristics. The results are shown in Figures 11 and 12.

The composition analysis (as shown in Figure 11) indicated that the infilling materials mainly included kaolinites, quartz, illite, smectite mixed layers, small quantities of feldspars, montmorillonites, chlorites, calcites, siderites, and other mineral compositions. The PSD results (Figure 12) showed that the particle sizes varied from 1μm to a few mm for both Samples 1 and 2, with majority spanning from 10μm to 1 mm. This could be a result of long-term weathering, but the weak and unfavorable rock properties observed in the field indicate that KCPs are very fractured and loose in structure, with low strength, high permeability, being prone to water erosion, and the water inrush risks.

3.3. Numerical Model Setup

Although diverse in geometry, KCPs are usually presented similarly to a dome shape [9]. Each of these approximates a 3D axisymmetric model, which can be simplified as a 2D symmetrical model.

In this work, according to the lithology of the mine, the bottom diameter, top diameter, and height are 15 m, 10 m, and 20 m, respectively (as shown in Figure 13). The pressure at the water inlet on the lower boundary was 2.0 MPa and the pressure at the water outlet on the top was 0.1 MPa; initial water pressure in the model was 0.1 MPa; initial particle volume concentration was 0.01; initial average fracture aperture was 0.1 mm; and initial porosity was 10%. The main parameters of the model were all determined according to the relevant references [3, 12], in which the studies were conducted on the same seam. Other used parameters are listed in Table 1. The model was divided into 5,000 grids using mapped mesh method embedded in the COMSOL software. The shape of the grids is rectangular.

Heterogeneity has been determined to be a prominent feature of rock. There are currently various approaches to obtaining the characteristics of the heterogeneous distribution of rock materials. For example, digital core techniques and mathematical statistics methods have both been utilized. Previous studies have shown that the heterogeneity of rock can be described by the Weibull distribution [4042]. Thus, the Weibull distribution was selected in this study due to its effectiveness and great simplicity to obtain the heterogeneous distribution of fracture apertures of the KCPs. The distribution probability density equation was as follows:where indicates the fracture aperture, denotes the initial fracture aperture, and is the uniformity index. The larger represents, the higher level of uniformity. The fracture aperture distribution obtained by the numerical generation method is shown in Figure 13(c).

Based on this model geometry, (9), (11), (14), (15), (17), and (18) will be solved simultaneously to investigate the dynamic change of water flow behaviors and the associated inrush risk under different flow conditions.

4. Results and Analysis

In order to analyze the characteristics of seepage at different times, in this study, six different moments were selected (i.e., 14 × 103, 15 × 103, 16 × 103, 17 × 103, 17.5 × 103, and 18 × 103 s). The specific results will be discussed below.

4.1. Distribution of the Fracture Apertures

Figure 14 illustrates the spatial changes of the fracture apertures at different times. Figure 15 shows the evolution of fracture apertures along O-O′′ and O′-A′ lines as marked in Figure 13(b), and the average fracture aperture-time curve was plotted in Figure 16. The results show the following:

The fracture aperture gradually increases under the effects of the erosion. However, the increase rate is higher in the upper part of the model than in the lower part. This is likely because the upper outlet boundary is narrower than the bottom inlet boundary, giving rise to higher flow velocity. Such higher velocity causes more significant erosion, and consequently the fracture aperture change is faster.

As the infiltration continues, the average fracture aperture increases by approximately 200%, from 0.1 mm to nearly 0.3 mm (Figure 16). At the initial stage, the fracture is subject to random distribution, with minor difference in the apertures. However, the variation diverges under the effects of the erosion. The fractures with smaller initial apertures become slowly enlarged by approximately 20%, while those with larger apertures continuously increase by more than 200%. As the fractures gradually dilate and interconnect with each other, a number of dominant seepage channels are formed, as illustrated in Figure 14(f).

4.2. Distribution of Water Seepage Velocity

The seepage velocities (i.e., Darcy velocity) in the model domain and along O-O′′ and O′-A′ lines were shown in Figures 17 and 18, respectively. The predicted failure time using the inverse velocity theory was plotted in Figure 19. The key observations include the following:

As shown in Figure 17, the changes in the flow velocity are consistent with the distribution of the fractures (Figure 14), and more fracture aperture changes are observed at the locations where higher seepage velocities take place. As shown in Figure 18, at 1.4 × 104 s, the minimum and maximum seepage velocities were approximately 3.5 × 10-6 m/s and 4.11 × 10-5 m/s. At 1.8 × 104 s, these two figures increased by 20% and 200%, to 4.2 × 10-6 m/s and 1.25 × 10-3 m/s, respectively.

There was a continuous increase in flow rate at the outlet, from approximately 3.7 m3/h initially, to approximately 100 m3/h at 1.8 × 104 s, increased by nearly 30 times. However, its increase showed a nonlinear trend, likely due to the fact that the increase rate was relatively slow in the initial stage but accelerated gradually until the water inrush occurred. Figure 19 was generated according to the inverse velocity theory, and the results showed that the water inrush would occur at 1.85 × 104 s for this particular flow condition. The trend has been reflected by the change in flow velocity shown in Figure 17. The figure shows that shortly before the inrush, the flow channels rapidly connected together and formed a main channel, resulting in a surge in water flow.

4.3. Distribution of the Particle Concentrations

Particle concentration is directly related to the erosion rate. Higher particle concentration indicates more severe erosion, which is more likely to result in greater fracture aperture change. Figure 20 shows the evolution of the average particle concentration along O-O′′ and O′-A′ lines. The average particle concentration and particle loss rate was plotted in Figure 21.

Figure 20 plots the spatial distribution of the particle concentrations for the whole model at different times. Results show that the particle concentration gradually increases from the bottom upwards. This is because, under the effects of the erosion, increasing amount of particles migrated into the water and gradually moved upwards through the fractures and pores. This is consistent with the change in the particle concentration at different locations, as shown in Figure 21.

Results in Figure 22 are well in line with the results shown in Figure 18. For example, at the beginning, the particle loss rate is relatively low; however, the rate constantly increases over time and culminates at the inrush. This is due to the fact that, along with the increases in the seepage velocities, the erosion effects on the solid particles are constantly intensified, which results in frequent migration and substantial particle loss.

4.4. Sensitivity Analyses

The outburst processes were affected or controlled by various factors, including the spatial heterogeneity of lithology of the KCPs, water pressure, absolute aperture, and distribution of the fractures, distribution law of concentration of the karst particles, and their sensitivity to water erosion. Based on the newly developed coupling model, this work further examined the water and particle migration characteristics in the KCPs under different hydrological and geological conditions, as well as the variations in the inrush occurrences. The scenarios to be studied are summarized in Table 2.

4.4.1. Impact of Heterogeneity on Fractures Evolution

Figure 23 illustrates the distribution law of the fractures, seepage velocities, and particle volume concentrations under identical flow pressures, but different fracture distribution conditions, at 17.5 × 103 seconds. The results show that more interconnected fractures are developed along with the increase in m values. For example, when , there are only two interconnected fractures. However, as m value increases, an increasing number of fractures were generated and interconnected. In the fourth case (), three interconnected fractures were generated. This finding indicates that the homogeneity can accelerate the propagation of the fractures. An increase in flow velocity and particle concentration with increasing m values was also observed, as shown in Figures 23(b) and 23(c).

4.4.2. Prediction of Water Inrush under Different Conditions

The analysis and prediction of the occurrences of water inrush in KCPs, based on the changes in seepage, are essential for the safety of mining processes. Some researchers have adopted catastrophe theories to predict the occurrences of water inrush. However, an approach that can be used to predict water inrush based on water flow has yet to be found in the literature. This study presents the first attempt to resolve this challenge by applying the inverse velocity theory, which has been widely approved effective in many fields such as slope stability prediction. The relevant results and analyses are shown in Figure 24–Figure 27.

As shown in Figure 24, water flow rate exhibits a nonlinear increase, and as the value of m increases, the seepage rate accelerates gradually. For example, at  s, the seepage rate was approximately 15 m3/h for , while it jumped to 25 m3/h for . Results from Figure 24(b) show that water inrush occurs sooner for a greater value. In other words, the inrush occurs earlier as the heterogeneity increases. The water inrush time follows the power function. This is due to the fact that the water is more prone to flowing through the existing fractures (which indicates a higher permeability) in the heterogeneous structure. Nevertheless, due to the heterogeneous distribution of the aperture, the fractures cannot be easily interconnected, which hampers the formation of the seepage channels. In regard to the homogeneous rock strata, since the flow of water was dominantly controlled by the pressure gradient, the fractures generated during the process of erosion could be easily interconnected, and as a result the water inrush occurred sooner.

Figure 25 illustrates the impact of the initial KCP’s water pressure on the flow rates and water inrush time. It can be seen that, as the pressure increases, the changes in the flow rates accelerate. For example, for  MPa, the flow rate was only approximately 5 m3/h, and it increased to 30 m3/h when  MPa. In addition, with increasing the aquifer water pressure, the water inrush time was also brought forward, as shown in Figure 25(a). The outburst time was approximately  s for p = 1 MPa, while only  s for  MPa. The relationship between the time of the inrush occurrences and the aquifer water pressure follows the power function, as illustrated in Figure 25.

For varying initial particle concentrations, the evolution of both the water flow rate change and the water inrush time shows a similar trend to that of different aquifer water pressures. The flow rate accelerates with increasing initial particle concentration, and consequently the water inrush occurs earlier, as shown in Figure 26. The outburst times also show a power function to the initial particle concentration.

The initial average fracture aperture also shows significant impact on water flow behavior. For instance, results from Figure 27 clearly illustrate that the flow rate increases with increasing initial average fracture aperture. Please note that the inrush time follows an exponential function to the initial average fracture aperture.

These numerical results can be explained by the fact that the initial seepage velocity is expected to be higher for the modelling cases with the higher initial aquifer pressure, initial fracture aperture, and/or initial particle concentration. The higher flow velocity would intensify the erosion of the fractures, enhance the permeability, and thereby accelerate the migration of the particles and interconnection of seepage channels, which causes water inrush.

5. Conclusions

In this study, a systematic approach was established to investigate the impact of controlling factors on water flow behavior for KCPs and to predict water inrush time by introducing the inverse velocity theory. In this approach, a suite of governing equations was developed to couple all the physics involved in this complex process, including water flow, rock erosion, and change in rock microstructure. These equations were then implemented into a finite element software COMSOL Multiphysics to simulate the fully coupled process. Based on this numerical model, a series of sensitivity studies were conducted, and some important findings have been drawn:(i)Water flow rate increases with increasing aquifer pressure, initial particle concentration, and fracture opening. This results in an increase in the permeability and finally an acceleration of the migration of the particles and interconnection of the seepage channels. This finding suggests the grouting of the KCPs could be an effective measure to mitigate water inrush risk.(ii)The inverse velocity theory was successfully applied to analyze the occurrences of water inrush under different conditions. The results show that the time of inrush has a power relationship to the rock heterogeneity, water pressure, and initial particle concentration and follows an exponential relationship to the initial fracture apertures.

The framework developed from this work not only helps enhance coal mining operations safety by better understanding and predicting water inrush risks, so that proper mitigation measures can be put in place timely. It can also be extended to other fields such as tunneling, backfill paste piping, tailing dam erosion, and other engineering applications, which are likely to encounter time-dependent erosion or deformation, or water inrush issues.

Please be noted that there are a number of other stochastic methods available in the literature that can be used to illustrate rock heterogeneity, such as the Monte-Carlo method or First-Order approach. The differences of these methods on the results of water inrush prediction may be worth comparing in the near future.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This project is supported by the National Natural Science Foundation of China (nos. 51304072, 51774110), the Program for Innovative Research Team in University of Ministry of Education of China (IRT_16R22), and the Postdoctoral Science Foundation of China (2017M612398).