Nonlinear seepage erosion model of water inrush considering particle size distribution of karst collapse column and its engineering applications

Jan 09, 2024

Scientific Reports volume 12, Article number: 17078 (2022) Cite this article

443 Accesses

1 Citations

Metrics details

Water inrush through karst collapse column is one of the great disasters which threaten coal mine safety production. The particle size distribution of karst collapse column is one of its most basic physical properties, which has a strong correlation with particle migration, and is an important basis for evaluating the water inrush risk of collapse column. The nonlinear flow tests of broken rock under different gradation conditions were carried out by a custom-built apparatus, and the relationship equation between nonlinear flow parameters (permeability and non-Darcy factor) and Talbol power exponent n were constructed. A nonlinear flow model with variable mass of water inrush from karst collapse column was established. The spatio-temporal evolution law of pressure, velocity, porosity and concentration under particle loss and the influence of particle gradation on the water inrush risk of karst collapse column at Fan gezhuang mine were discussed. During the water inrush, the flow state of fluids in karst collapse column gradually transitions from a weak inertial flow to a strong one, eventually becoming a turbulent flow. The flow model based on single flow state cannot reflect the essence of flow regime transition in water inrush. The larger n is, the stronger the water permeability of the karst collapse column, the faster the particles migrate and are lost, the faster the flow channel with high porosity develops, the shorter the time for the water inflow to reach its peak value, and the greater the risk of water inrush.

Karst collapse columns, a unique hidden vertical structure found in the coalfields of the Boreal Permo-Carboniferous system in China, are widely distributed in 45 coal mining areas of 20 coal fields1. The bottom of karst collapse columns is situated in a cave of soluble rocks. Generally, it may serve as a strong water-conducting channel that not only transmits hydraulic connections between an aquifer and a working face of an Ordovician limestone but also induces water inrush. For example, the ‘3·1’ major water inrush accident that took place in the Camel Hill Mine in Wuhai City, Inner Mongolia, on 1 March 2010, caused 32 deaths, 7 injuries and a direct economic loss of RMB 48.53 million2. On 10 September 2018, a water inrush accident from the seam floor of the 1313 coal face occurred in the Xiaoyun coal mine in Shandong Province. During this accident, the peak inrush water range reached 3,673 m3/h, resulting in a direct economic loss of up to RMB 25.66 million3.

Scholars at home and abroad have extensively investigated karst collapse column water inrush. Wang4 investigated the water inrush rules of karst collapse columns at the floor and the coal passing stratum by developing a laboratory table for water-conducting karst collapse column water inrush simulation. Zhang5 designed a three-dimensional large-sized simulation experimental model to reproduce the karst collapse column water inrush situation in Luotuoshan coal mine and acquire both the clinical hydraulic pressure and points of water inrush. Karst collapse column water inrush has rarely been studied using similarity simulation experiments, but the evolution law of non-linear seepage parameters has been frequently investigated in relation to fractured rock masses. For example, Moutsopulos6 discovered that both linear and non-linear resistance terms in the Forchheimer equation are inclined to decrease as large particle contents in porous media increase using non-linear seepage tests for accumulative granular porous media. Neild7 concluded that non-linear seepage parameters are sensitive to changes in pore structures based on research findings obtained by Dupuit and Forchheimer. The process of fractured rock mass seepage is usually combined with changes in the pore structure. Bai8, Yao9, Chen10, and Ma11,12 investigated variable mass flows of fractured rock masses in line with Darcy's law. The impacts of variable mass flows and particle size distribution of fractured rock masses on non-linear seepage parameters under non-linear flow conditions have rarely been reported.

In the numerical simulation, Huang13 not only used Darcy's law to describe the seepage behaviour of rock masses, but they also constructed a rock stress-rock damage-seepage coupling model based on the relationship between rock mass damages and permeability. On this basis, the influence of development height and water pressure of concealed collapse column on coal floor water inrush was studied. Yao14 constructed a mechanical model for deformations, seepage and erosional forces incurred by karst collapse column water inrush under multi-field coupling conditions by generally considering rock skeleton deformation, water seepage and filling particle migration during the seepage of karst collapse columns. Nevertheless, their model does not manifest non-linear seepage characteristics. Zhao15 built a model that integrates non-linear coupling seepage and conduit flows into one by combining a fluid–structure coupling theory, a flow state conversion theory and a strength-reduction approach of rock masses in order to analyse flow state conversions of bursting water in pressure-bearing caves. And highlighted the importance of whole-process analysis on flow state variations and water inrush to reveal the seepage-inducing water inrush mechanism. The author's team16,17,18 investigated uniform dynamic characteristics for three flow fields, namely water-bearing strata, water-conducting channels and roadway, from the perspective of flow state conversion. On this basis, Yang19 proposed a flow erosion model of water inrush in a fractured zone in the Jiangjiawan Mine, that couples the Darcy, Forchheimer, and Navier–Stokes fields under the theory of continuum mechanics, and the effects of rock disintegration and the coupled effects of flow and erosion were incorporated. However, the above model does not consider the influence of karst collapse column particle size distribution on water inrush.

