#### Abstract

A large number of statistics indicate that broken rock mass always transforms into a flowing channel and leads to water inrush disasters in mining engineering, such as fault, karst, and strongly weathered rock mass. During the process of water inrush, the structure of the broken rock mass is constantly changing due to seepage erosion under high-velocity flow. Therefore, it is of vital importance to quantitatively evaluate the flow behavior of the water inrush related to the seepage erosion in order to prevent or reduce the risks. This study described a coupled nonlinear flow model, which couples the high-velocity seepage, the small particle migration, and the evolution of the broken rock mass structure. The model was verified firstly for simulation of nonlinear flow behavior by comparing with the traditional one. Then, the proposed model was used to simulate the evolution of particle migration and seepage properties of the water inrush through broken rock mass by a numerical case. The simulation results generally agree well with the existing experimental results. The simulations indicate that small particle migration causes the unstable characteristics of the seepage and the heterogeneity properties of the broken rock mass, which lead to the nonlinear flow behavior of the water inrush in both time and space. From a different perspective, it also indicates that the proposed model is capable of simulating the interaction of high-velocity seepage, small particle migration, and evolution of broken rock mass structure in the process of water inrush.

#### 1. Introduction

Water inrush is considered to be a kind of serious geologic disaster in tunnel construction and underground mining. A large number of water inrushes have taken place throughout history in China, leading to serious casualties and some irreversible damage on the hydrogeological environment in many cases [1–7]. Although these kinds of hazards have been known for decades, the mechanism for water inrushes remains elusive [7]. This may be due to the complexity of the water inrushes, such as multiphysics coupling characteristics, nonlinear properties, and time-variant characteristics of parameters [1, 8]. It is believed that, with the depth of mining increasing, the water inrush will be a greater potential threat to the deep mining.

According to the literature, water inrushes always occur through broken rock mass which has the properties of weak water stability, low strength, and cataclastic structure, such as fault [2], karst [9, 10], fractured zone [11], and strongly weathered rock mass [12, 13]. The broken rock mass is prone to erosion during seepage because of uncemented or poor cemented properties of small particles inside the broken rock mass. Therefore, the small particle loss caused by seepage erosion is considered a key factor in the water inrush through broken rock mass, which has been confirmed by experiments [8, 14]. Therefore, it is of vital importance to model the mass transfer and flow properties of the water inrush through broken rock mass, not only for revealing the mechanism of water inrushes induced by seepage erosion but also for the real-time prediction of water inflow.

It is no doubt that the flow properties play an extremely important role in the evolution of water inrushes through broken rock mass. Many simplified methods have been employed to address the water inrush but avoid the problems of groundwater flow during underground mining, such as the water inrush coefficient method [15], vulnerability index approach [11], and key strata model [16]. These research studies are effective in solving corresponding practice problems of water inrushes, but they cannot accurately predict the water inrushes because of oversimplification of the geological conditions. Other models have been built with the flow properties taken into account, such as the flow-stress coupling model [17] and damage-based hydromechanical model [7, 18]. These models can be used to simulate the mining-induced groundwater inrushes related to geologic structures but cannot be applied for the simulation of water inrushes associated with seepage erosion. Furthermore, when the water inrush occurred, the flow velocity is generally high and the flow behavior shows obvious nonlinear characteristics [19]. Accordingly, previous models based on the linear Darcy’s law may be inconsistent with the actual and lead to large error.

In the practical engineering of water inrush related to seepage erosion, the small particles will be fluidized and carried out of the broken rock mass gradually. This will inevitably affect the flow properties of the broken rock mass such as the porosity and permeability and then lead to the change in seepage fields [20]. Meanwhile, the changed seepage field will react on the seepage erosion, which promotes or restricts the evolution of the broken rock mass structure. Therefore, there is a coupling effect between the seepage erosion and the permeability properties of the broken rock mass. Liu et al. [8] has conducted an experimental research on mass transfer and flow properties of water inrush in completely weathered granite. Their research studies present encouraging results in the parametric interaction during the water inrush related to seepage erosion. However, due to the limitation of experimental conditions, only a few parameters over time can be tested; it is difficult to accurately measure the spatial distribution of parameters inside the sample. In this respect, numerical simulation can be adopted to further studies.

