#### Abstract

With the increase in mining depth, potential water inrush in the limestone-fractured aquifers of Taiyuan formation and Ordovician formation becomes seriously threatened to the mining process. Grouting reinforcement, as a common prevention method, has been widely used in mine water inrush prevention. In this paper, a model of “seepage resistance” is proposed for the slurry diffusion in a single fracture and the theoretical expression on the seepage resistance corresponding to slurry with different flow patterns is conducted. The total seepage resistance model of water-displacing diffusion by slurry in series and parallel fractures is established. With the Lagrange interfacial tracing method, the control equations of water-displacing diffusion by slurry in fracture networks are derived, which get mutually validated with the theoretical results.

#### 1. Introduction

The total coal reserves are abundant in China which accounted for 94% of total primary energy sources and about one-seventh of the world’s total proven coal reserves. Although the proportion of total coal consumption in China’s total disposable energy consumption has decreased in recent years, the current situation of coal-based energy structure in China will not change in the short term [1]. With the exhaustion of mining of shallow coal resources, the mining of deep coal resources has become normal. However, potential water inrush in the limestone-fractured aquifers of Taiyuan formation and Ordovician formation becomes seriously threatened to the process of mining deep coal resources [2]. According to incomplete statistics, the reserves of the North China type coal field threatened by the fractured aquifers of Taiyuan formation and Ordovician formation in the mining process are as high as 10 billion tons [3]. For the mining of coal resources under loose aquifers, the gangue waste filling-mining method is usually used to reduce the development height of water-conducting fracture zone, so as to improve the upper mining limit and reduce the waste of coal resources. At the same time, in order to liberate a large number of coal resources threatened by confined fractured aquifers, the grouting reconstruction technology for the fractured aquifers is mainly adopted to prevent water inrush from the floor of coalbed in the mining process [4–7], which can not only displace the water in the fractured aquifer but also greatly increase the strength of the fractured rock masses after the fractures are blocked; thus, the purpose of turning the fractured aquifer into a relatively watertight layer is achieved.

In recent years, according to research on the diffusion mechanism of slurry in fractured rock masses by theoretical analysis, numerical simulation, laboratory test, and other research methods, many scholars at home and abroad have achieved many fruitful research results [8–17]. El Tani and Stille [18] derived the expression of the relationship between grouting time and diffusion range of slurry when Newtonian and Bingham slurry diffused in canal fractures and parallel plate fractures under the conditions of constant pressure grouting and constant flow grouting, which can provide a reference for grouting design of fractured rock mass. Mohajerani et al. [19] designed the slurry pressure display calculation program (EGFP) for calculating the diffusion distance of slurry with different flow patterns (Newtonian slurry, Bingham slurry, and power-law slurry) in fractured rock mass. According to the analysis of the calculation examples using EGFP, it can be seen that the largest advantage of this method is high calculation efficiency. Liu and Sun [20] developed a grouting simulator considering the coupling effect between the slurry and rock, to calculate the fracture initiation, propagation, deformation of solids, and the flow process of slurry in the fracture. This method takes into account not only the fluid-solid coupling effect in the flow process of slurry but also factors such as hardening of slurry and evolution of fracture grid. Jin et al. [21] adopted the method of orthogonal experiment to analyze the diffusion effect of chemical slurry in the fracture under the condition of dynamic water; the results show that the water-plugging rate is positively correlated with the fracture roughness coefficient JRC but negatively correlated with slurry gel time, initial water flow rate, and fracture opening, and the curve of water flow rate and seepage pressure can be divided into three forms. The research results can provide certain reference for fracture grouting in rock mass. Hao et al. [22] studied the diffusion mechanism of self-expanding polyurethane slurry in fractures by laboratory tests and established a density model of the diffusion of polyurethane slurry in fractures by fitting the experimental data obtained. Chen et al. [13] used the coupling of finite element method (FEM) and volume of fluid (VOF) method to simulate the fracture of grouting process. In this method, the coupling effect of water and slurry two-phase fluid flow, stress distribution, and mesh cell damage is considered; in addition, the hardening process of slurry in soil is described by Young’s modulus and viscosity which change with time. Finally, some examples of fracture grouting are simulated, and the effectiveness of the method is verified under the premise that the results obtained are similar to the experimental results.

