Abstract
For water inrush induced by fracture network flow, the critical velocity of the incipient motion of sand particles was obtained, and the functional relation between critical velocity and particle size was established through a series of tests on the nonlinear flow characteristics of a filling fracture network. The influence of the particle size distribution, hydrodynamic force, and geometric features of the fracture network on the characteristics of particle loss; distribution laws; and water-sand, two-phase migration was also explored. Moreover, the interactions amongst water, movable particles, the surface of the skeleton, and fracture walls, and the formation mechanism of the flow channel were qualitatively analyzed. In addition, the change rules of the mass loss characteristics and porosity of the samples with time were tested successfully. The calculation methods of the permeability and non-Darcy factor of the filling fracture network were also determined.
1. Introduction
Filling fracture networks are common in practical projects. Most materials used for filling fracture networks are particle accumulation bodies (e.g., breccia, sand, travertine formed by tectonism, weathering, or the action of groundwater, as well as fractured rock masses from the destruction of the protolith by earthquakes, mining disturbances, and other factors) [1–3]. During underground construction, fracture zones develop continuously and achieve a hydraulic connection between the nearby aquifer and the working face due to actions, such as disturbance and unloading, which may lead to nonlinear flow in the fracture network. Under the constant scouring action of the current, materials with small particle sizes filled in the fracture network are continuously lost; moreover, the permeability of the fracture network increases, subsequently inducing water inrush accidents. A schematic of water inrush is shown in Figure 1 [4]. Therefore, the nonlinear flow characteristics of filling fracture networks must be studied to reveal the mechanism of water inrush catastrophe induced by fracture network flow and the water-sand, two-phase migration law, and provide scientific basis for the prediction and prevention of water inrush disasters during underground excavation.