Numerical modeling is currently one of the commonly used methods in solving complex problems in rock seepage mechanics and engineering, such as coal bed methane recovery [21, 22], oil exploitation [23], groundwater resources assessment [24], and prevention of mine water hazard [1, 2, 7]. In this paper, a quantitative coupled nonlinear flow model was proposed to describe the flow behavior of water inrush through broken rock mass related to small particle migration. The model was used to investigate the flow mechanism of water inrush intending to gain an insight into the coupling mechanism between nonlinear flow and small particle migration.

#### 2. Theoretical Model Description

The accuracy and versatility of the theoretical model depend on factors considered and the refined description of the physical processes [22]. The flow process of water inrush associated with small particle migration is extremely complex because of coupling effect and phase transformation of small particles [8, 14]. For the process of water inrush, the coupled effect of the small particle migration and the nonlinear flow behavior may be important to be formulated. The formulation of this model may help deeply understand the evolution mechanism of the water inrush. It is well known that any theoretical model is established based on some assumptions and has applicability in certain cases. To build this nonlinear flow model accurately, four assumptions for the broken rock mass and the fluid must be presented firstly: (i) all pores in the broken rock mass are connected, and no sealed pore is taken into account; (ii) the broken rock mass is nondeformable during the small particle migration; (iii) the fluid mixed with fluidized particles is taken as single-phase Newtonian flow; and (iv) the broken rock mass is always saturated during the seepage. Hereby, the evolution of the flow, the small particle migration, and the porous medium are presented in mathematical language in the following sections.

##### 2.1. Continuity Equation for the Flow

From the assumptions above, the single-phase flow consists of water and fluidized particles. Therefore, the content of the small particles in the flow can be quantitatively described by volume concentration. Then, the density of the fluid with fluidized particles can be expressed aswhere is the fluid density; and are the density of water and small particles, respectively; and is the volume concentration of the fluidized particles in the flow. In general, the density of small particles is greater than that of water, so the fluid density will increase with the increase in the volume concentration. If no particles were fluidized, i.e., , the fluid density will reduce to the water density.

According to the fundamental theory of seepage flow in porous media [25], the continuity equation of the flow through the broken rock mass can be described aswhere is the real-time porosity of the broken rock mass; is the average flow velocity; is the time; and is the mass source in kg/(m^{3}·s), which can be interpreted as the mass newly added to the flow associated with small particle migration per unit time and volume; or put another way, it is the rate of small particle migration per unit volume, which can be presented in another form by some physical and mathematical approaches.

Consider a control volume, with volume and initial porosity , taken from the broken rock mass, because of small particle migration, the porosity of the control volume during a short time interval, , will change to . The pore volume increment of the control volume is . Therefore, the mass change of the control volume can be expressed as . By the principle of mass conservation, this must be equal to the increment of the fluidized particles in the flow within the control volume during given by . Hence, it can be obtained as

The mass source may now be expressed by

Accordingly, the continuity equation for the flow can be expressed as

##### 2.2. The Equation of Motion for the High-Velocity Flow