The previous research mainly focuses on exploring the diffusion mechanism of slurry in non-water-filled fractured rock mass, but there are few research studies on the two-phase flow problem of water-displacing diffusion of slurry in the water-filled fractured structure. Since the existence of water in the fracture and the size of water pressure will change the diffusion state of slurry and the herringbone fracture structure is the basic component unit of the relatively complex discrete fracture network, it is necessary to further explore the diffusion mechanism of slurry in series and parallel structures such as herringbone fractures, so as to lay a foundation for the research on the water-displacing diffusion mechanism of slurry in complex fracture networks and finally achieve the purpose of choosing reasonable grouting parameters and ensuring the safety of mining.

#### 2. Analysis of the Water-Displacing Mechanism in a Single Fracture

Exploring the mechanism of water-displacing diffusion of slurry in a single fracture is the basis of studying water-displacing diffusion of slurry in fracture networks, and also, the pure cement slurry often used in the project site can be divided into power-law slurry (0.5 < W/C < 0.7), Bingham slurry (0.7 < W/C < 1.0), and Newtonian slurry (W/C > 2.0) [23]. In this section, theoretical analysis will be adopted to systematically analyze the changing law of the inflow flow rate along with the grouting time in the process of water-displacing diffusion of slurry in a single fracture with different flow patterns under constant pressure grouting.

##### 2.1. Basic Assumptions

In the process of water-displacing diffusion of slurry in single fracture, the basic assumption is as follows:(1)The slurry is an isotropic homogeneous fluid and incompressible(2)The slurry has laminar flow in the process of water-displacing diffusion(3)Slurry and water are incompatible, and water is completely displaced by slurry(4)The fracture is saturated with water at the initial stage, and the slurry and water only diffuse in the fracture(5)The velocity of slurry at the upper and lower surface of the fracture is 0, and the trace length of the fracture is the limit water-displacing diffusion distance of slurry

##### 2.2. Theoretical Model

The schematic diagram of the water-displacing diffusion of each flow pattern of the slurry and the stress on the corresponding slurry microelement in a single fracture are shown in Figure 1. As shown in Figure 1, the equivalent hydraulic gap width of the fracture is *b*, the trace length is *L*, the grouting pressure is *P*_{1}, and the hydrostatic pressure is .

According to the force balance of the slurry microelement along the *x*-axis, it can be seen that

The rheological equation of Bingham slurry can be expressed aswhere is the plastic viscosity of Bingham slurry; *τ* is the shear stress; and *τ*_{0} is the yield stress of Bingham slurry.

From the velocity boundary conditions of the flow core and nonflow core areas, it can be seen that

By introducing (4) into (3), we obtain the velocity distribution:

According to the axisymmetric relation, the average velocity of the slurry can be given by

By introducing (5) into (6), we obtain

The starting pressure gradient of Bingham fluid and the height of flow core area (*h*) relation expression are as follows:

By introducing (8) into (7) and since *λ* is far less than *dp/dx*, ignoring an infinitesimal of higher order, the average slurry velocity can be calculated as

The expression of flow rate *q* of Bingham slurry in the fracture iswhere is the fracture thickness.

The rheological equation of power-law slurry can be expressed as follows:where *c* is the consistency coefficient of power-law slurry and *n* is the rheological index.

By introducing (10) into (1) and integrating, we obtain

It can be known from the hypothesis that ; then, by introducing it into (12), the following equation can be obtained:

According to the axisymmetric relation, the average velocity of the slurry can be given by

By introducing (13) into (14), we obtain

From the generalized Darcy’s Law satisfied by power-law slurry, it can be observed that

The equivalent apparent viscosity of power-law slurry in a single smooth fracture is given by

The average shear rate can be expressed as

Simultaneously, from (17) and (18), it can be obtained that

The expression of the inflow flow *q* of power-law slurry in the fracture is

#### 3. Analysis of the Water-Displacing Mechanism in Fracture Networks

##### 3.1. The Model of Slurry (Water) Diffusion Seepage Resistance in a Single Fracture