By conducting numerous experiments and theoretical studies, Tian et al. [5] and Su and Zhan [6] concluded that the fracture width, porosity, and characteristic particle size of filling materials considerably affect the permeability of filling fractures. On this basis, many scholars have performed extended studies on filling fractures. Tao [7] derived a general formula of sectional average velocity, which can not only describe nonfilling fracture flow but also reflect filling fracture flow. On this basis, Yu and Tao [8] derived the expression of coefficients A and B in the Forchheimer flow equation through theoretical analysis. This expression is suitable for filling and nonfilling fracture flow and for describing linear and nonlinear flow. Liu et al. [9] held that particulate packing matter is loose and that small particles easily get lost under high hydraulic gradient conditions because their physical and mechanical properties considerably change after water saturation. Therefore, the flow state of water in the filling fracture also needs to be verified through a hydraulic test. Many scholars have also studied the influence of the effects of coupled seepage and stress and the hydraulic effects of the filling medium on the flow characteristics of filled fractures [10, 11]. The discharge capacity of filling fractures is not only associated with fracture opening and the packing porosity and particle size of filling materials, it is also affected by the particle size distribution and arrangement of the filling materials. However, the influence of particle size distribution on nonlinear flow parameters is rarely reported in the field of porous media.
At present, studies on water-sand, two-phase flow have focused on sand movement in rivers; these studies have also discussed issues, such as the incipient motion characteristics of sand [12], the transport of unbalanced sand [13], the motion characteristics of sand-carrying current [14], the motion law of suspended load [15], and the motion law of bed load [16], and made remarkable progress. The critical incipient condition of sand is the first problem encountered in the study of engineering sand. Drag force is the main driving force of sand motion [17]. However, in consideration of the correlation between the velocity field and shear field, an increasing number of scholars have shown interest in deducing incipient velocity to facilitate practical application without considering the incipient drag force of sand [18]. The incipient law of sand has become complicated, and no unified formula has been developed to describe the incipient velocity of particles due to the influence of many factors, such as the position, arrangement, and relative exposure of sand particles in riverbeds [12, 13]. Moreover, under the influence of the randomness of water flow movement, sand movement in water flow, and sand particle position at the bed surface, the transition from the resting phase to the incipient motion of sand particles is random and is a typical probability event [19]. The water-sand, two-phase flow in rock fractures differs from sand movement in rivers. The content of sand in rivers is relatively low. Particle movement is mainly affected by hydrodynamic conditions and is irrelevant to collisions between particles [20]. Wang et al. [21] believed that the interactions between particles are of four kinds: static support, relative slipping, collision, and diffusion. And, the particle collision and particle-wall collision cannot be neglected because the sand concentration in fractures is high. Therefore, the water-sand, two-phase flow is a particle flow considering the interstitial fluid action.
For filling fracture networks, the mixing law of the sand at the intersection is also worthy of attention. At present, the main mixing models for water and sand at the intersection of two-dimensional fractures are the streamline path and complete mixing models [22, 23]. Wilson and Witherspoon [24] used crisscrossing circular tubes of equal diameter to simulate intersected fracture flow and proposed applicable conditions for the two mixing modes at the fracture intersection. Robinson and Gale [25] found by numerical simulation that the intersection angle of fractures has a certain influence on the mixing characteristics of fluids. The mixing degree of the solutes reaches the maximum when the intersection angle is within 67.5°–112.5° and the minimum when the angle is any other value. Park et al. [26] studied the mixing characteristics of sand at fracture intersections from two perspectives: local flow conditions and geometric parameters of fracture intersections. The results showed that geometric parameters remarkably affect the mixing characteristics. Although most of the studies above were based on the cubic law and rarely considered the Navier-Stokes flow, these studies still laid a solid foundation for understanding the characteristics of water-sand, two-phase nonlinear flow in fracture networks. In addition, a team from the China University of Mining and Technology conducted an in-depth study on the variable mass flow of fractured rock masses [27, 28].
Water-sand flow through rock fractures is a liquid-solid, two-phase flow in a broad sense. Previous studies could not observe the entire process of particle loss due to the limitations of the test equipment. Studies on critical incipient velocity, the migration law of particles, the formation process of dominant channels, and the characteristics of water-sand, two-phase migration under high hydraulic gradient are lacking because the accurate separation phase measurement of water and sand could not be achieved. Therefore, this study used a self-developed test device to perform a nonlinear flow test on filling fracture networks, considering the characteristics of particle loss; moreover, it analyzed the law of incipient motion, the migration of filling particles in fracture networks, and the porosity change rules caused by particle loss, thereby revealing the catastrophe mechanism of water inrush caused by flow in filling fracture networks.
2. Experimental Details
2.1. Experimental Equipment
Through extensive practice and testing, a set of planar experimental system for hydraulic sediment transport that satisfied the experiment requirements of high concentration sand-water mixture flow motion at higher pressure and greater flow rate is designed. The matching water-sand separation phase velocimeter solved the problem of complete solid-liquid phase separation and could measure the water-sand separation phase velocity in a real-time and accurate manner. The system mainly consists of seven parts: a water supply system, a fracture network flow model, a rotary table, a water-sand separation phase velocimeter, a back pressure device, an online saturation device, and a data acquisition system (see Figure 2).(1)Power supply system: the water supply modes in this study are divided into two kinds: constant pressure and constant current. Constant-pressure water supply can provide the maximum constant upstream water pressure of the sample, i.e., 0.7 MPa (70 m water column); moreover, the generated hydraulic gradient exceeds most of the test equipment, with a maximum water flow of 48,000 ml/min. Constant-current water supply can provide a steady current of 20–24,000 ml/min. During the experiment, a 0.005 MPa water supply pressure, provided by the constant-pressure pump, is increased in each step until the water pressure reaches 0.7 MPa. Water-sand, two-phase flow is a test that occurs at a high flow rate, using more than 1,000 ml of liquid. At present, for the core flow test in petroleum engineering, the filling pump (volume of single pump: 1,500 ml) with the largest volume is the ISCO pump manufactured in the United States. The volume of the double-cylinder filling pump used in this study is 24,000 ml and can meet the requirements of most tests on rock fracture flow.(2)Fracture network flow model: the dimensions of the mountable part are 600 mm ∗ 300 mm ∗ 40 mm. Slots for perforated plates are provided at both ends of the sample, which can be fixed by perforated plates. The diameter of the pipeline connecting the outlet pump and the sample is generally smaller than the size of the sample; thus, the velocity of the liquid flowing in the pipeline is greater than the flow velocity of the liquid in the sample. If the infusion tube is directly connected to the sample, then the velocity at the upstream end of the sample on the cross section becomes uneven. Therefore, an inlet mechanism is required to realize the uniform dispersion of high-speed water flow on the cross section of the sample. The inlet mechanism is composed of equal-diameter inlet pipes, a diverging section, and a perforated plate. Sand and other particle materials are typical energy-dissipating materials; thus, the outlet mechanism is an outlet that contains a converging channel. Under stable flow, the water-sand mixed fluid flowing out of the sample improves the sand-carrying flow capacity of water and prevents sand from accumulating at the outlet of the sample, resulting in the flow of the water-sand mixed fluid through the filling fracture network; this result is due to the continuous contraction of the channel section and the increasing velocity.(3)Water-sand separation phase velocimeter: at present, solid-liquid phase separation can be performed in two ways: one is to separate water from solid-particle materials; the other is to separate solid particles from water. Although neither of the two technical means can completely separate the two phases, the velocity of sand and water can be measured separately by mass replacement. Therefore, the real-time measurement of velocity was replaced by the real-time weighing of mass in this study. The water-sand separation phase velocimeter is shown in Figure 2(b). It can measure the phase separation of water-sand mixed fluid with a flow of up to 1,200 ml/min and a sand volume of over 50%; in addition, the measurement accuracy can reach up to 95%. Xu et al. [29] provide other technical measures to reduce measurement errors and improve the accuracy of phase separation velocity measurement (e.g., calibration of the volume of the constant-volume separator, a method to reduce the effect of water-sand mixture on electronic balance, a method to separate the suspended particles from the spilled fluid completely, and a method to eliminate the additional force of the connection).(4)Back pressure device: if the outlet is directly connected to the atmosphere, then airwaves cause flow instability. In such a case, even the use of high-precision sensors cannot improve the accuracy of pressure measurement. Therefore, a certain back pressure must be applied at the outlet to make the pressure slightly higher than atmospheric pressure (e.g. 3 kPa) and eliminate the fluctuations in the outlet velocity at the test section. However, for water-sand mixed fluid, the sand cannot pass through the back pressure valve. The valve body is damaged even if it can pass through the back pressure valve. Therefore, the back pressure for the flow of water-sand mixed fluid is generated by the height of the liquid column.(5)Online saturation device: according to the principle of the water-sand separation phase velocimeter, if there are bubbles on the surface of sand particles, then the volume of replaced water becomes larger than that of sand. On the other hand, the mothball bubbles in the samples, the pipelines, and the connecting valves have great influence on the water flow and sand migration. To ensure test precision, the saturation of the sample may be improved to more than 95% by displacing the bubbles on the surface of the wetted fracture and equipping the installed sample with an online saturation device.(6)Data acquisition system: the model of the pressure sensor is Q/JL015-2009; the range is 0–1.38 MPa, and the precision is 1%. In accordance with the water flow at the outlet, balances produced by Sartorius Stedim Biotech GmbH in Germany (i.e. TE 2101-L, with a range of 2,100 g and a precision of 0.1 g; TE 12000-L, with a range of 12,000 g and a precision of 1 g; with a range of 0–150 kg and a precision of 0.01 kg) were used for measurement. The sampling interval was 3 s. The flow was recorded and observed in real time through the data acquisition system.

