Predicting Erosion-Induced Water Inrush of Karst Collapse Pillars Using Inverse Velocity Theory

1School of Safety Science and Engineering, Henan Polytechnic University, Jiaozuo, Henan 454000, China 2State Key Laboratory of Mining Disaster Prevention and Control Co-Founded by Shandong Province and the Ministry of Science and Technology, Shandong University of Science and Technology, Shandong 266590, China 3Collaborative Innovation Center of Coal Work Safety, Henan Province, Jiaozuo 454000, China 4School of Mechanical and Mining Engineering, The University of Queensland, Brisbane, QLD 4072, Australia


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 m 3 /min, resulting in the mine and other three neighboring mines submerged in a very short time, and caused a direct capital loss of around US$90 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

Working face
Figure 1: Illustration of a Karst Collapse Pillar: front view of a typical cone shaped KCP (not to scale) [23] and different shapes of rock layers represent different types of rocks.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 mininginduced damage on KCPs and the surrounding rocks and on the formation of the fracture zone and analyzed mininginduced KCP groundwater inrush risk.This work was later extended to study the effect of coal mining operations on KCPs related groundwater inrush using FLAC 3D , 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 [13][14][15].
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 [20][21][22], 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.

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 [24][25][26], 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.

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.

Definition.
Based on the aforementioned assumptions, a representative element for a typical KCP microstructure is shown in Figure 5.The length  (m) and volume  (m 3 ) of the element are expressed as where  (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  ≪ .

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   ( = 1, 2, 3), 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/m 3 ) denotes the density of solid particles, and  (m 2 /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 as where ṁ (kg/m 3 /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

Water Mass Conservation Equations.
For the water flow, the mass flowing out of the element in the direction of   ( = 1, 2, 3) can be expressed as where   (kg/m 3 ) 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  1 and  2 are constant, is the absolute value of the seepage velocity (m/s), and  max and  max 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.

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 (m 2 ).For this study, the flow velocity is relatively low; thereby, the acceleration term can be neglected, giving where  →  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 (m 2 ), which can be expressed as [32] in which 2 and   =  3 /12 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 (1) the sizes of the infilling particles are significantly smaller than the fracture apertures (m versus mm); (2) 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 (3) 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].

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  < 2, and it is linear when  = 2, while it is convex when  > 2. 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. 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.

Numerical Model Implementation
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 m 3 /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 m 3 /h to 150 m 3 /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.

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.

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  0 was 0.1 mm; and initial porosity  0 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 [40][41][42].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,  0 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).

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: (1) 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.
(2) 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).18, respectively.The predicted failure time using the inverse velocity theory was plotted in Figure 19.The key observations include the following:

Distribution of Water Seepage Velocity
(1) 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 (2) There was a continuous increase in flow rate at the outlet, from approximately 3.7 m 3 /h initially, to approximately 100 m 3 /h at 1.8 × 10 4 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 × 10 4 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.

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.

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.

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  = 1.5 × 10 4 s, the seepage rate was approximately 15 m 3 /h for  = 4, while it jumped to 25 m 3 /h for  = 10.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 occurredsooner.
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  = 2 MPa, the flow rate was only approximately 5 m 3 /h, and it increased to 30 m 3 /h when   = 3 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  = 4.0 × 10 4 s for p = 1 MPa, while only 1.0 × 10 4 s for  = 4 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.

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.

Figure 2 :
Figure 2: The distribution of major coal mines in China and the coal fields in Shanxi province.

Figure 5 :
Figure 5: Sketch of a representative element of a typical KCP.

Figure 6 :
Figure 6: Illustration of a three-dimensional characteristic element in a KCP.  and   are mass change of solids due to convection and diffusion, respectively, defined as   =     ,   = −((  )/  ).

Failure a > 2 Figure 8 :
Figure 8: Inverse velocity trend lines to predict the time of failure [27].

Figure 9 :
Figure 9: Verification of the inverse velocity theory in predicting water inrush.

Figure 11 :Figure 12 :
Figure 11: (a) Mineral components of infilling materials for Sample 1; and (b) mineral components of infilling materials for Sample 2.

Figure 13 :
Figure 13: Model geometry and setup: (a) illustration of KCP and rock strata; (b) 2D simulation model; and (c) distribution of aperture opening based on Weibull function.
. 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

Figure 15 :Figure 16 :
Figure 15: Evolution of fracture aperture along the lines of (a) O-O  axis and (b) O  -A  axis.

Figure 23 :
Figure 23: Distribution of fracture, seepage velocity, and particle concentration at different Weibull distribution conditions.(a) Impact of Weibull distribution on fractures evolution; (b) impact of Weibull distribution on water velocity evolution; (c) impact of Weibull distribution on particle concentration evolution.

Figure 24 :Figure 25 :
Figure 24: Impact of heterogeneity on water flow rate (a) and inrush time (b).

Figure 26 :
Figure 26: Impact of initial particle concentration on water flow rate (a) and inrush time (b).

Figure 27 :
Figure 27: Impact of absolute aperture on water flow rate (a) and inrush time (b).