In the process of constant pressure grouting in the fracture network, with the increase in the diffusion range of slurry, the flow of slurry inflow in the fracture network will continuously decrease, which is similar to the principle that the current value will continuously decrease when the total resistance of the constant-voltage circuit system increases with time. Based on this, the concept of resistance in the circuit was introduced into the water-displacing diffusion of slurry in the fracture network and the expressions of the resistance of different flow patterns of slurry and water in the diffusion process in a single fracture were analyzed, and then, the expressions of the total resistance of the fracture system at different grouting moments were obtained.

The flow of different patterns of slurry and water in series and parallel fracture structures under constant pressure is taken as the research model, and its schematic diagram is shown in Figure 2. *L*_{1} and *L*_{2} are the horizontal projection length of two series fractures; *p*_{1} and *p* are the pressure difference at their inlet and outlet, respectively; and *b*_{1} and *b*_{2} are their equivalent hydraulic gap widths.

**(a)**

**(b)**

Taking the water-displacing diffusion process of Bingham slurry in the series structural fractures (as shown in Figure 2(a)) as the research object, when the water in the fracture system is completely displaced by Bingham slurry, the flow rate of Bingham slurry in the fractures of *b*_{1} and *b*_{2} at this moment can be expressed as follows:

According to the expression of the total differential pressure *p* in the fracture network and the conservation of slurry flow at the serial fracture joints, it can be seen that

From (21), (22), and (23), the following formula can be obtained:

Similarly, when the slurry-water interface is at the entrance of the fractures in series structure, the water flow rate in the fractures of *b*_{1} and *b*_{2} can be expressed as follows:

Simultaneously, from (25), (26), and (27), it can be obtained that

When the water in the parallel fracture (as shown in Figure 2(b)) is completely displaced by Bingham slurry, the relationship between the total slurry inflow *q* of the parallel fracture system and the flow rates *q*_{1} and *q*_{2} of the two branches of fractures can be obtained as follows:

Simultaneously, from (21), (22), and (30), it can be obtained that

Similarly, at the initial moment of grouting, the total water flow rate in the parallel fracture system can be expressed as follows:

By comparing equations (24), (29), (27), and (30), it can be found that the pressure difference *p* between the inlet and outlet of the fractured structure of different flow patterns is similar to the pressure difference in the circuit system and the flow rate *q* is similar to the current. Based on this, the corresponding resistance expression can be obtained, which is called seepage resistance. For the corresponding seepage resistance of Bingham slurry, power-law slurry, Newtonian slurry, and water diffusion in a single fracture, we call them, respectively, as Bingham slurry seepage resistance, power-law slurry seepage resistance, Newtonian slurry seepage resistance, and aqueous phase seepage resistance and they are, respectively, represented by *R*_{Bg}, *R*_{Mg}, *R*_{Ng}, and . The expression of seepage resistance of different flow patterns of fluids can be obtained from the formula deduced, as shown in Table 1.

##### 3.2. Total Seepage Resistance Model of Water-Displacing Diffusion of Slurry

When the slurry displaces water and diffuses in the fracture network system, the seepage resistance in the whole system will constantly change with the continuous increase in the diffusion range of slurry. Under the constant pressure difference (voltage) of the system, the flow rate (current) of the slurry will constantly change. Therefore, for the fracture network with the series and parallel form, the series and parallel relation of slurry seepage resistance and water phase seepage resistance can be used to represent the total seepage resistance of the system during the process of water-displacing diffusion of slurry, and then, the change of slurry flow (slurry pressure) of the corresponding system at different slurry diffusion positions under the condition of constant pressure (constant flow rate) grouting is obtained, so as to provide certain theoretical reference for the selection of grouting parameters. Earlier, we proposed the concept of total seepage resistance model of water-displacing diffusion of slurry in the fracture network structure. Next, the total seepage resistance method of the system will be used to analyze the variation rule of the inflow flow of the system in the process of water-displacing diffusion of slurry in the herringbone fracture with series and parallel structure under the condition of constant pressure grouting. The selected herringbone fracture structure schematic diagram is shown in Figure 3, and the corresponding parameters are shown in Table 2.

