Erosion Effect on Non-Darcy Hydraulic Characteristics of Limestone and Mudstone Mixture

In order to study the water inrush mechanism of faulted rock mass, a series of erosion seepage experiments were carried out to test the non-Darcy hydraulic properties of a limestone and mudstone mixture. The test results are employed to verify and modify a classic permeability prediction model. Based on the evolution of mass loss rate, the erosion process can be divided into four stages, i.e., particle rearrangement, severe erosion, mild erosion, and stable seepage. During the erosion, the sample height decreases gradually, the water pressure at outlet grows first and then keeps stable, and the porosity and permeability decrease slightly first, then grow gradually and finally keep stable. The hydraulic properties show a more significant variation in the sample within high mudstone particle contents. Through the comparison of test results and the predicted results by the classic Carman-Kozeny model, it is found that the accuracy of the model is greatly affected by lithology. Based on this investigation, a revised model is proposed which introduces a proportional coefficient related to the rock composition, so as to increase the predicted precision in variable rock composition. The forecast accuracy of the revised model is much higher than the classical one in permeability although it decreases with the increase of mudstone content.


Introduction
Fault is a sort of common encountered geological body which is generated by tectonic movement [1]. The interior of a fault is composed of broken rock masses with different lithologies, and the rock stratum on both sides of this geological structure generally has large dislocations [2][3][4]. Faults create discontinuities in the stratum and are usually regarded as water channels that connected the aquifer and underground space. During the mining process, under the water pressure, and mining disturbance, the groundwater is prone to pass through the fault rocks and arrive in the working face, causing water inrush disasters [5][6][7][8]. In this process, a part of fine rock particles inside the fault rocks will migrate with the groundwater. This phenomenon is called the erosion behavior, which is catastrophic in underground construction: First, under the effect of erosion, rock particles and mud would enter the underground space, leading to equipment damage and personnel being buried [9,10]. Second, the erosion effect causes the increase of rock mass porosity and permeability, leading to more groundwater influx [11,12]; Last but not the least, the erosion effect reduces the strength of the fault rock mass, resulting in structural instability and damage [13,14]. Therefore, it is of great significance to study the hydraulic characteristics of fault rock masses during the erosion process to prevent and control fault water inrush disasters [15,16].
In the process of groundwater seepage, the fluid pressure gradient and the fluid velocity are considered to have a certain relationship. Back in the 19 th century, Darcy discovered that the fluid pressure gradient and the fluid velocity obey a linear relationship. Subsequently, this linear seepage law (Darcy's Law) has been widely applied in a variety of studies [17,18]. However, Darcy's law assumes that the fluid regime is dominated by viscous resistance. When the fluid flows at a higher flow rate (e.g., when groundwater flows in a fractured rock mass), the inertial resistance can no longer be ignored, and the water pressure gradient and the flow velocity gradually shows a nonlinear relationship [19,20]. In this case, doubling the water pressure can not obtain a doubled fluid velocity, which results in the inaccuracies of Darcy's law in the calculations and evaluations of water inrush [21,22]. To end this, Darcy's law is usually replaced by Forchheimer equation when the seepage was researched in a high flow rate, which includes a quadratic term of fluid velocity, reflecting the effect of inertial resistance on the fluid regime [23,24]. Many experimental studies [25][26][27][28] have verified the applicability of Forchheimer equation in different conditions, and the Forchheimer coefficient (the coefficient for the quadratic term) for various media forms (e.g., porous media and fractures) is also estimated.
From then on, based on the Forchheimer equation, a lot of scholars have conducted testing or numerical studies on the non-Darcy hydraulic properties of broken rocks. Through fluid testing of sandstone fractures, Zimmerman et al. [28] found that the fluid regime follows the Forchheimer equation when the Reynolds number is higher than 20. And this result is consistent with that of their high-resolution Navier-Stokes simulations. Zeng and Grigg [27] modified the Forchheimer coefficient and adopted this modified parameter as the standard for identifying non-Darcy seepage in porous media. Based on seepage experiments in three different rocks, the critical Forchheimer coefficient is calibrated. Zhou et al. [22] proposed a semiempirical prediction model for tunnel water inflow based on the Forchheimer equation. And the effectiveness of the semiempirical model is verified through numerical calculation results. Chen et al. [21] proposed a framework that combines the estimation of aquifer properties with numerical simulation to more accurately evaluate the flow rate and seepage erosion risk caused by tunnel construction in karst aquifers. This research results show that the predicted flow rate based on the Forchheimer equation agrees better with measured value than that in line with Darcy's law. Ma et al. [29] studied the nonlinear seepage characteristics of granular gangue during the filling process through laboratory and in situ tests and evaluated the reutilization of gangue in protecting overlying aquifers.
However, the above studies did not consider the erosion behavior during seepage. In fact, the erosion effect is a common encountered phenomenon during seepage. Since the 1980s, a series of experiments and numerical studies conducted indepth research on the sand production phenomenon (oil and natural gas flows that drive the movement of fine particles in the rock mass) caused by the oil extraction stage [30,31]. Vardoulakis et al. [32,33] proposed an erosion model of broken rock mass based on the equivalent continuum theory. In this model, the broken rock mass is divided into three phases, which are rock skeleton, water, and fine rock particles. The skeleton particles do not shift during the seepage process, and the fine particles and water flow out at the same speed under the action of water pressure. After entering this century, the erosion problem in the process of underground water inrush has gradually attracted attention [34,35]. In order to study the erosion phenomenon in the process of water inrush, Ma et al. [17,36] developed a set of erosion seepage devices to realize the free migration of fine rock particles in the broken rock mass medium. The porosity and permeability of the rock increase rapidly and then gradually stabilize. Based on the test results, the influence mechanism of erosion on the nonlinear seepage characteristics of fluids is also analyzed. Yao et al. [37] studied the evolution characteristics of the water inrush channel under the influence of particle migration, and based on the established fluid-solid coupling control equation, predicted the water inrush time under different geological conditions. Although there have been many studies simulating the process of water inrush from faults, most of these studies only study the erosion effect of a single kind of rock and do not research the hydraulic characteristics of the mixture by different lithology rocks. In fact, almost all fault rocks are composed of two or more lithologies [38,39]. The interaction between different lithologies and water is very different, so the hydraulic characteristics of mixed lithologies will be different from that of single lithology [40].
In the no. 4 normal fault of the Buliangou mine, the fault rock mass is a mixture of mudstone and limestone (see Figure 1). This fault is in the northeast to southwest direction, with a dip angle of 60-70°and a drop range of 30-50 m. It is located in the middle of the mine field with an extended length of 815 m. In line with the previous research [41], the erosion effect of mudstone after encountering water is significantly stronger than that of limestone. Therefore, studying the influence of different rock components and their mixing ratios inside faults on their hydraulic characteristics is of great significance for preventing water inrush hazards that may occur in faults.
In this paper, limestone and mudstone mixtures in different proportions were selected as experimental materials to test the non-Darcy hydraulic evolution characteristics and porosity evolution of the samples under erosion. Based on the experimental results, a method for predicting the permeability of a limestone and mudstone mixture is proposed. The research results can provide a basis for the prevention and treatment of fault water inrush disasters.

Testing
Material. The test materials are limestone and mudstone sampled from Buliangou Mine. The dry density of limestone is 2.42-2.71 g/cm 3 , and the dry density of mudstone is 2.35-2.61 g/cm 3 . Since the difference in dry density between the two is small, we approximate both as 2.5 g/cm 3 . Rocks are crushed into 0-20 mm crushed particles first and then divide the crushed particles into five groups by diameter, namely, Group A (0-2 mm), Group B (2-5 mm), Group C (5-10 mm), Group D (10-15 mm), and Group E (15-20 mm). Before the test, these particles are mixed according to the continuous gradation ratio (Equation (1)) and take into account the different ratios of the two rocks (denoted by the mass content of mudstone).
where α is the group number of crushed particles, d α is the largest particle diameter in the Group α, d max is the largest diameter of all groups, P α is the mass proportion of the groups which size is smaller than d α , and β is gradation parameter denoting the gradation properties of mixed rock samples, which is set as 0.5 in this study. The sample mass is set as 4300 g and the mixing scheme is shown in Table 1, and the continuous grading curve of the sample is shown in Figure 2.

Testing System.
In order to carry out erosion tests on the limestone mudstone mixture, the following test system is designed to carry out experiments on the limestone mudstone mixture samples. As shown in Figure 3, the test system is composed of the fluid supply system, the data acquisition system, the particle collection system, and the seepage erosion unit.

Fluid Supply System.
It is mainly composed of the oil pump, the water pump, the double-acting hydraulic cylinder, the connected pipelines and valves. The double-acting hydraulic cylinder consists of two cavities, and the volume of the cavities is controlled by the piston. Before the test, valve 1 is closed and valve 2 is opened, the water is injected into the lower cavity of the hydraulic cylinder through the water pump. After the test begins, valve 2 is closed and valve 1 is opened. When the oil pump is started, hydraulic oil is injected into the upper cavity, driving the piston to compress downward, so that provides stable water pressure to the erosion unit.

Erosion Unit.
It comprises the piston, the cylinder, the upper and lower seepage plate, the conical bottom plate, the base and other parts. At the beginning of the test, the fluid supply device provided water pressure to the piston inlet.

Geofluids
Small holes of 2 mm in diameter are arranged on the upper seepage plate to disperse the water flow. And small holes with a diameter of 10 mm are arranged on the lower seepage plate to facilitate the flow of broken rock particles with a diameter of less than 10 mm. The conical bottom plate can collect the fine particles from erosion outflow and then lead them to the particle collection system.

Data Acquisition
System. It is composed of a linear variable differential transformer (LVDT), a water pressure gauge, a flow gauge, a recorder and a computer. LVDT is arranged at piston and floor to monitor the collapsed height of rock mass during erosion. Two water pressure gauges are arranged at the inlet and outlet of the erosion unit to monitor the water pressure at the water inlet and outlet, respectively. The flow meter is arranged at the entrance of the erosion unit to measure the fluid flow during the test. Each instrument is connected with the recorder and computer to record and output the measured data in real time.

Particle Collection
System. It is composed of a particle collector, a filter screen, an electronic scale, and a water tank.
This system can collect and weigh the fine particles flowing from the rock sample in real time. The reading of the electronic scale is zero at the beginning of the test, so the reading of the electronic scale records the mass change of the collection tank during the test. The fine particles will discharge the same volume of water when they flow into the collection tank. Because the density of fine particles is greater than the density of water, based on the mass change of the collection tank, we can calculate the volume of the fine particles by the density difference between the fine particles and the water.   Figure 3: Test system. Note: (a) fluid supply system, (b) erosion unit, (c) data acquisition system, and (d) particle collection system.
Test preparation: First, the seepage pipeline is connected, and the data acquisition system is adjusted to make it work normally. After this, the double-acting hydraulic cylinder is filled with water. Then, the broken sample is put into the erosion unit, and the water is injected into the pipeline to exclude air in the pipeline and saturate the sample. Subsequently, the piston is placed and the initial height of the sample is recorded.
Procedure 1 Water pressure loading: the oil pump is started, and water is injected into the erosion unit, the data from the electronic scales in the particle collection system and the data measured by the acquisition system and are monitored.

Procedure 2
Posttest treatment: oil pump is turned off first, and then stop the water flow. Finally, the sample is taken out from the erosion unit.

Procedure 3 4 Geofluids
where ρ r is rock dry density, ρ w is water density, and the accumulated mass loss is The mass loss at time interval i, ΔM i , can be calculated by And the mass loss rate _ M i , which denotes the mass loss in unit time, can be acquired by the following equation: 2.4.2. Sample Height. As shown in Figure 4, the dimensions of each part of the test instrument are piston height h 1 = 100 mm, cylinder upper edge thickness h 2 = 15 mm, permeability plate thickness h 3 = 10 mm, and cylinder height h 4 = 350 mm. According to the geometric relationship in Figure 4, we can get where h s is the sample height and h is the measured height through LVDT. Substituting the dimension data of each part, the sample height at time t i can be obtained by the measured value of LVDT: where M 0 is the initial sample mass and M 0 = 4300 g and r is the radius of the cylinder and r = 50 mm.
In line with Forchheimer's equation [27,36], there is a well-known relationship between water pressure gradient and flow velocity: where μ is the viscosity of fluid, k is the permeability, and β is the non-Darcy factor. Assuming that the water pressure is linearly distributed in the sample, the water pressure gradient can be expressed as where p A and p B are the water pressure at fluid inlet and outlet, respectively. According to the test principle, p A is a constant and p A = 1:5 MPa, and p B is a variable changing with time, which can be measured by the pressure gauge. Combining Equations (10) and (11), there is For each moment t i , Equation (12) can be written as According to previous studies [29,42], the permeability at time interval i could be approximated to the average value of that at time interval i and i + 1, so that the permeability 5 Geofluids can be calculated by Equation (14).  Figure 5. According to the variation trend of the mass loss rate, the entire erosion process can be divided into four stages, i.e., particle rearrangement, severe erosion, mild erosion, and steady seepage.

Test Results and Discussions
(1) Particle rearrangement: in the initial stage of the test, under the action of water flow, fine particles are dragged into the pores of large particles to form a more stable structure. At this stage, fewer particles are discharged. The amount of quality loss is small In the comparison of samples under different mudstone contents, it can be found that in samples with higher mudstone content, the particle rearrangement stage is shorter and the peak mass loss rate is more significant, and it enters the stable seepage stage soon after the severe erosion, which means there is a shorter mild erosion stage. This result manifests that the erosion effect of mudstone is stronger than that of limestone. Under the action of water flow, mudstone particles are more prone to migrate.

Sample
Height. The evolution of the sample height during the seepage is shown in Figure 6. The height of all the samples has decreased under the action of erosion. The sample heights show a sharp decline in the early stage of the erosion test, which is because the rock particles are redistributed under the effect of water flow, and lots of fine particles are filled in the pores between the large particles, contributing to the rapid decreases in sample height. Then, the fine particles move and under the effect of the water flow, and finally flowed out from the bottom, so that the sample structure collapsed and the height declines. In these stages (severe erosion and mild erosion), the drop rate of the sample height gradually slows down. Finally, and sample height remained stable, which indicates that the erosion effect has terminated.
For the height before erosion, it can be observed that the heights of the samples with different rock composition ratios are the same. This shows that lithology has almost no effect on the natural accumulation height of the samples. As the erosion process progresses, it is found that the more mudstone content is, the height of the sample drops faster. This shows that under the action of water flow, the higher the mudstone content, the more significant rearrangement and erosion in the sample.   Figure 7 illustrates the porosity evolution of the sample. It can be found that in the first stage of erosion, i.e., the particle rearrangement stage, the porosity of the sample keeps stable or decreased slightly. This tiny decline in porosity is because there are a lot of small rock particles filled into the pore of large rock particles under the effect of water flow. In the second and third stages, a substantial increase in the porosity of samples can be observed. This corresponds to the phenomenon discussed earlier in which particles flow out of the sample in large quantities. Finally, the porosity of each sample has fallen into a stagnant stage, and the pore structure of the sample stays in a stable state, which indicates that most of the flowing particles have left the sample, and the migration of fine particles has not been observed. Comparing the porosity evolution curves of different samples, it is found that in samples with higher mudstone content the porosity has increased significantly, and the final porosity of the rock mass is also greater. This indicates that the erosion effect is stronger in these samples.

Water
Pressure at the Outlet. The evolution curves of water pressure at the outlet are drawn in Figure 8. It is noted that in previous studies, the water pressure at the outlet of the sample is often regarded as connected to the atmospheric pressure and assumed to be 0. The results of this test show that the water pressure at the outlet of the sample is much greater than 0. Such an effect of water pressure at the outlet can not be ignored, especially when the water pressure at the inlet is small, improperly neglecting the hydraulic power at the outlet will harm the precision of the test results. In addition, the water pressure at the outlet is not a stable value: with the progress of the erosion, it can be found that all the samples have the phenomenon of rising water pressure at the outlet. This is because, under the erosion effect, the fine particles in the sample gradually flow out, and the resistance to the fluid decreases, resulting in a drop in the water pressure difference between the two ends of the sample.
Comparing different samples, it can be found that although the outlet water pressure of each sample is almost the same at the initial moment; in the samples with more mudstone content, the outlet water pressure has increased on a larger scale. This is due to the more significant erosion effects that appear in these samples.  7 Geofluids 3.1.5. Permeability. Figure 9 describes the variation in the permeability of the sample. It can be observed that the permeability of the sample has a similar variation trend comparing with the porosity of the sample, that is, the permeability of the sample decreases slightly with the progress of erosion, then gradually increases and finally trend plateau. The highest permeability can be observed in samples with higher mudstone content, and the time required for the increase in permeability in these samples is shorter than that in other samples, which means more severe permeability changes and a higher risk of water inrush.

Permeability Prediction Model Based on Porosity.
In the process of underground water inrush, an important issue is to determine the permeability change of the rock mass. Unfortunately, due to the influence of the geological environment, direct measurement of permeability of fault broken rocks is a challenging task [43,44]. Therefore, many previous studies have used other properties of the rock mass to indirectly predict its permeability [45,46]. For example, in the well-known Carman-Kozeny model [47,48], permeability is described as a cubic function of porosity. Based on this model, the permeability can be calculated when the porosity of the rock mass is obtained. On the basis of non-Darcy hydraulic characteristics of the limestone and mudstone mixture, the Carman-Kozeny model is employed to predict the permeability in this research, and the predicted results are compared with the testing results to verify the reliability of the equation.
The Carman-Kozeny model is as follows: where k R is constant, and in the initial time, there is If k 0 and ϕ 0 are known, combining Equations (15) and (16), we get To evaluate the model's accuracy, the average percentage error (APE) is induced to describe the difference between testing and predicted values in each data point: where k p i and k t i are the testing and predicted values of the plot i (i = 1, 2 ⋯ n). And the mean average percentage error (MAPE) is employed to describe the mean difference between testing and predicted values in a certain sample. 8 Geofluids Based on Equation (17), the predicted curves are obtained as shown in Figure 10. It is obvious that the trends of the predicted and measured results are semblable; that is, as the porosity increases, the permeability of the rock mass gradually increases. As the porosity increases, the prediction curve gradually deviates from the measured value point, which indicates that the prediction performance becomes worse when the porosity is large. By observing the APE distribution, the APE shows a high dispersion under different porosity conditions. According to the MAPE value of the model, it is indicated that the accuracy of the model is poor, the maximum MAPE value reached 32.25% (sample 5), which means that the model is not suitable for predicting permeability.
Comparing the predicted results in different samples (see Figure 10(b)), it is seen that in the sample with a 1 : 1 ratio of mudstone to sandstone, the prediction accuracy of the model is the highest, and the MAPE value is only 5.2%, while for samples with other mixed ratios, the MAPE is much higher. This phenomenon shows that the accuracy of the model is greatly affected by the lithology of the sample, and there is a certain relationship between the predictive performance of the model and the composition of the rock, which is worthy of further study.

Modified Prediction Model Based on Porosity and
Mudstone Content. Based on the analysis in Section 3.2.1, to further study the relationship between the predictive performance of the model and the rock mixed ratios of the sample, the fitting lines between testing values and predicted values in different samples are drawn and shown in Figure 11. In Figure 11, the abscissa represents the predicted value, and the ordinate represents the testing value. When the predicted value is the same as the testing value, their fitted straight line is the equality line, and there is no error in the model at this time. When the fitted line deviates from the equality line, it means that the model has an error in this prediction, and the closer the fitting slope is to 1, the better the model's predictive effect. From this figure, it is clear that for samples with high mudstone content, the fitting line's slope exceeds 1, which means the predicted value is lower than the testing value and the model underestimates the testing results. On the contrary, in the samples containing more  10 Geofluids limestone, the fitting line's slope is lower than 1; that is, the model overestimated the experimental value. In addition, according to the R 2 value of fitting curves, it can be found that the fitting accuracy of each sample is high, indicating that the predicted value of the model and the experimental value present an obvious linear relationship. Due to the linear relationship between the predicted value and the testing value, this research tries to add a scale factor to modify the Carman-Kozeny model. By introducing the variable of mudstone content, the modified model could consider the effect of variation in the proportions of different rock components on permeability prediction. Figure 12 shows the relationship between the slope of the fitting line ξ and the mudstone content c. It can be found that the slope of the fitted line and the mudstone content present an exponential distribution, and the fitting equation is given below: Based on the above discussion, we propose the following modified formula to predict the influence between permeability and porosity of limestone mudstone mixture: The above prediction equation fully considers the influence of mudstone content. According to the modified model (Equation (21)), we have predicted the permeability of the limestone mudstone mixture. In order to describe the accuracy of the predicted surface, the MAPE of the total sample is introduced and can be obtained by Equation (22): where j is the number of samples. The prediction results of the modified model are shown in Figure 13. In the comparison of the predicted value and the measured value, it can be found that the testing value point is very close to the predicted surface. Through the predicted surface, it is observed that as the porosity and the mudstone content increase, the permeability of the rock mass increases. As shown in Figure 13(b), although the dispersion of APE is still high, according to MAPE results, the accuracy of the modified model has been significantly improved. The MAPE of all samples is within 6%, and the MAPE of the total sample is 4.46%. Comparing the results of different rock samples, it can be found that as the mudstone content increases, the accuracy of the model gradually decreases (except when comparing samples 3 and 4). This may be due to the formation of more argillaceous material in the sample with higher mudstone content, which blocked the seepage channel, resulting in much higher uncertainty during the erosion process.

Conclusions
In this article, a series of erosion seepage experiments are carried out to investigate the non-Darcy hydraulic characteristics (e.g. mass loss rate, sample height, porosity, water pressure at the outlet, and permeability) of the limestone and mudstone mixture were studied. The test results are adopted to evaluate the accuracy of a classic permeability prediction model and then put forward a modified model considering the lithology of the medium. The main conclusions are as follows.
The erosion process can be divided into four stages: particle rearrangement, severe erosion, mild erosion, and stable seepage. The higher the mudstone content in the sample, the stronger the erosion effect in the sample. Under the action of erosion, the height of the sample decreases rapidly and then remains stable. The higher the mudstone content, the faster the height of the sample decreases under the effect of erosion, and the greater the decline. The porosity of the sample increases first and then remains stable. The higher the mudstone content, the more significant the increase in porosity. The water pressure at the outlet gradually increases with the erosion effect. The higher the mudstone content, the more obvious the increase. Permeability and porosity growth have a similar trend; that is, it increases first and then stabilizes under erosion. A higher porosity growth is observed in samples with high mudstone content.
Through the same test results and Carman-Kozeny model prediction results, it is found that the prediction performance of the classic model is unaccepted, the maximum MAPE is more than 32%, and the accuracy of the model is greatly affected by lithology. Further analysis of the model results shows that there is an obvious linear relationship between the predicted value and the test value. Based on this phenomenon, a proportional coefficient related to the rock composition is introduced to improve the Carman-Kozeny model. By comparing the test and calculation results, the prediction accuracy of the modified model is much higher than that of the classic one, and the MAPE values are within 6%. In addition, the predicted accuracy decreases with the increase of mudstone content. This phenomenon may be related to the clogging effect of muddy substances on the water conducting pathway.

Data Availability
The experimental data used to support the findings of this study are included in the article.

Conflicts of Interest
The authors declare no competing financial interest.