The design of the perforated plate should meet the following principles: (1) Large and small holes should be arranged on the stainless-steel base plate. (2) The total surface porosity should be equal to or close to that of the coarse skeleton of the filling medium. (3) The diameter of large holes should be greater than or equal to six times the maximum diameter of sand particles that migrate, to ensure that the sand can flow out of the sample. Meanwhile, the heterogeneity of particles should be considered. The diameter of large holes should be equal to or smaller than 1/3-1/2 of the size of coarse particles to prevent coarse particles from moving [30–32]. (4) On the premise of ensuring strength, the stainless-steel perforated plate should be as thin as possible to minimize the additional pressure loss caused by the perforated plate. The perforated plate designed in accordance with the above principles is shown in Figure 2.
2.2. Experimental Materials and Scheme
Studying the rule of particle loss and the formation process of the flow channel from the perspective of particle size distribution is the development trend of future studies [33]. White quartz sand with particle sizes of 2.0–2.36, 2.36–4.75, and 4.75–9.0 mm was configured to form skeleton particles in different proportions. Yellow aeolian sand with particle sizes of 0.075–0.15, 0.15–0.3, 0.3–0.6, 0.6–1.0, and 1.00–2.0 mm was configured to form movable particles in different proportions. The test materials are shown in Figure 3. The particle size distribution curve for the sample is shown in Figure 4. Given that the maximum size of sand particles was 9.0 mm, the fracture network with a fracture opening of 40 mm was selected. The filled sample is shown in Figure 5. The test plan is shown in Tables 1 and 2.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)