In Table 2, *b*_{1}, *b*_{2}, and *b*_{3} are equivalent hydraulic gap width values of fracture 1, fracture 2, and fracture 3, respectively, and *p*_{1} and *p*_{2} are slurry pressure and hydrostatic pressure, respectively. The changes of total seepage resistance and corresponding slurry flow in the system under the condition of constant pressure grouting and water-displacing diffusion of slurry in the model in Figure 3 were analyzed based on Bingham slurry with different yield stress values *τ*_{0} and power-law slurry with different rheological index values *n* as the research objects.

The viscosity of Bingham slurry is *μ*_{g} = 0.01 Pa·s, and yield stress *τ*_{0} is 1 Pa, 2 Pa, and 3 Pa, respectively; then, the slurry enters the herringbone fracture system from the injection end, the total seepage resistance of the system changes constantly, and the slurry flow rate changes constantly under the condition that the pressure difference of the system is constant. When the water-displacing diffusion of slurry occurs along the horizontal direction of distance *L*_{x} < 10 m, the expression of the total seepage resistance *R* of the system is as shown inwhere *R*_{g1} is the corresponding slurry seepage resistance when the water-displacing diffusion distance of slurry is *L*_{x} and *R*_{w1} is the seepage resistance of the water phase in the horizontal main fracture. *R*_{w2} and *R*_{w3} are the corresponding aqueous phase seepage resistance of the secondary fractures, respectively, and the expressions of the seepage resistance of the slurry and the aqueous phase are shown in Figure 1.

According to the pressure continuity conditions at the slurry-water interface and the conservation conditions of system flow, it can be observed thatwhere △*p* represents the total pressure difference of the fracture system and △*p*_{1} represents the difference between the slurry-water interface pressure and grouting pressure at the slurry diffusion position.

Simultaneously, equations (31) and (32), combined with the expressions of the seepage resistance of Bingham slurry and water phase in Table 1, can be used to obtain the values of the total seepage resistance *R* and the inflow flow *q* of the corresponding system at this position. Similarly, by the same serial-parallel relation, the value of the total seepage resistance *R* and the inflow flow *q* of the system can be obtained when the water-displacing diffusion distance of the slurry is along the horizontal direction *L*_{x} > 10 m.

As shown in Figure 4, with the increase in water-displacing diffusion distance of slurry in the herringbone fracture, the total system seepage resistance increases; when the water-displacing diffusion distance of slurry is small, different seepage resistance yield stress values corresponding to the system total resistance were similar. With the increase in water-displacing diffusion distance of slurry, the total seepage resistance of the slurry system with a higher yield stress value and difference value increases. As shown in Figure 4, when the horizontal diffusion distance of slurry is 13 m, the total seepage resistance *R* (as shown in point A in Figure 4) of the Bingham slurry system whose yield stress *τ*_{0} = 3 Pa is approximately 6 × 10^{10} N·s·m^{−5} and the corresponding slope was a straight line as the growth rate at this point approached infinity.

As shown in Figure 5, when the horizontal diffusion distance of slurry is 13 m, the flow (as shown in point A in Figure 5) of the Bingham slurry system whose yield stress *τ*_{0} = 3 Pa is close to 0 m^{3}/s, the flow of slurry with different yield stress values decreases with the increase in water-displacing diffusion distance, and the reduction rate decreases continuously. When the diffusion distance of slurry is the same, the flow of the slurry with higher yield stress is smaller and the difference increases with the increase in the diffusion distance.

The above analysis shows that, in the case of constant pressure grouting, the total seepage resistance of the system increases with the increase in the slurry diffusion distance, and when the slurry diffusion distance is the same, the larger the yield stress is, the greater the total seepage resistance of the slurry corresponding system is, the smaller the slurry flow rate is, and the more difficult the water-displacing diffusion of slurry is.

#### 4. Verification of the Total Seepage Resistance Model of Water-Displacing Diffusion of Slurry

In order to verify the correctness of the theoretical model of total seepage resistance of the system, the numerical simulation method was used to analyze the variation rules of the slurry-water interface pressure and the flow of slurry inflow at different diffusion positions during the water-displacing diffusion process of the Bingham slurry in the series and parallel fracture structure under the condition of constant pressure grouting, so as to verify the theoretical results mutually.