In summary, there are two flaws in the existing literature on the disastrous mechanism of karst collapse column seepage. First, the sensitivity of non-linear seepage parameters to particle size distribution remains unclear; second, the flow state conversion mechanism caused by particle migration during the non-linear seepage of karst collapse column water inrush is still unknown. Therefore, a non-linear seepage testing system for porous media independently developed at the Northeastern University was used in this study to perform non-linear seepage tests on fractured rock masses under different particle size distribution conditions. An equation that expresses the relationship between non-linear seepage parameters and the Talbol power exponent was established in this manner. Thereafter, a non-linear seepage model was constructed for the variable mass of karst collapse column water inrush to investigate the influence of rapid and non-linear seepages and particle size distribution on the dangers of karst collapse column water inrush caused by particle loss. Finally, the disastrous mechanism of water inrush caused by karst collapse column seepage was revealed, providing a research foundation for early warning, prevention and control of water inrush in mines, as well as reasonable prediction of the water inrush range, explorations of the evolvement characteristics of seepage fields under complex hydrogeological conditions, etc.

The experimental apparatus was performed at Northeastern University, China, used a custom-built apparatus to model high-velocity seepage in a porous medium20. It mainly consists of four parts: the experimental unit, the water supply system, the data measurement equipment, and the recording equipment. Figure 1 illustrates the system connections and principles. Solid blue represents water, red parts represent samples, and blue dashed lines represent circuits.

Experimental schematic diagram.

Karst collapse columns contain internally fractured rock masses of various grain sizes, which form a high-dimensional parameter space. To investigate the influence of particle size distribution on permeability and the non-Darcy factors, the Talbol gradation theory was selected to describe the size distribution of fillings in karst collapse columns21:

where pi is the mass fraction of particles in a particle size range of group i (%), di is the particle diameter (mm), dm is the maximum particle diameter (mm) and n is the Talbol power exponent, which is dimensionless.

To reduce the influence of the size effect of particles on the experimental results, the maximum diameter of particulate matter should be less than 1/5–1/6 of D, which is the inner diameter of the sample-loading barrel used in the experiment. Since D = 60 mm, the maximum particle size of the selected sand grains must be less than 10 mm. In Eq. (1), the Talbol power exponent of the test samples with five different particle size distributions was set at 0.5, 0.75, 1.0, 1.25 and 1.5. Table 1 shows the mass fractions of particles of various size ranges.

The test tube (height of 0.2 m, diameter of 0.06 m and wall thickness of 0.007 m) was made of organic glass, which permitted the visualization of the test process. During filling, particles should be filled by layers and compacted uniformly by a mallet to ensure the uniformity of the sample. Before starting the experiments, the vacuum pump was used to remove the air from the inside of the sample, and then the negative pressure in the test tube is used to slowly passed the water through the tube from the bottom to the top. This method can make the sample saturation up to 95%, effectively avoid the influence of air on the experimental results.

Water was pumped from the bottom of the experimental column, which flowed upwards through the test tube, and then flowed into the mass measurement device. The quality method was adopted for measuring flow rate experimentally. Three experiments were performed under each pressure condition, and the average value of the experimental results was used for analysis. A total of 5 samples, about 180 experiments need to be conducted. The water pressure upstream and downstream of the sample was measured by a pressure collector. The real-time experimental data were recording by a data collection software in the computer. It should be noted that, the measurements were carried out after the flow stabilized for at least 1 min.

The corresponding relationship between the pressure gradient and flow rate of the samples with five different particle size distributions was obtained. Thereafter, data fitting was completed using the Forchheimer equation, as shown in Fig. 2. The result of the corresponding relationship fitting was satisfactory, and R2 = 0.99 in all cases. For the same sample, the pressure gradient tended to grow non-linearly as the flow rate gradually increased, and a larger slope of the curve indicates that the non-linearity of the pressure gradient and flow rate is more apparent. When the pressure gradient changed within a range of 0–1.2 MPa/m and the Talbol power exponent increased, the superior fovea of the corresponding curves became more prominent. This indicates that porous media with a relatively high Talbol power exponent are likely to produce non-linear seepages when subjected to the same pressure gradient.

The relationship between pressure gradient and flow velocity under different power exponent of Talbol.