Groundwater flow in porous media is commonly described by the well-known Darcy’s law, which is limited to laminar flow through porous media at low Reynolds numbers (say ). As increases beyond this range, the pressure gradient versus flow velocity displays nonlinearity and the deviation from Darcy’s law will be more and more clear with the increase in flow velocity. Reviews of the subject are presented early by Scheidegger [26]. Many research studies indicated that the flow behavior of the water inrush through broken rock mass shows pronounced nonlinearity and can be expressed by the Forchheimer equation [1, 8, 19, 27–29]. Considering the effects of fluid viscous resistance and inertial resistance, the equation of motion for the high-velocity flow of water inrush can be expressed aswhere is the fluid pressure; is the fluid dynamic viscosity; is the permeability of the porous medium; and is called the non-Darcy coefficient or inertial coefficient with dimension of , which is a fundamental property of porous media. If the nonlinear effect is negligible, i.e., , the Forchheimer equation will reduce to the linear Darcy equation. In this equation, three parameters, i.e., the fluid dynamic viscosity, the permeability, and the non-Darcy coefficient, will vary with the small particle migration.

The permeability and the non-Darcy coefficient are two properties of the porous medium and also the key parameters of the Forchheimer equation. According to the classical Kozeny–Carman equation [30], the permeability can be expressed by porosity, as shown inwhere is the coefficient that depends on the properties of the porous medium, such as particle size, shape, and arrangement, which can be measured by tests.

Previous studies have reported that the non-Darcy coefficient has a negative correlation with the permeability. Many empirical equations were proposed to obtain the non-Darcy coefficient according to the permeability, which have been summarized by Wang et al. [31], Li et al. [32], and Yao et al. [33]. In this paper, a simple formula just associated with the root of the permeability was used to calculate the non-Darcy coefficient [1, 34–36], as presented inwhere is the form drag coefficient, also known as the Forchheimer coefficient. This equation is dimensionally homogeneous and contains only one coefficient, which is easy and convenient to measure.

According to the theoretical formula for the relative viscosity of muddy water with uniform small spherical particles derived by Einstein [37], the viscosity of the inrush flow with small particles can be roughly expressed aswhere is the dynamic viscosity of the fluid without fluidized particles.

##### 2.3. Transmission Equation for the Fluidized Particles

The transmission equation for the fluidized particles can be expressed by the convective-diffusive equation, which has been used for the descriptions of underground erosion and contaminant transmission [38–40]. However, for high-velocity flow through porous media, the effect of convection on the particle migration is much more significant than the effect of diffusion. In this case, the effect of diffusion can be neglected completely. And hence, the transmission equation for the fluidized particles will take the following form:

The right-hand side of Equation (10), , denotes the rate of newly generated fluidized particles expressed by the volume concentration. By introducing Equation (4) into Equation (10), the transmission equation for the fluidized particles can be written as

##### 2.4. Evolution Equation for the Porosity

Before the water inrush, small particles are part of the porous medium, but after that, they are part of the flow of the water inrush due to fluidization. Therefore, the variation of the porosity is directly caused by the fluidization. Based on the filtration tests on porous media with nonviscous particles, Sakthivadivel [41] proposed an equation to describe the evolution of porosity during water flow-induced particle migration, which can be expressed aswhere is the coefficient. It denotes the ability of particle fluidization by water and can be measured by tests. This equation has been used to describe the flow behavior of piping erosion in dam foundation and foundation pit engineering.

Prior research has shown that, in general, only some small particles and not all particles in the porous media can be fluidized. Therefore, the porosity of the broken rock mass has an upper limitless than 1. For this reason, the evolution equation for the porosity based on Equation (12) can be represented as

The equations mentioned above indicate that the pore structure of the broken rock mass, the properties of the fluid medium, and the seepage field are interacting and coupled in the flow model. It is difficult to solve the model analytically because of the coupling effect. However, when the boundary and initial conditions are specified, the model can be numerically calculated by the finite element method based on the Finite Element Language And it’s Compiler (FELAC), which has been used in the simulations of various engineering problems [42–44]. Figure 1 illustrates the flowchart of the mathematical model.

#### 3. Validation of the Model for Simulating Nonlinear Flow Behavior

It is well known that, if a water inrush disaster occurred in the underground mining engineering, there must be three essential conditions at the same time: the source, the channel, and the space for the water inrush. In fact, the evolution of the water inrush channel which connects the aquifer and the working face is the key factor for the process of water inrush. Therefore, the development of the flow-induced water inrush channel will be simulated based on the nonlinear flow model above.