The Lagrange interface tracking method was used to track the moving distance of the abrupt interface. The process of water-displacing diffusion of slurry was regarded as piston displacement. Putting the calculation area for Ω, slurry area for Ω_{g}, and water area for , the two-phase function *S*(*x*, *t*) expression is

For one-dimensional flows, the Lagrange interface tracking control equation iswhere *u*(s) is the velocity of slurry or water.

According to Darcy’s Law, the expression is obtained as follows:where *µ*(s) is the viscosity value of aqueous phase or slurry phase and *K*(s) is the phase permeability.

The continuity equation of passive single-phase fluid flow in a single fracture is as follows:

For incompressible fluid, (36) can be rewritten as

Therefore, the continuity equation of incompressible two-phase fluid flow in a single fracture can be written as follows:

For the water-displacing diffusion of slurry in the fracture network structure, the overall mass conservation equation of the system can be expressed as

Equations (5)–(1) can be rewritten as

In the process of water-displacing diffusion of Bingham slurry, the corresponding phase permeability *K*(s) is expressed aswhere *τ*_{0} is the yield stress and *b* is the equivalent hydraulic gap width of the fracture.

In the process of water-displacing diffusion of power-law slurry, the corresponding phase permeability *K*(s) can be expressed aswhere *n* is the rheological index and *b* is the equivalent hydraulic gap width of the fracture.

From the previous paragraphs, equations (33), (34), (35), (40), and (41) are the numerical control equations of water-displacing diffusion of Bingham slurry in fracture network. Equations (33), (34), (35), (40), and (41) are the numerical control equations of water-displacing diffusion of power-law slurry in fracture network. By using COMSOL Multiphysics numerical software of low-dimensional PDE module, the control equations are defined directly on all of the fracture networks, and at the same time, in order to prevent the mutation of the phase function *s* between 0 and 1 from the reduction in the convergence of the partial differential equation, the phase function *s* was defined in the form of step function (step(*s*)), so that it could smoothly transition between 0 and 1. The parameter position of the step function was 2.5 × 10^{−5}, the size of the transition zone was 5 × 10^{−5}, and the order of continuous derivative was 2.

The Lagrange interface tracking method is adopted to track the interface. The expression corresponding to the change in the position of the slurry-water interface with time can be expressed as [24]where *ξ(x)*_{n} is the position of the slurry-water interface of the time step; *ξ(x)*_{n+1} is the slurry-water interface location of the next step; *u*_{n} is the interface velocity of slurry and water corresponding to the time step; and *t* is the calculation time step.

To maintain the stability of the calculation, the calculated time step *t* should also satisfy the following relation [25]:where *x* is the shortest fracture length and *l* is the shortest distance between adjacent nodes in the fracture network structure. We use COMSOL Multiphysics numerical software to define the above governing equations and compare the numerical results with the theoretical results in order to verify the correctness of the theoretical model. In the calculation, the form function is Lagrange type and the order of corresponding element is 2. The initial value and boundary conditions of the numerical solution are as follows: at the initial moment, the fracture network is saturated with water, that is, when *t* = 0, *s* = 0, and each inlet of the linear source is Dirichlet pressure boundary conditions. In other words, for the whole model, when *x* = 0, *s* = 1 and *P* = *P*_{1}. The outlet of the model is the boundary of flux source and Dirichlet pressure, that is, when *x* = *L*, *P* = and ∂S/∂*x* = 0. The numerical model and grid division of water-displacing diffusion of slurry in fractures are shown in Figure 6. The length and width of the model are 20 m and 15 m, respectively. Assuming that Bingham slurry viscosity value *μ*_{g} = 0.01 Pa·s, yield stress value of slurry *τ*_{0} = 1 Pa, and the rest of the parameters are as shown in Table 2, the diffusion range of water-displacing diffusion of slurry with the change in the relationship with grouting time is as shown in Figure 7.

**(a)**

**(b)**

**(c)**

**(d)**