The porosity of the five samples was 0.28. This phenomenon simply means that the pores of samples of the same volume have the same total volume. The pores of particles in the sample had diverse forms because of differences in particle size distributions, further resulting in uncertainties in seepage paths and resistance levels of fluids in the sample. As shown in Figs. 3 and 4, non-linear seepage parameters may change along with the Talbol power exponent under corresponding testing conditions. As shown in these figures, the corresponding permeability increased by two orders of magnitudes, from a magnitude of 10−12 to 10−10, when the Talbol power exponent increased from 0.5 to 1.5. In terms of the non-Darcy factors, they decreased by four orders of magnitudes, from 108 to 104. As the Talbol power exponent increased, the rate of increase in permeability and the non-Darcy factors increased and decreased, respectively. As shown in Table 1, as the Talbol power exponent increased, the large and small particle contents increased and decreased, respectively. When the total porosity of the samples was the same, porous media formed by small particle accumulation had a large number of pores and a low mean pore diameter. However, porous media formed by large particle accumulation had a few pores and a large average pore diameter. Based on the comparison, the latter had a low tortuosity, which shortened the path required by lateral diversions of fluids and reduced flow resistance.

Variation curve of \(\frac{{k(1 - \varepsilon )^{2} }}{{\varepsilon^{3} d^{2} }}\) with n.

Variation curve of \(\frac{{\beta \varepsilon^{3} d}}{1 - \varepsilon }\) with n.

Scholars have extensively investigated the values selected for permeability k and the non-Darcy factor β, as presented in Table 2. Both the non-Darcy factor and permeability are inherent attributes of a medium. When the particle size gradation of samples change, the compositions of particles in the sample may become diverse and uncertain. There are no equations in Table 2 that can express the relationship between permeability, non-Darcy factor and particle size distribution. Thus, the Talbol power exponent was introduced to slightly modify and improve them.

where k is the permeability (m2), β is the non-Darcy factor (m−1), ε is the porosity, which is dimensionless, \(\overline{d}\) is the mean particle size (m) and α1, α2, α3 and α4 are the fitting coefficients obtained by fitting the experimental data.

Equations (2) and (3) were used to find the corresponding relationship between the non-linear seepage parameter and the Talbol power exponent via data fitting (see Figs. 3 and 4). Here, the correlation coefficient (R2) was 0.99 in all cases, and the α1, α2, α3 and α4 values were 5 × 10−5, 2.8, 206.6 and − 8.3, respectively.

A karst collapse column is composed of three phases: an aqueous phase, a rock skeleton phase and a movable fine particle phase.

The flow rates of movable fine particles and the water phase are identical at any moment, and they move together, neglecting flow velocity losses and energy loss caused by particle collisions during water flows carrying sand.

The rock skeleton phase in a karst collapse column is rigid; the skeleton remains unchanged during the loss of movable fine particles, and both water and the water–sand mixed fluid are incompressible single-phase Newtonian fluids.

The porosity of the model remains effective; in other words, connected pores in a karst collapse column are completely filled with the water and movable fine particle phases, and pores that are not connected are viewed as the skeleton of this column.

The model-solving area is completely saturated.

The permeability of Ordovician limestone aquifers is generally low, and the pressure gradient has a linear relationship with the flow velocity, and the seepage itself conforms to Darcy's law, which can be expressed as the equation below:

where kD represents the permeability of the Ordovician limestone aquifer (m2), μw is the coefficient of dynamic viscosity of water (Pa·s), vD is the velocity (m/s), PD is the pressure (Pa), g is the gravitational acceleration (m/s2) and ρw is the water density (kg/m3), D is only the subscript, representing the Darcy field.

The equation of continuity for non-steady seepages is as follows:

where t is the time (s), and q is the source (sink) strength. The q values were set to be positive and negative for the sources and sinks, respectively (s−1). In the case of passive sink flows, there exists q = 0. Moreover, ∇ is the divergence operator.

The equation of continuity for mixed fluids

Seepage produced by karst collapse column water inrush has a very complicated mechanism. It is a mass-variable non-linear seepage22. Provided that water–sand mixed fluids are regarded as the incompressible single-phase Newtonian fluids, the equation of continuity can be written as follows19:

where ρm is the density of a mixture composed of migrated particles and water in the channel at moment t (kg/m3), and ρs is the density of solid particles (kg/m3), F is only the subscript, representing the Forchheimer field.

The equation of motion for mixed fluids

where μm is the dynamic viscosity of water–sand mixed fluids (Pa s).

A porosity evolution equation

where λ is the suffusion coefficient (m−1), and its dimension is a reciprocal value of length, which can be obtained through laboratory tests. In addition, vFc is the critical flow velocity of particle initiation (m/s), and εmax represents the maximum porosity.

The concentration transmission equation for migrated particles.

Assume that the mixed fluids are incompressible. Furthermore, the diffusion effects of particles in mixed fluids are ignored. The concentration transmission equation for migrated particles was obtained using the law of conservation of mass of fluids and the conventional constitutive equation of seepage suffusion23,24,25:

where cF is the concentration of migrated particles at moment t or the volume fraction of solid particles in mixed fluids.

The equation for changes in the density of mixed fluids

The relationship between the dynamic viscosities of mixed fluids and water26

Equations (2) and (3) were used to describe the relationship between permeability, non-Darcy factor and particle gradation of the karst collapse column fillings. Equations (2), (3), (8)–(11) are all auxiliary equations.

The advection–dispersion equation