2.3. Testing Procedure
As shown in Figure 6, the test procedures for the nonlinear flow characteristics of the filling fracture network include five steps: sample preparation and filling, online saturation, system debugging, data acquisition, and unloading.(1)Sample preparation and filling: the sand required for each particle size range was weighed in accordance with the specific particle size distribution and then properly recorded. An appropriate amount of water was added to wet the particle surface. The particles were mixed well and filled to the fracture network of organic glass. During filling, particles should be filled in layers to ensure the uniformity of the sample. Particles of the same mass should be filled each time and then compacted by layers.(2)Sample saturation: the saturation method and precautions were the same as before.(3)System debugging: before the start of the test, the following conditions should be checked: whether the pipeline connection is correct, whether the valve is open, and whether the pressure sensor, electronic balance, and data collector are operating normally; in addition, the camera should be set up.(4)Data acquisition: the total length of the sample was 60 cm. Five key data acquisition areas and phenomenon observation areas should be successively arranged at 0, 15, 30, 45, and 60 cm. Data measurement software was used to monitor and record data, including water-sand separation phase velocity and pressure distribution; photos or videos should be taken and recorded, respectively.(5)Unloading: after the test, the outlet stop valve of the plunger pump should be closed, and the pipeline interface between the plane model, the plunger pump, and the water-sand separation phase velocimeter should be disconnected. The bolts, hollow steel cover plate, organic glass plate, and other parts should be removed successively, and the sample partition should be removed, screened, and weighted after drying.