As shown in Figure 7, with the increase in grouting time, the diffusion range of water-displacing diffusion of slurry increases continuously, but the growth rate of diffusion range of slurry decreases continuously within the same time interval. The above results show that, in the process of constant pressure grouting, with the increase in grouting time, the advance velocity of slurry-water interface decreases continuously, that is, the inflow flow decreases continuously and the difficulty of water-displacing diffusion of slurry increases continuously. If the above phenomenon is explained by the concept of total seepage resistance of the system, it can be understood that the total seepage resistance of the system keeps increasing with the continuous increase in the diffusion range of water-displacing mechanism. Therefore, under the condition of constant pressure grouting, the value of slurry flow will keep decreasing.

A monitoring line was arranged along the horizontal direction of the model to monitor the position of the slurry-water interface at different times, as well as the numerical values of the slurry-water interface pressure and slurry flow at different times, and the numerical solution was compared with the theoretical solution. The results are shown in Figure 8.

**(a)**

**(b)**

As shown in Figure 8, with the increase in the diffusion range of slurry, the interface pressure of slurry and the value of system slurry flow continuously decrease. It can be seen from the comparison between the theoretical solution and the numerical solution in Figure 8(a) that the difference between theoretical solution and the numerical solution of the slurry-water interface pressure obtained at different diffusion distances is small and so is the difference between the theoretical solution and the numerical solution corresponding to the pressure drop of slurry. Therefore, under the same other conditions, the obtained difference of the slurry flow of the system is also small. Through the above analysis, it can be seen that the pressure of slurry-water interface and the inlet flow value of the system obtained by using the theoretical model of total seepage resistance of the system are basically the same as using the numerical solution; that is, the numerical results and the theoretical results are mutually verified.

#### 5. Conclusions

(1)On the basis of reasonable assumptions, the theoretical model of water-displacing diffusion of slurry of different flow patterns in a single crack under constant pressure was established; on this basis, the concept of seepage resistance of each flow pattern of the slurry when diffusion occurs in a single fracture is put forward and the expression of corresponding seepage resistance is derived.(2)Based on the concept of seepage resistance of each flow pattern of slurry diffusion in a single fracture and combined with the relationship between series and parallel fractures, the concept of the total seepage resistance model of the corresponding system of water-displacing diffusion of slurry in a fracture network structure with series and parallel combination is proposed, meanwhile, taking the herringbone fracture with a series and parallel structure as the research object for example analysis to analyze the change rules of total seepage resistance and slurry flow in the corresponding system with Bingham slurry with different yield stress values at different diffusion locations. The results show the following: for Bingham slurry with different yield stress values*τ*

_{0}of 1 Pa, 2 Pa, and 3 Pa,respectively, when the slurry diffusion distance along the horizontal direction is 13 m, the total seepage resistance

*R*of the corresponding system increases from 3.04 × 10

^{8}N·s·m

^{−5}at the initial time of grouting to 3.78 × 10

^{9}N·s·m

^{−5}, 7.12 × 10

^{9}N·s·m

^{−5}, and 6.1 × 10

^{10}N·s·m

^{−5}, respectively, and the slurry flow

*q*decreases from 4.93×10

^{−4}m

^{3}/s at the initial time of grouting to 3.97 × 10

^{−5}m

^{3}/s, 2.11 × 10

^{−5}m

^{3}/s, and 2.47×10

^{−6 }m

^{3}/s, respectively.(3)The numerical control equations of water-displacing diffusion of different flow patterns of slurry in the fracture network were established, and the numerical analysis of water-displacing diffusion process of Bingham slurry in herringbone fracture structure is carried out. On the basis of mutual verification between theoretical and numerical conclusions, the correctness of the total seepage resistance model of water-displacing diffusion of slurry in the fracture network is verified. Based on the analysis of the diffusion law of different flow patterns in series and parallel fractures, the concept and expression of seepage resistance in the process of slurry flooding and diffusion in fractures are derived, which can be observed as the theory of slurry flooding and diffusion in complex fracture networks. The analysis provides a certain reference, and the related research results can provide theoretical support for the selection of grouting parameters during grouting in fractured aquifers.

#### Data Availability

The data used to support the findings of this study are available from the first author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was supported by Doctoral Research Start-Up Fund Project of Suzhou University (2021BSK014), School Level Key Projects of Suzhou University (2021yzd04), the Green Mine Research Center of Suzhou University (2021XJPT53), and Key Projects of Natural Science Research in Colleges and Universities of Anhui Province (KJ2021A1112).