As shown in Figure 2, a two-dimensional numerical model is established for the complex flow simulations through the broken rock mass. The simulation domain is 180 mm high and 60 mm in width. The fluid flows into the broken rock mass from the bottom and flows out from the top. The inlet pressure is denoted by , while denotes the outlet pressure, obviously, . Both sides of the domain are set as the impermeable boundary, which can be mathematically expressed as . The initial and the maximum porosity of the broken rock mass are denoted by and , respectively. The initial concentration of the flow in the saturated broken rock mass is denoted by , which should be set to a small value but not equal to zero according to Equation (13). The physical and seepage parameters used in the model are listed in Table 1. In addition, for convenience of discussion, five measuring points are marked at equal interval spacing on the vertical central axis of the numerical model. The five points along the flow direction are , , , , and , respectively, which are illustrated in Figure 2.

In order to validate the model for simulating the process of high-velocity water inrush, another numerical model with the equation of motion described by the typical Darcy’s law is also established which fails to take into account the nonlinear properties of the water inrush. Two cases are devised in the simulations to more fully comprehend the differences between the two models. For case I, it is considered a low pressure gradient, in which the inlet pressure is set as 450 Pa, while the outlet pressure is 0. It was simulated with a time interval of three minutes per step. As for case II, a relatively high pressure gradient is considered, so the hydraulic pressure of 0.18 MPa is applied on the inlet, and the same pressure as case I is applied on the outlet. This simulation is discretized into 120 steps with a time interval of five seconds per step.

Figure 3 shows the distribution of pressure, average flow velocity, porosity, and concentration on the central measuring point, , calculated by the two models. It is well known that the flow velocity depends on the pressure gradient. For case I, the flow velocity is so small that the inertial resistance can be neglect completely compared with the viscous resistance. Therefore, the results of pressure, average flow velocity, porosity, and concentration simulated by the nonlinear flow model are nearly the same as those computed by the linear Darcy model. However, as for case II, the pressure gradient is 1.0 MPa/m and the flow velocity can reach about 0.03 m/s obtained by the coupled nonlinear flow model. In this case, the inertial resistance or the nonlinear flow property has to be considered because of the high-velocity flow, which is currently recognized [28, 45]. Otherwise, the flow velocity will be ridiculously high, with a maximum velocity of 0.162 m/s calculated by the linear Darcy model, which is completely unrealistic and also unreasonable. Therefore, it is reasonable to conclude that the established model is effective for simulating the nonlinear flow behavior of the water inrush through broken rock mass. According to the simulations, the evolution of nonlinear flow behavior during the water inrush will be analyzed in detail in the following sections.

**(a)**

**(b)**

**(c)**

**(d)**

#### 4. Evolution of Nonlinear Flow Behavior of Water Inrush

##### 4.1. Evolution of Concentration of Fluidized Particles

The concentration of fluidized particles reflects the effect of the mass transfer and seepage processes. Figure 4 shows the spatial-temporal distribution of the concentration of fluidized particles. From time perspective, as shown in Figure 4(a), the concentration of fluidized particles inside the broken rock mass first increases, then decreases, and gradually tends to be stable over time. The increase in concentration is related to the large number of small particles and the low permeability. This is because the numerous small particles are fluidized but cannot be transferred in time. However, the concentration attenuation is because of both the increase in seepage ability and the reduction in the content of the small particles in the broken rock mass. From a spatial perspective, as shown in Figure 4(b), for any noninitial moment before stabilization, the concentration of fluidized particles increases gradually from the inlet to the outlet, and the closer to the outlet, the greater the change of the concentration. For instance, on the outlet, m, the maximum concentration can reach 0.56, but for the middle position, m, the maximum concentration can only reach 0.26. This is the reasonable results of the fluidized small particles flowing continually to the outlet. After all of the transportable small particles are fluidized and migrated, it will change into the steady seepage.