2.4. Calculation of Nonlinear Flow Parameters
2.4.1. Particle Loss
As shown in Figure 2(b), initially the water-sand mixed fluid enters vessel 1 filled with water, then the sand particles sink down to the bottom, and pure water falls into vessel 2. Electronic balance 1 reads the difference in the masses of the aeolian sand and the water it replaces, and electronic balance 2 reads the mass of the water with the same volume as the water-sand mixed fluid. And, the volume of the aeolian sand is equal to the volume of the displaced water, thereforewhere is the record of electronic balance 1, is the record of electronic balance 2, is the mass of the aeolian sand flowing out of the sample, is the mass of the water displaced by the aeolian sand, is the mass of the water flowing through the sample, is the volume of the water displaced by the aeolian sand, is the volume of the aeolian sand, is the volume of the water flowing through the sample, is the density of aeolian sand, and is the density of water.
Equations (1) and (2), lead to
The mass of aeolian sand and water flowing through the sample at a certain interval of time (sample time) can be calculated by the following formulas, respectively:
Therefore, the total mass of the aeolian sand and water at (i = 1, 2, 3, …, n) can, respectively, be calculated by
The volume concentration of the fluidised particles in the flow, can be calculated aswhere cti is the volume concentration of the fluidised particles in the flow, and its physical significance is the ratio of particle volume and fluid volume flowing through unit volume sample in unit time.
2.4.2. Porosity Evolution
The initial porosity of the sample can be written aswhere is the total volume of the fracture network, is the volume of the aeolian sand filling in the fracture network, is the volume of the quartz sand filling in the fracture network, is the mass of the aeolian sand filling in the fracture network, is the mass of the quartz sand filling in the fracture network, and is the density of the quartz sand.
According to the total mass of particle loss, the porosity ϕti of the sample can be obtained as follows:
The change rate ϕ′ of porosity in time can be redefined as:
2.4.3. Permeability and Non-Darcy Factor
Yu and Tao [8] obtained the expressions of the permeability kff and non-Darcy factor βff of the filling fracture through theoretical derivation:where is the permeability of filling fracture network, b is the fracture aperture, d20 is the grain size of 20% passing through the corresponding sieve by weight, is the non-Darcy factor of filling fracture network, A and B are dimensionless coefficients.
When the motion state of the fluid in the filling fracture is laminar flow and b/d20 > 10, A = 5.1 [6], then
The following relation exists between the permeability of the spherical particle through porous medium and the size and porosity of particles [6, 34]:where is the permeability of porous media.
It can be obtained by comparing formulae (11) and (12).
Formula (14) shows that the difference in the permeability of the filling fracture network and porous media is mainly decided by b/d20, i.e., the wall effect on the flow characteristics in porous media. When b/d20 is great enough, the wall effect can be ignored. The specific range needs to be further determined through tests.
3. Experimental Results and Discussion
3.1. Migration Law of Particles
Figure 7 shows the corresponding relation between the hydraulic gradient and the velocity of the sample under the conditions of different particle size distributions and different fracture intersection angles. In this test, the upstream water pressure was gradually increased. The process of particle loss (flow erosion) can be divided into the following four stages according to the figure: ① Incipient motion of particles: when the velocity does not reach the starting critical velocity of aeolian sand, the aeolian sand particles are motionless, the specimen porosity remains constant, and the velocity increases smoothly with the increase of the hydraulic gradient (See the red data points). Generally, the incipient motion of particles can be divided into three states: particles about to move, a few particles moving, and numerous particles moving [34]. Drag force is the main driving force of sand motion and increases as the water velocity increases. When it increases to a certain value, it can overcome the effective gravity of particles and the bite force between particles; it can also break the network of forces formed between particles through a force chain (or arch). The incipient motion of particles starts. At this moment, the hydraulic gradient is the critical hydraulic gradient of particle loss. ② Particle migration: before particle migration, a marked increase in flow velocity occurs; such an increase is caused by the increase in porosity due to the volumetric dilatation of aeolian sand. Moreover, the larger the diameter of the pore channel, the smaller the fracture network angle, the more the particles with a particle size of 0.075–0.15 mm amongst the movable particles, and the smaller the hydraulic gradient in the case of mutation. Incipient particles tend to migrate along the direction of water flow. When the diameter of the adjacent pore is larger than the size of incipient particles, these particles enter this pore; when the diameter of the adjacent pore is smaller than the size of the incipient particles, these particles can no longer migrate and become part of the skeleton. So, the volatility of the data points in this period is relatively large (see the green data points). ③ Intersection of local channels: under the continuous action of water drag force, the finest particles are constantly washed away, forming intersecting channels with a small pore diameter and generating a preferential flow. Even though these channels are blocked by particles with greater particle size, which seemingly blocks the migration of particles, the blocking particles deflect the current flowing through that point, taking away the surrounding fine particles and forming a new flow channel. In such a case, a marked increase in velocity occurs (see the green data points). ④ Formation of stable flow channel: under the alternate action of blocking intersections, the following occurs: small channels are intersected to form large ones, partial particles are rearranged, areas with locally dense particles are released, sand particles escape on a large scale, and the sand gushing area continues to erode upstream until all fine particles in the sample are washed out to form a stable flow channel (see the blue data points).

(a)

(b)

(c)

(d)

(e)
3.2. Relation between Critical Incipient Velocity and Size of Particles
The flow, migration, and loss of movable fine particles occur along with the pore water. Water-sand, two-phase interactions occur throughout this process. The kinetic energy of water is continuously converted into the kinetic energy of aeolian sand in this process. The critical incipient conditions for aeolian sand of different particle sizes are listed in Table 3.
Under the effect of adhesion stress and water pulsation, the incipient motion of fine particles is usually in the form of particle swarm. These particles still move with water flow in the form of individual particles after incipient motion. However, the incipient motion of coarse particles is in the form of individual particles at any time [35]. Given that flow erosion is the carrier of fine particles through pores or (and) fracture channels, the forces on individual particles must be analyzed. Figure 8 shows the relation between the critical incipient velocity of particles and . As shown in the figure, the critical incipient velocity of particles increases with the increase in particle size. The equation can be obtained by data fitting: .