The transport properties of particles in a roadway can be described using the advection–dispersion equation below27:

where Kdep and Kdis are the sedimentation and dispersion coefficients, respectively (s−1). Since mixed fluids move rapidly in a roadway, the corresponding diffusion effects are insignificant when compared to those of advection, and Kdis = 0. N is only the subscript, representing the Navier–Stokes field.

Navier–Stokes equations

The continuity equation for mixed fluids

Neighbouring fields must meet three conditions in order to combine three flow fields: pressure balance, velocity continuity and concentration continuity. The model can be solved if the Forchheimer equation set is closed. There are 8 uncertain model parameters in Eqs. (2), (3), (6)–(11): PF, vF, cF, ε, k, β, ρF and μm. The system of equation is sufficient to solve all the equation, and properly determine the values for the 8 uncertain model parameters. A weak form of fluid mechanics equation was constructed based on the virtual displacement principle using the FELAC2.2 software (Finite Element Language And it's Compiler). The convective term was dispersed using a finite volume method, while the remaining terms were dispersed using a finite element method that was performed to numerically solve the model. For details, please refer to the literatures Ref.19 and Ref.28.

The Fangezhuang Mine is located to the southeast of the Kaiping Coalfield in Tangshan City, Hebei Province. Its stratum is oriented NE–SW, inclining towards NW. Furthermore, the corresponding stratigraphic dip is generally within a range of 8°–24°. Folds and fractures are developed in this area. According to the structural characteristics, the minefield is divided into 3 structural areas, namely Tatuo syncline area in the north, central monocline structural area and Bigezhuang syncline area in the south29. On 2 June 1984, a catastrophic water inrush from the karst collapse column of the Ordovician limestone occurred at the 2171 working face of the Kailuan Fangezhuang Mine. This is rare in the history of mining. According to calculations based on the flooding volume, the average capacity of water inrush was up to 2053 m3/min in peak hours. The water inrush channels of the 2171 working face have been proven to be Ordovician limestone karst collapse columns that have an extremely strong water-conducting capability and are hidden within the coal face based on the characteristics of the water inrush procedure and capacity, as well as data collected from extensive drilling, geophysical prospecting and hydrogeological tests, etc., during water control30. Figure 5 depicts the position and spatial forms of the karst karst collapse columns. The total volume of the karst collapse columns is 861,000 m3. There is a large cave (volume: 39,000 m3) that is 8–32 m high above a location equivalent to the 7 s coal seam. However, the karst collapse column damages the integrity of a 280-m thick stratum that runs from the Ordovician limestone to the roof of the 5 s coal seam, forming a water-conducting channel well. Thus, high-pressure karst water of the Ordovician limestone can interface with the water-bearing stratum of sandstone at the 5 s coal seam roof. The karst collapse column has a long axis of about 67 m and a short axis of about 46 m, with an area of about 2875m2. The diameter of its spatial form decreases in a segment from 12 to 14 s coal seams. Most of the soft rocks used as fillings in the karst collapse column have aged and softened, as evidenced by accumulations generated by core drilling and those rushing in the roadway. The specific porosity values for the upper, middle and lower segments were calculated to be 0.21, 0.62 and 0.047, respectively.

The location and spatial distribution of the Karst collapse column.

Ordovician limestones are buried shallowly in the eastern and northern parts of the field but deep in its west and south. Outside the field, it is a subcrop in direct contact with the Quaternary unconsolidated aquifer. The Ordovician limestones have been determined to be an inter-connecting karst water-bearing body based on field pumping tests and long-term dynamic observations on the Ordovician limestone aquifer. However, the Ordovician limestones have dramatically non-uniform water yield properties. To be specific, some boreholes in the north of the well field have a unit water inrush capacity of 6.59 L/(s m) and a permeability of up to 31.87 m/d; however, some boreholes in the southern part have a unit water inrush capacity of less than 0.01 L/(s m). In normal cases, no direct water-filling relationship is formed between the aquifer and the mine. Owing to the existence of karst collapse column, the Ordovician limestone water is directly channelled into the coal measure strata, making it a direct water supply source for mines. Furthermore, the hydraulic pressure of the Ordovician limestone aquifer is 9 MPa29.

A numerical simulation model for karst collapse column water inrush was developed based on the forms of karst collapse column provided in relevant geological data. The model is composed of three parts. Darcy's equation was used to express the seepage of the lower Ordovician limestone aquifer, the Forchheimer equations were used to describe the mixed fluid flow of the middle karst collapse column and Navier–Stokes equations were used to describe the upper free flows of roadways. The initial concentrations of migrating particles inside the karst collapse column and the roadway were set at 0.01 and 0, respectively, to account for migrated particle sedimentation in the roadway inside the karst collapse column. Hydraulic pressure of 9 MPa, which is equivalent to that of the Ordovician limestone aquifer, was applied towards the lower part of the model. In addition, the roadway export directly communicated with the air, and the relative pressure was 0. The other exterior boundaries of the model were all confined. Figure 6 depicts the geometry size and boundary conditions of the model. The solution domain was segmented into 35,000 structured grids using triangular elements. Table 3 lists the relevant calculated parameters.

Numerical model of the of the Karst collapse column.

Figure 7 depicts the spatio-temporal evolution process of hydraulic pressure for karst collapse column water inrush in the Fangezhuang Mine. Hydraulic pressure distribution in the aquifer and karst collapse column may change noticeably at different stages of water inrush. Such a phenomenon mainly occurs in the lower part of a karst collapse column and its neighbouring aquifer. In addition, the pressure decreases to the greatest extent at the karst collapse column, but its decline rate becomes gradually small and even remains unchanged in positions far from the karst collapse column. Therefore, we selected hydraulic pressure variations observed in the Ordovician limestone aquifer as an early warning of water inrush.

Space–time evolution process of pressure (unit: MPa).

The curves in Fig. 8 represent variations in pressure distribution from the Ordovician limestone aquifer to the roadway. As shown in this figure, as the time of water inrush progresses, the pressure on a contact position between the aquifer and the karst collapse column changes significantly, decreasing by 2.19 MPa from 8.31 MPa at the beginning to 6.12 MPa. At other locations, the pressure variation range is relatively narrow. The pressure distribution curves are consistent at t = 1200 s and 1800s, indicating that the fluid pressure field distribution in the karst collapse column reached a steady state. Therefore, hydraulic pressure constantly changes in the entire flow region during water inrush.

Curve of water pressure from Ordovician limestone aquifer to tunnel (A1-A2-A4-A5).

Figure 9 presents the spatio-temporal evolution of the flow velocity, whereas Fig. 10 depicts the flow velocity distribution curves from the Ordovician limestone aquifer to the roadway. As can be observed, there is a change in the velocity of fluids due to water inrush flood from the aquifer to the roadway through the karst collapse column. According to the spatial analysis results, the flow velocity increases from 4.32 × 10−4 to 1.99 × 10−3 m/s in a segment from A1 to A2 of the aquifer (t = 1800s), thereby becoming 4.6 times greater. After the fluids flow into the karst collapse column, the decompression action of this column causes a hydraulic disturbance in the aquifer. Consequently, the seepage balance state of the original aquifer is destroyed, and the flow velocity abruptly increases in a short time by an order of magnitude from 1.99 × 10−3 to 2.15 × 10−2 m/s at the narrowest location of the karst collapse column (i.e. A3). This is an early warning of a sudden increase in the flow velocity of the karst collapse column during a transient process of water inrush. Once it enters the roadway, the flow velocity further increases to 6.25 × 10−2 m/s. The flow velocity is two orders of magnitude higher than that at A1 (4.32 × 10−4 m/s), which is the export boundary of the aquifer. In addition, the flow velocity is transformed from uniform changes to ‘step changes’. According to the temporal analysis results, the steady flow velocity is 3.90 × 10−3 m/s in the roadway at t = 1 s, but it increases to 6.25 × 10−2 m/s at t = 1,800 s, which is a 16-fold increase. Therefore, karst collapse column water inrush can be regarded as a dynamic process in which seepage changes gradually at first before changing abruptly, and ‘step changes’ increases in flow velocity are the most visual reflection of water inrush occurrence and development.

Space–time evolution process of velocity (unit: m/s).

Curve of velocity of flow from Ordovician limestone aquifer to tunnel (A1-A2-A4-A5).

Figure 11 depicts the spatio-temporal evolution of porosity within the karst collapse column. Figure 12 shows that the porosity changes along with time along the A2A4 surveying line. Spatial analysis revealed that the porosity does not change uniformly. Internally filled small particles are activated at the narrowest site of the karst collapse column and its export, resulting in the formation of fluidised particles. The porosity of the above two positions increases sharply as sand-carrying water flows into the roadway. Moreover, their interconnection results in the formation of a preferential migration passage, decreasing the resistance to sand-carrying flows. Thus, conditions beneficial for flow acceleration can be created, enabling high-velocity flows, which facilitate the migration of large particles. This is an interactive variable mass process. The preferred migration passage continues expanding towards the lower part of the karst collapse column until all fine particles are lost and the porosity reaches its peak value. The temporal analysis results revealed that the porosity variations are minor at an initial time of t = 1 s. The preferred seepage channel had a prototype form at t = 180 s. However, the preferred seepage channel was noticeably formed at t = 600 s. Since then, the channel has constantly been expanding circumferentially because of the action of heavy water flows inside it, resulting in the initiation of particles that had been stationary. At t = 1800 s, almost all porosity values reach their peaks in the entire karst collapse column. This can be mutually verified with the spatio-temporal evolution law of pressure, velocity and concentration. When the porosity of the karst collapse column increases from 0.16 to 0.45, the Karst collapse colum ‘activates’ and gradually becomes a water-conducting channel, allowing the Ordovician limestone aquifer to communicate hydraulically with the roadway. In this process, the migration of small particles causes a sharp increase in porosity, which is a direct cause of water inrush.

Space–time evolution process of porosity.

Porosity distribution of measuring line A2A4.

Figure 13 illustrates the spatio-temporal evolution of concentrations for migrating particles in a karst collapse column and roadway, whereas Fig. 14 depicts the time-varying concentration curves at typical sites. The concentration has a variation law entirely consistent with that of porosity. Particles at the narrowest site (A3) of the karst collapse column are first fluidised, followed by those at the export (A4) and the bottom (A2) of the karst collapse column successively. At the monitoring point (A2), the particle concentration increases rapidly at first to its peak value and then sharply decreases, eventually becoming 0. Only rock skeletons that cannot be fluidised remain there. At A3 and A4, the particle concentration increases to its peak and begins to fluctuate around this peak for the following reason. Instead of no particle migration, the particles flowing in have the same volume as the particles flowing out. The migrating particles flood the roadway space as water flows, and their concentration gradually decreases to 0 during sedimentation. Thus, many sediments and rock fragments can be found at the working face and the roadway space after water inrush.

Space–time evolution process of concentration.

Concentration history curves of different typical points.

As shown in Fig. 15, the inlet pressure of a karst collapse column and the flow rate at its export are both very sensitive to particle size distribution. When the hydraulic pressure of an aquifer remains unchanged, the Talbol power exponent increases from 0.5 to 1.5. In this context, the steady hydraulic pressure at the inlet of the karst collapse column decreases by 2.72 MPa from 6.09 to 3.37 MPa. This indicates that an increase in the Talbol power exponent may improve the pressure relief effect as long as the water recharge amount of the aquifer remains constant. At the export of the karst collapse column, the flow rate increases as the Talbol power exponent increases, and the discharge per unit width increases from 0.14 to 0.27 m2/s, which is a 93% increase. As shown in Fig. 16, a high Talbol power exponent indicates that the water-percolating capacity of the karst collapse column is high. In addition, particles may be lost more rapidly, resulting in a dramatic increase in the particle concentration in the mixed fluid. The formation of a steady seepage channel with a high porosity also takes a relatively short time, and the water inflow reaches its peak value quickly. The time is shortened by 20% from 720 to 420 s. In such situations, the risk of water inrush becomes more dramatic. Thus, variations in particle size distribution are believed to have a significant effect on variations in water inrush pressure and flow rates at the export in karst collapse columns that are transition areas for aquifer seepage and free flows of roadway water inrush. In practice, water inrush is caused by karst collapse column ‘activation’, which is incurred by particle migration and loss under actions of sufficient water recharge to the aquifer and constant high hydraulic pressure.

Water inflow and pressure at the entrance of the karst collapse column under different values of n: (a) pressure; (b) water inflow.

History of porosity and solid volume concentration at monitoring point A3 under different values of n: (a) porosity; (b) solid volume concentration.

A non-linear parameter denoted as E was introduced and expressed as follows to quantitatively describe the non-linear flow characteristics of fluids in a karst collapse column32:

Equation (15) expresses the ratio of a pressure drop caused by inertia force to another pressure drop caused by both inertia force and viscosity. Generally, E = 0.1 is selected as the critical point at which the linearity of a fluid flow becomes non-linearity33. Shi22 emphasised that E should be 0.5, which is the critical point at which weak inertial flows change into strong ones, while E = 0.9 can be used as the critical point at which strong inertial flows become turbulent flows. As shown in Fig. 17, the curve can be roughly divided into three segments based on the E values: (1) weak inertial flows when 0.1 ≤ E < 0.5, (2) strong inertial flows when 0.5 ≤ E < 0.9 and (3) turbulent flows when 0.9 ≤ E < 1.0. Provided that the hydraulic pressure of the aquifer remains constant, a karst collapse column gradually evolves into a highly porous and strongly permeable water-conducting channel during fine particle migration. The flow state of fluids in the column does not remain unchanged; instead, it gradually transforms from a weak inertial flow to a strong inertial flow, eventually becoming a turbulent flow. A seepage model based on a single flow state cannot reflect the nature of water inrush flow state transitions. The Forchheimer equations can effectively represent the intermediate state for water flows changing from a Darcy flow in an Ordovician limestone aquifer to a turbulent flow in the roadway. Moreover, these equations have the potential to quantitatively reveal the flow state transition mechanism caused by particle loss during the entire karst collapse column water inrush process. In the flow region of the karst collapse column, the non-linear parameter (E) is above 0.4, which is significantly greater than the critical value of linear flows (i.e. 0.1). Thus, the model built in this study highlights the concept of non-linear flows.

History of E at monitoring point A3 under different values of n.

The Talbol power exponent n is introduced to properly modify and improve the traditional computational formula, which is \(k = 5 \times 10^{ - 5} \frac{{\varepsilon^{3} }}{{(1 - \varepsilon )^{2} }}\overline{d}^{2} n^{2.8}\) and \(\beta = 206.6\frac{(1 - \varepsilon )}{{\varepsilon^{3} \overline{d} }}n^{ - 8.3}\).

When the hydraulic pressure of an aquifer remains unchanged, the Talbol power exponent increases from 0.5 to 1.5. In this context, the steady hydraulic pressure at the inlet of the karst collapse column decreases by 45% from 6.09 to 3.37 MPa. And the discharge per unit width increases from 0.14 to 0.27 m2/s, which is a 93% increase. The formation of a steady seepage channel with a high porosity also takes a relatively short time, and the water inflow reaches its peak value quickly. The time is shortened by about 42% from 720 to 420 s.

During the water inrush, the flow state of fluids in karst collapse column gradually transitions from a weak inertial flow to a strong one, eventually becoming a turbulent flow. The Forchheimer equations can reveal the intermediate state of water flows transitioning from Darcy flows in an Ordovician limestone aquifer to turbulent flows in a roadway. Additionally, the equations have the potential to quantitatively uncover the flow state transition mechanism caused by particle losses throughout the entire karst collapse column water inrush process.

The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.

Yin, S. X. & Wu, Q. Wang Shangxu Studies on characters and forming mechanism of karstic collapse columns at mine area of north China. Chin. J. Rock Mech. Eng. 23(1), 120–123 (2004).

Google Scholar

Yang, T. H. et al. A non-linear flow model based on flow translation and its application in the mechanism analysis of water inrush through collapse pillar. J. China Coal Soc. 42(2), 315–321 (2017).

Google Scholar

Wu, C. L., Teng, Z. J., Zhang, R. T. & Mu, H. P. Treatment of water gushing from the collapse column in Xiaoyun Coal Mine of Shandong province. Drill. Eng. 48(3), 161–169 (2021).

Google Scholar

Wang, J. C., Li, J. B. & Xu, G. M. Development and application of simulation test system for water inrush from the water-conducting collapse column. J. Min. Saf. Eng. 27(3), 315–321 (2010).

Google Scholar

Zhang, W. Z. Three-dimension large-scale simulation on experimental research on water inrush of collapse column. J. Taiyuan Univ. Technol. 46(6), 685–690 (2015).

Google Scholar

Moutsopoulos, K. N., Papaspyros, I. N. E. & Tsihrintzis, V. A. Experimental investigation of inertial flow processes in porous media. J. Hydrol. 374(3–4), 242–254 (2009).

Article ADS Google Scholar

Nield, D. A. & Bejan, A. Convection in Porous Media (Wiley-Blackwell, 2013).

Book Google Scholar

Bai, H. B. Seepage characteristics of top stratum of ordovician system and its application study as key aquifuge. Chin. J. Rock Mech. Eng. 30(06), 1 (2011).

CAS Google Scholar

Yao, B. H. Research on variable mass fluid-solid coupling dynamic theory of broken rockmass and application. PhD Thesis (China University of Mining and Technology, 2012).

Chen, Z. Q., Wang, L. Z., Kong, H. L. & Yao, B. H. Study on nonlinear seepage test of variable mass of broken rock. Saf. Coal Mines. 45(02), 15–17 (2014).

Google Scholar

Ma, D. et al. Impact of particle transfer on flow properties of crushed mudstones. Environ. Earth Sci. 75(07), 593–612 (2016).

Article Google Scholar

Ma, D., Miao, X. X., Jiang, G. H., Bai, H. B. & Chen, Z. Q. An experimental investigation of permeability measurement of water flow in crushed rocks. Transp. Porous Media 105(03), 571–595 (2014).

Article Google Scholar

Huang, Y. & Lin, Z. B. Numerical analysis on water inrush catastrophe of concealed collapse column of floor under mining. Saf. Coal Mines. 52(5), 193–200 (2021).

Google Scholar

Yao, B. H., Wang, L. C., Wei, J. P., Li, Z. H. & Liu, X. J. A deformation-seepage-erosion coupling model for water outburst of Karst collapse pillar and its application. J. China Coal Soc. 43(7), 2007–2013 (2018).

Google Scholar

Zhao, Y. L. et al. Solid-fluid coupling-strength reduction method for karst cave water inrush before roadway based on flow state conversion theory. Chin. J. Rock Mech. Eng. 33(9), 1852–1862 (2014).

Google Scholar

Yang, B. Non-Darcy flow model for water-inrush in mine and its initial application. Master thesis (Northeastern University, 2015).

Yang, T. H., Shi, W. H., Li, S. C., Yang, X. & Yang, B. State of the art and trends of water-inrush mechanism of nonlinear flow in fractured rock mass. J. China Coal Soc. 41(7), 1598–1609 (2016).

CAS Google Scholar

Shi, W. H., Yang, T. H., Liu, H. L. & Yang, B. Numerical modeling of non-darcy flow behavior of groundwater outburst through fault using the forchheimer equation. J. Hydrol. Eng. 23(2), 1–9 (2018).

Article CAS Google Scholar

Yang, B., Yang, T. H. & Hu, J. Numerical simulation of non-darcy flow caused by cross-fracture water inrush, considering particle loss. Mine Water Environ. 40(2), 466–478 (2021).

Article Google Scholar

Yang, B. et al. Impact of particle-size distribution on flow properties of a packed column. J. Hydrol. Eng. 24(3), 04018070.1-04018070.11 (2019).

Article Google Scholar

Talbot, A. N. & Richart, F. E. The strength of concrete, its relation to the cement aggregates and water. Illinois Univ. Eng. Exp. Sta. Bull. 137, 1–118 (1923).

Google Scholar

Shi, W. H. A non-Darcy flow model for water inrush through broken rock mass and its engineering application. PhD Thesis (Northeastern Univ, 2018).

Cividini, A. & Gioda, G. Finite-element approach to the erosion and transport of fine particles in granular soils. Int. J. Geomech. 4(3), 191–198 (2004).

Article Google Scholar

Papanastasiou, P. & Vardoulakis, I. Sand erosion with a porosity diffusion law. Comput. Geotech. 32(1), 47–58 (2005).

Article Google Scholar

Stavropoulou, M., Papanastasiou, P. & Vardoulakis, I. Coupled wellbore erosion and stability analysis. Int. J. Numer. Anal. Meth. Geomech. 22(09), 749–769 (2015).

3.0.CO;2-K" data-track-action="article reference" href="https://doi.org/10.1002%2F%28SICI%291096-9853%28199809%2922%3A9%3C749%3A%3AAID-NAG944%3E3.0.CO%3B2-K" aria-label="Article reference 25" data-doi="10.1002/(SICI)1096-9853(199809)22:93.0.CO;2-K">Article Google Scholar

Fei, X. J. Viscous coefficient (stiffness coefficient) of high concentration muddy water. Shuili Xuebao. 03, 59–65 (1982).

Google Scholar

Zhang, P. Y., Bai, B. & Jiang, S. C. Coupled effects of hydrodynamic forces and pore structure on suspended particle transport and deposition in a saturated porous medium. Rock Soil Mech. 37(5), 1307–1316 (2016).

CAS Google Scholar

Shi, W. H. et al. Non-Darcy flow model and numerical simulation for water-inrush in fractured rock mass. Chin. J. Rock Mech. Eng. 35(3), 446–455 (2016).

Google Scholar

Meng, Z. P., Yi, W., Lan, H. & Wang, M. Water inrush characteristics of fangezhuang coalmine field in kailuan and its geological condition analysis of water inrush from coal seam floor. Chin. J. Rock Mech. Eng. 28(2), 228–237 (2009).

Google Scholar

Zhu, W. C., Wei, C. H., Zhang, F. Z. & Yang, T. H. Investigation of water inrush from karst subsidence columnby usinga coupled hydromechanical model. Chin. J. Underground Space Eng. 37(5), 928–933 (2016).

Google Scholar

Yin, S. X., Lian, H. Q., Liu, D. M. & Yin, H. C. 70 years of investigation on Karst collapse column in North China Coalfield: Cause of origin, mechanism and prevention. Coal Sci. Technol. 47(11), 1–29 (2019).

Google Scholar

Zimmerman, R. W., Al-Yaarubi, A., Pain, C. C. & Grattoni, C. A. Non-linear regimes of fluid flow in rock fractures. Int. J. Rock Mech. Min. Sci. 41(supp-S1), 163–169 (2004).

Article Google Scholar

Javadi, M., Sharifzadeh, M., Shahriar, K. & Mitani, Y. Critical Reynolds number for nonlinear flow through rough-walled fractures: The role of shear processes. Water Resour. Res. 50(2), 1789–1804 (2014).

Article ADS Google Scholar

Download references

The current thesis was supported by the Educational Commission of Liaoning Province of China (LJKZ0322) and the Youth Foundation of University of Science and Technology Liaoning (2020QN10).

School of Civil Engineering, University of Science and Technology Liaoning, Anshan, 114051, China

Bin Yang

School of Civil Engineering, Suzhou University of Science and Technology, Suzhou, 215009, China

Wenhao Shi

School of Civil and Transportation Engineering, Guangdong University of Technology, Guangzhou, 510006, China

Xin Yang

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

B.Y. and X.Y. conceived and designed the experiments; B.Y. and W.S. established the numerical model; B.Y. wrote the main manuscript text.

Correspondence to Bin Yang.

The authors declare no competing interests.

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

Yang, B., Shi, W. & Yang, X. Nonlinear seepage erosion model of water inrush considering particle size distribution of karst collapse column and its engineering applications. Sci Rep 12, 17078 (2022). https://doi.org/10.1038/s41598-022-21623-4

Download citation

Received: 12 June 2022

Accepted: 29 September 2022

Published: 12 October 2022

DOI: https://doi.org/10.1038/s41598-022-21623-4

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.