**(a)**

**(b)**

**(c)**

Similar tests have been carried out by Liu et al. [8, 12]. They studied the mass transfer and flow properties of water inrush in completely weathered granite by a self-designed large-scale triaxial testing system. The sample cylinder has a diameter and height of 100 mm and 300 mm, respectively. The sample was made by mixing completely weathered granite particles which were classified into five groups (0-0.25, 0.25-2.0, 2.0-3.0, 3.0-5.0, and 5.0-10.0 mm) according to the particle size. The sample was first saturated and confined by the axial and confining loads. The test started when a water pressure was gradually applied to the sample. The pressure and flow velocity were monitored, and the eroded particles and water flow were collected every 30 s. As reported by Liu et al. [8, 12], the fluidized particle concentration can be divided into two stages, namely, increased at first and decreased at last, which is shown in Figure 4(c). Just from the evolution trend of the concentration, the simulation results are in good agreement with the experimental results. Undeniably, the value of the concentration calculated is different with the experimental results. The difference may be caused by two kinds of reasons: one is that the parameters employed for the simulations vary from the experimental materials and conditions, such as particle properties, porosity, and hydraulic gradient; the other is that only one experimental concentration can be obtained by collecting the eroded particles, while the calculated concentration is dispersed in space and simulated at each position of the model.

The fluid is mixed with fluidized particles. Therefore, the fluid properties such as fluid dynamic viscosity and fluid density highly depend on the concentration of fluidized particles. Figure 5 gives the spatial-temporal distribution of fluid dynamic viscosity. Because of the linear relationship between the concentration and the fluid dynamic, the change tendencies between them are exactly the same.

**(a)**

**(b)**

##### 4.2. Evolution of Porosity of Broken Rock Mass

The spatial-temporal distribution of porosity is shown in Figure 6. As can be seen from Figure 6(a), the porosity of the measuring points all increases constantly and then tends to be stable over time. It is quite clear that the evolution of the porosity is associated with small particle migration. When the small particles are fluidized and carried out of the broken rock mass by water, the pore volume of the broken rock mass will increase, thus leading to the increasing porosity. If all transportable particles were fluidized and discharged by seepage erosion, the porosity will reach its maximum value. Therefore, it is revealed that the small particle migration is the key factor that induces the increase in porosity as well as the evolution of seepage parameters of the broken rock mass. As shown in Figure 6(b), although the broken rock mass is set to the homogeneous initial and final porosity, the porosity close to the outlet changes much more rapidly than that near the inlet. For instance, on the outlet, m, the porosity changes from 0.16 to 0.44 for 66 s, but for the middle position, m, 186 s is needed for the same increment of porosity. The results indicate that the small particles in the broken rock mass near the channel outlet of flow are easy to be fluidized. Furthermore, the increased porosity near the outlet will provide a good connected channel for the particle migration.

**(a)**

**(b)**

**(c)**

As reported by Liu et al. [8, 12] and Ma et al. [14], the porosity increases with small particle migration. During these tests, the porosity experiences a significant increase to milder increase and finally remains constant, despite the time needed to reach the steady value is quite different. Therefore, the evolution trend of porosity calculated is consistent with the experiments, as shown in Figure 6(c). Generally, the sample with high original porosity will lead to a big rate of porosity change. It is because the sample with lower porosity has a higher compactness and a smaller flow channel [12]. Furthermore, if the sample contains a high percentage of small particles, it would result in a higher rate of porosity change because of higher potential for particle discharge in the initial seepage [14].

##### 4.3. Evolution of Nonlinear Seepage Parameters