3.3. Formation Mechanism of the Flow Channel
After the test, the remaining sand sample in the plane model was successively removed four times and then dried for particle analysis using a high-frequency sieve shaker. Figure 9 shows the changes in the particle size distribution of the sample before and after the test. FF1–FF5 represents the initial particle size distribution of the sample. The subscripts represent the samples at different heights after the test. From the top of the sample, 0–15 cm is denoted as the first layer (FF11–FF51), 15–30 cm is denoted as the second layer (FF12–FF52), 30–45 cm is denoted as the third layer (FF13–FF53), and 45–60 cm is denoted as the fourth layer (FF14–FF54). As shown in Figures 9 and 10, the particles greater than 2.36 mm formed a skeleton of the filling medium and no particle loss occurred; thus, the curve of particle size distribution for those particles remains unchanged before and after the test, whereas a significant difference is observed in the curve for particles smaller than 2.36 mm before and after the test and at different depths due to the loss of internal movable particles. The content of movable particles in the five groups of samples is approximately 27%, and the loss percentage is approximately 13%–18%, accounting for 50%–70% of the content of movable particles; the loss is relatively sufficient. The closer to the downstream of the sample, the less the particles are lost. This result is not due to the small number of migrating particles in this layer but due to the long distance of particle migration in the lower part of the sample and the small flow force. When flowing through the upper pore space, retention and filling effects tend to happen due to size limit, thus influencing the loss amount of upper filling particles. These results are consistent with the conclusion in the study conducted by Zhou et al. [36] and Yao et al. [37] on the microscopic mechanism of sand piping.

(a)

(b)

(c)

(d)

(e)

3.4. Spatiotemporal Evolution Law of Porosity
According to Figure 11, in the first 40 s of the test, the porosity of each sample showed an accelerated increase, indicating that the mass loss rate of particles increases continuously during this period; when t = 40–80 s, the porosity of each sample showed decelerated increase, indicating that the number of particles getting lost decreases; after t = 80 s, porosity tended to stabilize, indicating that no particle is lost. Comparing samples FF6, FF7, and FF8 revealed that the greater the nonuniform coefficient and the curvature coefficient, the greater the increase in porosity, which is consistent with the phenomenon observed in [36]. The reason for this is that in the case of the same initial porosity, fewer pores are formed by skeleton particles in samples with a large nonuniform coefficient; the larger the pore diameter, the smaller the resistance to the sand-carrying flow. Therefore, the particles get rapidly lost, and the lost particles account for a large proportion of the content of movable particles. Although samples with a small nonuniform coefficient have high content of fine particles, fine particles are not easily lost due to their small pore diameter.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)
As shown in Figure 11 and Table 4, comparing samples F11–F14 reveals that the smaller the size of migrating particles, the faster and the greater the increase in porosity. Comparing samples FF6 and FF10 reveals that the greater the angle between the fracture and the vertical direction, the smaller the increase in porosity. In sample FF6, particles turned by 60° at the intersections of the fracture network during migration. However, the particles in sample FF10 turned by 90°; thus, more kinetic energy and relatively less mass were lost. Comparing samples FF6 and FF7 reveals that the higher the proportion of particles with a size of 0.075–0.15 mm, the greater the increase in porosity within the same period. The reason is that the loss of fine particles creates expansion space for the movement of large particles and leads to the rapid loss of large particles. Nevertheless, the particle size distribution of movable particles does not affect the final loss amount.
The movable particles in the filling medium constantly get lost due to the constant erosion of the current, directly leading to the change in the porosity of the filling fracture network. Sakthivadivel [38] once proposed the evolution equation of porosity of loose media under current erosion:where λ is the erosion coefficient, the dimension of which is the reciprocal of the length and can be obtained by laboratory test. The method to obtain parameter λ is little complicated. It is assumed that the parameter is only related to the ratio and initial porosity of samples, that is, the parameter can be determined by testing multiple groups of samples and using regression analysis at the same ratio and initial porosity.
The applicability of equation (15) for these experimental results has been discussed by many scholars, such as Yao [39], Ma et al. [40], Yang [41], and Shi and Yang [42]. The particle loss in filling fracture networks is different from the particle migration in loose medium. Not all the filling particles in the fracture network can be fluidized. Coarse particles keep still as the skeleton. Only a small part of small filling particles is fluidized. Therefore, the porosity of the filling fracture network has a maximum value . The porosity evolution equation can be corrected as follows [39]:
Equation (16) is true only in the case of particle loss, i.e., the flow velocity of water reaches the critical incipient velocity of particles . Therefore [43],
As shown by the fitting results in Figure 11, the overall trend of the porosity change rate of the nine groups of samples was consistent, i.e., increasing to the extreme and then decreasing slowly and then tending to 0. The difference is that each sample fluctuates to different extents at different time points. This result can be attributed to the alternate action of blocking intersections during the migration of aeolian sand particles. The fitting results of porosity evolution equation (16) are consistent with the test results, indicating that equation (16) can accurately describe the evolution law of the porosity of the filling fracture network under current erosion.
Significant differences in the erosion coefficient λ were observed amongst the nine groups of samples, indicating that the geometry of the pore channel, the particle size distribution, size of movable particles, and the angle of the fracture network have a great influence on the value of the coefficient λ. Figure 12 shows the correspondence between erosion coefficient and the average particle size of migrating particles. According to the figure, the relation between the erosion coefficient and the average particle size of migrating particles is a power function. The erosion coefficient decreases nonlinearly with the increase in average particle size.