The spatial-temporal distribution of permeability and non-Darcy coefficient are presented in Figures 7 and 8, respectively. Since there is a power function relationship between permeability and porosity, the permeability as well as the porosity increases continually and then tends to be stable. As shown in Figure 7(a), when the porosity increases from 0.16 to 0.44, the corresponding permeability will change from to , showing an increase of about two orders of magnitude. However, the non-Darcy coefficient is inversed to the square root of permeability, so the spatial-temporal distribution of the non-Darcy coefficient is completely opposite to that of permeability, as shown in Figure 8. When the porosity increases from 0.16 to 0.44, the corresponding non-Darcy coefficient will decrease by an order of magnitude. It is known that the seepage behavioral way depends on the characteristics of the porous media under the given external condition. The permeability and the non-Darcy coefficient are two key parameters for the non-Darcy flow through porous media. The high permeability and the low non-Darcy coefficient mean the low viscous resistance and the low inertial resistance, respectively. Therefore, the evolution of the two parameters shown in Figures 7 and 8 is a sign of growing risk of water inrushes.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

The experimental results of the permeability by Liu et al. [8, 12] and Ma et al. [14] were illustrated in Figure 7(c). It shows that, no matter if the seepage is Darcy or non-Darcy, the permeability obtained by experiments is closely related to the porosity. The increasing permeability means that major channels for flow and particle migration are forming inside the sample. When the flow reaches steady-state condition, no more particles migrate out of the samples and the permeability remains constant; this is also in agreement with the calculated results.

##### 4.4. Evolution of Pressure and Flow Velocity

Figure 9 presents the spatial-temporal distribution of pressure. Although the pressures on the inlet and the outlet are fixed, the distribution inside the broken rock mass is not stable. The pressures distributed inside the broken rock mass all first decrease and then gradually go back to their initial state. The distribution of pressure is closely related to the pore structure of the broken rock mass. At the initial and final stable states, the pore is homogeneously distributed inside the broken rock mass. In other words, the broken rock mass is a homogeneous and isotropic porous medium. Therefore, the pressure is linearly distributed along the flow. For any unsteady state, the broken rock mass is a heterogeneity porous medium, which brings about the nonlinear distribution of pressure in the broken rock mass. Furthermore, the nonlinearity of the pressure distribution will be much more obvious with the increasing heterogeneity of porous media. As seen from Figure 9(b), the nonlinearity of the partial distribution of pressure at about s is the most strong. That is, the porous medium is highly heterogeneous around this time. This can be further confirmed by the porosity distribution from Figure 6(b).

**(a)**

**(b)**

The flow velocity is shown in Figure 10. It should be noted that the flow velocity is the average flow rate through the cross section of the broken rock mass or called specific discharge. Because the cross section of the broken rock mass is constant, the flow velocities at different measuring points are the same. However, the porosity of the porous medium is increasing continuously with small particle migration, so the permeability is enhanced. Although the pressure gradient at both ends of the broken rock mass is constant, the flow velocity still increases. For this computational condition, as the porosity increases from 0.16 to 0.44, the flow velocity shows up to about a tenfold increase, from 0.0029 m/s to 0.027 m/s.

**(a)**

**(b)**

Liu et al. [8] divided the increasing water inflow into three stages, the linear flow stage with a small amount of particle migration, the rapid flow stage with a large amount of particle migration, and the stable flow stage with few (or no) particle migration, as shown in Figure 10(b). It can be seen that the latter two stages can be simulated well by the coupled nonlinear flow model. However, the linear flow stage cannot be calculated, because the flow in the model is assumed nonlinear which is more decisive, long-lasting, and nonnegligible for the water inrush compared with the linear flow. This is also the main reason of the rapid change of the flow behavior, such as the concentration, the porosity, the permeability, and the velocity.

#### 5. Discussion