3.5. Water-Sand, Two-Phase Migration Characteristics
According to Figure 13, during particle migration, the volumetric flow rate of water is significantly greater than that of sand, and the water-sand separation phase velocity and concentration show fluctuations. The transient fluctuations are caused by the alternate action of the migration and blocking of aeolian sand particles; moreover, with the constant migration and loss of aeolian sand particles, the overall macro average porosity of the filling fracture network keeps increasing. However, the local porosity decreases when migrating particles block the pore throat, resulting in the temporary decrease in the mean velocity of the current; in addition, the fluid velocity increases rapidly again when the blocking particles are washed away. Therefore, the water-sand separation phase velocity and concentration show considerable fluctuations. By contrast, the water-sand separation phase velocity in sample FF6 can maintain a relatively smooth motion. Ma et al. [40] indicated that this behaviour may be related with the rearrangement of the skeleton. By conducting numerous experimental studies, Yang et al. [44] concluded that, at the beginning of the test, water is the dominant phase. As the water flow velocity increases, a large number of aeolian sand particles flow with the water flow, and the aeolian sand becomes the dominant phase; meanwhile, a large amount of movable aeolian sand particles occupies the flow pathways, leading to the reduction of water flow velocity. As the flow rate of water decrease, the particle-carrying ability decreases accordingly. However, this may be related to the rheological properties of the mixed fluid, because extremely small changes in the concentration of aeolian sand can lead to large variations in velocity [45]. Therefore, these fluctuations need to be further studied from the perspective of meso-mechanics.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)
The test results above show that the sand-carrying flow under constant pressure gradient is unstable. With the appearance of sand-carrying flow, the velocity of water phase and sand phase increases continuously, and the violent sand gushing is caused by the development of sand-carrying flow. The reason for the instability of sand-carrying flow is as follows: with the loss of sand particles, the porosity keeps increasing, but the resistance to sand-carrying flow decreases, thus creating conditions conducive to the acceleration of water flow. Consequently, larger sand particles migrate until the sand particles are completely lost and the porosity stabilizes; this migration is an interactive mass change process. Therefore, the loss of sand particles and the acceleration of water flow are mutually reinforcing positive feedback processes and are the key factors to inducing flow mutation in fracture networks.
3.6. Influence of Filling Medium on the Flow Characteristics of the Fracture Network
Nonlinear flow tests were conducted on porous media with different particle size distributions; see Reference [4] for details on these tests. And the sand samples numbered 2, 4, 5, 6, and 10 in the literature were filled in the fracture network with a fracture aperture of 10 mm and an angle of 30°, respectively. A stainless-steel filter screen with a sieve mesh aperture of 0.074 mm was fastened to the middle of the porous plates. Because the smallest sand grain diameter was 0.075 mm, the filter screen guaranteed that the sand was not washed away. To analyze the influence of the particle size distribution of the filling medium on the flow characteristics of the filling fracture network quantitatively, the relative permeability error is defined as follows:where is the relative permeability error.
Figure 14 shows the results of the error analysis under different b/d20 conditions. When b/d20 = 7, δ = 16%; when b/d20 = 33, δ = 3%. The results show that with the continuous increase in b/d20, the relative error between the permeability of the filling fracture network and that of the filling medium decreases gradually; when b/d20 > 33, the permeability of the filling medium is almost equal to that of the filling fracture network. Meanwhile, given that a power function relation exists between the permeability and the non-Darcy factor, the nonlinear flow parameters of porous media can be used to characterize the permeability of filling fractures. The results also show that for water, the permeability of the filling fracture network is lower than that of the filling material itself; this result is consistent with the conclusion of Reference [5].

4. Conclusion
In the test herein, a filling fracture network was used as the study object. A self-developed, plane hydraulic sand transport test system was used to study the influence of the particle size distribution of the filling medium, the intersection angle of fracture networks, and other factors on flow mutation. The following conclusions were drawn:(1)The process of particle migration can be divided into four stages: ① incipient motion of particles, ② particle migration, ③ intersection of local channels, and ④ formation of stable flow channels. The formation of intersecting channels has randomness, abruptness, and intermittency. The alternate action of blocking intersections redistributes the hydraulic conditions of the entire flow field, making the mechanism of water inrush disaster complicated. The water inrush in the filling fracture network, which contains a complex water-sand, two-phase interaction process, is the macroscopic embodiment of the combinations of countless microscopic processes.(2)Under laboratory conditions, the critical incipient velocity of aeolian sand with a particle size of 0.075–0.15 mm is 1.96–2.35 mm/s, and the corresponding hydraulic gradient is approximately 3.40–5.90; the critical incipient velocity of aeolian sand with a particle size of 0.15–0.3 mm is 2.21–2.70 mm/s, and the corresponding hydraulic gradient is approximately 4.25–6.80; the critical incipient velocity of aeolian sand with a particle size of 0.3–0.6 mm is 2.28–3.20 mm/s, and the corresponding hydraulic gradient is approximately 5.10–7.60. The critical incipient velocity of particles increases with the increase in particle size. The fitting equation can be obtained by analyzing the critical incipient velocity and size of particles. This equation is only applicable to the model in this study.(3)The process of flow erosion is a multifield, multiphase, coupled, highly nonlinear, and dynamic process. Under the action of the flow force, fluidized particles constantly migrate and get lost, resulting in increasing porosity and creating conditions conducive to the acceleration of water flow. Thus, particle migration and water flow acceleration are a mutually reinforcing positive feedback process. The fitting results of the porosity evolution equation are consistent with the test results, indicating that this equation can accurately describe the evolution law of the porosity of the filling fracture network under current erosion. Meanwhile, the relation between the erosion coefficient and the average particle size of migrating particles is a power function. The erosion coefficient decreases nonlinearly with the increase in average particle size.(4)The influence of factors, such as the particle size distribution of the skeleton, the particle size distribution of movable particles, the intersection angle of fractures, and different upstream water pressures on the flow characteristics of the filling fracture network, was investigated. Results showed that in the case of the same particle size distribution of the skeleton, the higher the proportion of the medium and small particles in movable particles, the greater the amount of particle loss within the same period; however, the particle size distribution of the movable particle itself does not affect the final loss amount. When other conditions are the same, the greater the angle between the fracture and the vertical direction, the smaller the particle size; the higher the upstream water pressure, the faster the particle loss and the relatively more stable the water-sand separation phase velocity.(5)With the continuous increase in b/d20, the relative error between the permeability of the filling fracture network and that of porous media decreases gradually; when b/d20 > 33, the influence of the wall effect of fracture networks can be ignored. That is, the flow parameters of porous media can be used to characterize the permeability of filling fracture networks.(6)It is concluded that the migration and loss of fine particles in the fracture network are the origin of the water inrush disaster. From the perspective of water inrush prevention, on the one hand, the filling particles can be solidified by pumping coagulants (such as concrete grout) inside the fracture network to prevent the occurrence of particle fluidization. On the other hand, the relationship of water power between confined aquifer and working face can also be blocked by plugging measures, such as the advanced pipe grouting and the reserving of waterproof surrounding rock.
Data Availability
The underlying data supporting the results of our study are included within the manuscript.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The current thesis was supported by the State Key Program of National Natural Science Foundation of China (U1710253), the National Natural Science Foundation of China (Grant No. 51574059), and the Youth Foundation of University of Science and Technology Liaoning (Grant No. 2020QN10).