Water inrush through broken mass is a very complex process of disaster, which may relate to both mechanical erosion of particles and chemical erosion of rock minerals. This may be very common in a water inrush through fault which is considered a typical hazardous medium of water inrush [14]. In this kind of disaster, the water usually becomes muddy with the developing fault channel and slowly clears after the channel is formed [2]. This phenomenon resulted from small particles that are fluidized and carried out of the fault. The fluidized particles flow with water, but it is not soluble in water. Therefore, the number of fluidized particles can be described by the volume concentration other than the molar concentration. In other words, the proposed model just describes the water inrush caused by the mechanical erosion. In this situation, if the pore throat is not large enough or the flow velocity is not high enough, the fluidized particles will deposit and plug the pore, which usually does not lead to the water inrush disaster.

Development of the pore structure is the key factor for the evolution of water inrush. However, in the theoretical model, the permeability of the channel related to the particle migration is the main macroparameter for the risk assessment of water inrush. Under the previous assumption (iii) that the fluid mixed with fluidized particles is taken as the single-phase Newtonian flow, the permeability may be considered only depending on the porosity of the rock mass. But in fact, the permeability will be affected by the particle concentration because of particle separation, which can be understood as the filtration. For this condition, it should be more appropriate to use the two-phase flow theory, not only for the evolution of the permeability but also for the pressure dividing the two-phase flow.

Numerous water inrush disasters indicate that the water inrush through broken rock mass always occurs in a short time, but the flow still starts from low flow rate and evolves into high flow rate ineluctability. The corresponding seepage behavior is linear to nonlinear [1, 8]. Although the proposed model is more reasonable in characterizing the nonlinearity of water inrush compared with the traditional seepage model [18, 22], and the simulations also show clear differences between the two models under high pressure gradient condition, the proposed model only contains one kind of flow regime, which is not fit perfectly in theory. When the pressure gradient is not too high or the small particles in the broken rock mass have a cementing property, the evolution of water inrush may be relatively slow which will lead to a long time linear flow. In this case, the model may cause a large error due to its nonlinearity.

#### 6. Conclusions

It is of vital importance to quantitatively evaluate the flow behavior of water inrush related to the nonlinear seepage erosion in order to evaluate the risk of this kind of disaster. In this work, a coupled nonlinear flow model was established to simulate the particle migration and seepage properties of the water inrush through broken rock mass. In this model, groundwater and fluidized particles are treated as one-phase mixed fluid during the water inrush, and volume concentration is employed to describe the mixed fluid properties. Under these assumptions, both the evolution of nonlinear flow and the coupling effect between seepage erosion and permeability properties of the broken rock mass are taken into account. Theoretically, the proposed model is capable of simulating the interaction of high-velocity seepage, small particle migration, and evolution of broken rock mass structure in the process of water inrush.

The model was then simulated through a case study. The simulation results generally agree well with indoor experiments conducted by Liu et al. [8], which shows that the loss of small particles is the main reason for water inrush evolution. And thus, the water inrush through broken rock mass is a mutually promoting process of the nonlinear seepage and the erosion. Moreover, the numerical simulations indicate that the small particle migration causes the unstable characteristics of the seepage and the heterogeneity properties of the broken rock mass, which essentially induced the temporal and spatial evolution of the porosity. Consequently, this will lead to the nonlinear distribution of pressure in the broken rock mass, although the pressure gradient between the inlet and the outlet is constant. It was further confirmed that the flow behavior of water inrush is nonlinear in both time and space. Therefore, spatial evolution of the broken rock mass structure should also be paid more attention when conducting the simulation experiments. On the whole, it is concluded that the model proposed in the paper is of significant help to clearly understanding the coupling mechanism between nonlinear flow and small particle migration during the water inrush.

#### Data Availability

All data, models, and code generated or used during the study appear in the submitted article.

#### Conflicts of Interest

No conflict of interest exits in the submission of this manuscript.

#### Acknowledgments

This research is financially supported by A Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions, Joint Funds of the National Natural Science Foundation of China (U1710253), fund of Suzhou University of Science and Technology (XKQ2018005), and Science and Technology Project of Jiangsu Construction System (2018ZD033).