#### Abstract

Rock fractures always influence the hydrological properties of a rock mass. To investigate the seepage characteristics of a rock mass with partly filled fractures, a mathematical model is established. In this model, the clear fluid in fractures is governed by the Navier-Stokes equation, and the fluid both in the porous medium and rock matrix are subjected to the Brinkman-Extended Darcy equation. The analytic solution of an equivalent permeability coefficient for a rock mass with partly filled fractures is solved, and it could be reduced to some special known results. Comparisons with experimental data show good agreement, thus verifying the validity of the present computations.

#### 1. Introduction

Natural rock masses are composed of intact rock blocks with numerous fractures. Rock fractures can greatly influence the hydrological properties and control the mechanical behaviour of a rock mass. The hydraulic and mechanical properties of rock fractures are of considerable interest in several areas of rock engineering, such as near-field modelling of hydraulic processes around radioactive waste repositories [1, 2], underground tunnels [3, 4] and storage caverns [5, 6], rock slope stability [7, 8], petroleum and geothermal reservoirs [9], and solute transport [10]. Single fractures, the elementary unit of a fractured rock mass, have become a research focus. Lomize [11] presented the open Cubic Law, which stems from the assumption that a fracture consists of the region bounded by two smooth, parallel plates. Louis [12] carried out penetration tests on rock with a single fracture and proved the validity of the Cubic Law when assuming laminar fluid flow. The classical view of a rock fracture is a pair of smooth, parallel plates.

At present, detailed studies on the permeability of a rock mass with an unfilled fracture have been conducted, and they focus on the influences of stress [13–15], fracture aperture [16–18], and fracture roughness [19–21]. In addition, Yeo et al. [22] found that with the increase of shear displacement, a fracture became heterogeneous, anisotropic, and more permeable in the direction perpendicular to the shear displacement rather than in the direction parallel to the displacement. Arbogast and Brunson [23] presented two computational models and determined that vug connectivity was the most critical variable in predicting macroscopic properties among layered vugs, meandering vug channels, constricted vug channels, and disconnected vugs. Yang et al. [24] conducted indoor tests of fluid flowing through a single rough fracture to investigate seepage mechanisms and the effect of the structure roughness to the fluid.

Actually, natural fractures always contain infilling material, which affects the overall permeability of the fractured rock mass. Takemura et al. [25] conducted a series of experiments and revealed that the permeability of a rock mass with filled fractures is related to fracture aperture and porosity of the infilling material, and they proposed a semiempirical permeability formula. Deng and Martinez [26] solved the velocity distribution of two-dimensional fluid flow in fractured and porous media using the stress-jumped boundary condition proposed by Ochoa-Tapia and Whitaker [27]. Li et al. [28] set up a Darcy-Stokes model and derived the generalized Cubic Law and the expression of equivalent permeability of a rock mass with an infilled single fracture in tensor form, considering a permeable rock mass.

In a natural geological body, there are not only completely filled fractures and empty fractures but also partly filled fractures. However, there are few studies that have investigated seepage in a rock mass with partly filled fractures, despite much research on the seepage of fractured rock masses having been conducted. In this paper, a nonlinear mathematical model is proposed. In this model, the motion of fluid in a clear fracture is governed by the Navier-Stokes equation, and seepage flow in infilling material and rock matrix is subjected to the Brinkman-Extended Darcy equation. Through analysis of the results, the flow velocities of fluid in open and filled fractures and the rock matrix were derived, as well as an analytical solution for the equivalent permeability coefficient () for a rock mass with partly filled fractures. The analytical solution could be reduced to some special known results. Finally, simulation tests were conducted by the authors for flow in the modelled geometry using a specially designed device in which porous stone was used to simulate the rock mass and the fracture was filled with sand to different heights. The modelling results show good agreement with experimental data previously obtained by the authors for flow in.

#### 2. Mathematical Model

In reality, a rock mass is not intact, as substantial fissures or fractures are present that were generated over a long period of geological history. Some fractures are filled with debris and sediment, some are open, and others are partly filled. Clearly, the permeability of a rock mass varies with the degree of fracture fill, which may also exert a significant impact on the safety of underground tunnels and the stability of rock slopes. To get a better understanding of the permeability characteristics of a rock mass with partly filled fractures, a periodical seepage analytical model was designed using a Cartesian coordinate system (), as shown in Figure 1, and a theoretical analysis was carried out based on the following assumptions: (a)Fluid flow is infinite along the direction only(b)The fluid is a fully developed Newtonian fluid(c)The fluid is incompressible, and the volume force along the direction is ignored(d)The fluid is viscous with laminar flow(e)The fracture/crack extends infinitely in the direction(f)Velocity and shear stress are continuous at the interface between different media(g)The fluid in fractures, fillings, and the rock matrix satisfies the continuity equation(h)Seepage through crack fillings and the rock matrix satisfies the Brinkman-Extended Darcy equationwhere is the seepage velocity in soil (); is the fluid pressure (); is the kinematic viscosity of water (); is the density of water (); is the mass force (); is the porosity of soil (dimensionless); is the permeability of soil (); is the gravitational acceleration constant; is the migration derivative; and is the Hamilton operator.

The clear fluid in fractures can be described by the Navier-Stokes equation: where is the flow velocity of clear fluid (), and other parameters are the same as above.

##### 2.1. Interface Conditions

The description of the interface conditions between a clear fluid and a porous medium has received considerable attention. In summary, three primary categories of conditions at the interface were found in the literature; these can be classified into two main conditions: slip and no slip [27, 29, 30].

Beavers and Joseph [29] presented the BJS boundary condition, where velocity is not continuous at the interface with a different medium and velocity meets the following condition: where is the flow velocity of fluid in the fracture () and is the Darcy flow velocity in a porous material (). The value of is in the range of 0.1~4.

Neale and Nader [30] have questioned the validity of the above condition. They proposed that both the velocity and shear stress must be continuous at an interface, using the following equations: where denotes the effective dynamic viscosity () and .

Ochoa-Tapia and Whitaker [27] presented the stress-jump condition, which states that the velocity is continuous but the shear stress is not. The stress-jump condition can be expressed as follows: where is an empirical constant that can be positive or negative.

Furthermore, we find that the interface description of each category is not a closed system of equations. The success of each model is dependent on the accuracy of a free integration constant, i.e., , , or , which must be independently determined. As and are empirical coefficients and there are no general methods to derive them, the boundary condition proposed by Neale and Nader was adapted in the present study.

##### 2.2. Flow Velocity Field

Based on the theoretical model and associated assumptions, a period unit was selected for analysis, as shown in Figure 1. is the length of fracture, is the thickness of one-period rock, and is the fracture aperture. The flow velocities through crack fillings and the rock matrix are and , respectively, while denotes the velocity of a clear fluid. The flow velocities , , and are discussed separately.

###### 2.2.1. Flow Velocity

The clear fluid in the fracture is governed by the continuity equation: where , , and are velocity components in the directions , , and in units of .

The motion equation satisfies the Navier-Stokes equation:

As the fluid velocities are 0 in the and directions, namely, , therefore, . By substituting into equation (6), is obtained. As is constant in the direction, namely, , therefore, . By substituting these equations into equation (7), equation (8) can be obtained as follows:

According to the condition that the gradation variation of along the direction is a constant, the expression of can be written as , where denotes the water pressure difference along the direction () and is the seepage distance along the direction ().

Together with the two conditions of and , the second partial derivative of with respect to is directly transformed to the second derivative. Then, substituting into equation (8) leads to

By substituting into equation (9), the following result is obtained:

The solution for equation (10) gives an expression for : where and are undetermined parameters.

###### 2.2.2. Flow Velocity

The fluid in fillings satisfies the continuity equation: where , , and are velocity components in the , , and directions in units of .

The motion equation satisfies the Brinkman-Extended Darcy equation: where and are the porosity and permeability of the crack-filling material, respectively.

Similar to the above analysis for equation (7), equation (13) can be simplified and rewritten as follows:

By solving equation (14), an expression for is obtained: where and are undetermined parameters.

###### 2.2.3. Flow Velocity

Fluid in the rock matrix satisfies the continuity equation: where , , and are components of velocity in the , , and directions in units of .

The motion equation satisfies the Brinkman-Extended Darcy equation: where is the porosity of the rock matrix and is the permeability of the rock matrix in units of .

Similar to the above analysis for equation (7), equation (17) can be simplified and rewritten as follows:

By solving equation (18), an expression for is obtained: where and are undetermined parameters.

##### 2.3. Boundary Conditions

According to the mathematical model, the following boundary conditions are satisfied: (a)At the interface between different media, satisfies the condition that flow velocity and shear stress are continuous(b)At the interface between different media, satisfies the condition that flow velocity and shear stress are continuous(c)At the interface between different media, satisfies the condition that flow velocity and shear stress are continuous

These boundary conditions could be expressed as follows: where denotes the thickness of one-period rock mass (); is the width of the fracture (); and is the ratio of the thickness of the crack-filling material to that of the fracture (dimensionless).

Substitute these boundary conditions into equations (11), (15), and (19), we get

Parameters , , , , , and can be obtained from equation (21). The detailed derivations will be described in the appendix.

##### 2.4. Equivalent Permeability Coefficient

The seepage field in a layered rock mass is a three-dimensional anisotropic field. In the horizontal plane, the flow velocities are the same. The permeability in the horizontal plane can be solved using the theory described above.

###### 2.4.1. For along Direction

As shown in Figure 1, the average flow velocity () in the direction can be expressed as follows:

Darcy’s law states that where is the permeability coefficient in the direction in units of and is the hydraulic gradient of a rock mass in the direction (dimensionless).

Additionally, we have where denotes the water head difference (). where is the unit weight of water (). where is the permeability of a rock mass along the direction in units of .

The equivalent permeability of a rock mass can be obtained by solving equations (23), (24), (25), and (26):

The analytic solution for the equivalent permeability coefficient () of a rock mass with partially filled fractures is derived as follows:

###### 2.4.2. For along Direction

As shown in Figure 2, water discharge from each layer should be equal along the direction, which can be described as follows: which can be expressed as follows: where , , and are the permeability coefficients of the water layer, infill layer, and rock layer, respectively, and , , and are the hydraulic slopes of the water layer, infill layer, and rock layer, respectively.

As the flow velocity in each layer is equal, it can be described as follows:

And because the total water head loss should be equal to the sum of the water head loss in each layer, we have where , , and are the thicknesses of water layer, infill layer, and rock layer, respectively, and is the total water head loss.

According to equations (30), (31), and (32), the permeability coefficient of a periodically partially filled rock mass along the direction is derived:

As in equation (33), the expression of can be rewritten as follows:

And the permeability of a periodically partially filled rock mass along the direction is as follows: where and .

Equation (35) has also been reported previously [31].

#### 3. Validity

##### 3.1. Experiment

###### 3.1.1. Testing Device

To verify the above analysis, a penetration device was designed; it is shown in Figure 3. The device consists of three parts: a water inlet device, a seepage device, and a water collection device. The seepage device is an organic box with dimensions of and is shown in Figure 3.

In this study, a porous stone produced by Hongzhou Experimental Instrument Co. Ltd. was used as a substitute for the rock matrix. For the sampling procedure, waterproof glue was applied to the inner surface of the organic box to avoid permeation, and the porous stone was fixed by applying waterproof glue (, along the seepage direction). Then, a foam board with a fixed thickness was put in the middle of the organic glass box, and the fracture was then filled with sand with a specified density and particle size of 0.25~2 mm. After the fracture was completely filled, the foam board was slowly pulled out and two porous discs were fixed at both the left and right sides of the porous stone. Finally, the penetration device was connected as shown in Figure 3.

The main steps during the test procedure are as follows: (1)Turning on the switch of the water inlet device and water will flow into the seepage device and slowly get the testing sample saturated(2)Clamping the discharging tube by a water stop clamp and adjusting the altitude of the water inlet device to get the designated water head difference, while the seepage device was full of water(3)Removing the water stop clamp and beginning to record the time for the specified amount of seepage after the flow became steady

###### 3.1.2. Determination of Permeability and Porosity of the Sand Fill

Constant head permeability tests were carried out as per GB/T 50123-1999 [32]. The testing device is shown in Figure 4. The inner diameter of the sampling pipe was 10 cm, and the distance between piezometer tube I and III was 20 cm.

The permeability coefficient can be determined using
where is the seepage volume per unit time (). is the height of sand between the centres of two pressure holes (10 cm), is the area of the specimen section (78.54 cm^{2}), and is the average water head difference, where .

The dry density of sand () was maintained at 16.14 kg/cm^{3}. The constant head permeability test results are shown in Table 1.

According to ( is taken to be at 20°C), the permeability of the fill material can be obtained, i.e., .

The porosity of the sand fill can be obtained as follows [33]:
where is the porosity of sand and is the dry bulk density of sand (). , is the bulk density of water (9.8 kN/m^{3}), is the specific gravity of the sand fill, and .

By substituting these parameters into equation (37), the porosity of the sand fill can be derived, i.e., .

###### 3.1.3. Determination of the Permeability Properties of a Rock Mass with a Partly Filled Fracture

The experiment was carried out using the procedure presented above. In the experiment, three kinds of foam were used, with thicknesses of 1 mm, 2 mm, and 3 mm. The results for these tests are shown in Tables 2, 3, and 4, respectively.

By substituting the parameters of the sand fill (, ), porous stone (, ), and the size of sample into equation (28), the results for seepage using the same geometry can be derived. The calculated values for the permeability coefficient with different degrees of fill are shown in Table 5. Table 5 also lists the test results and the computation results. From Table 5, it is obvious that the test values are slightly larger than the calculated values, and the errors between the two values approximately equal 8.0%. However, the discrepancies between the two results are small and are within an order of magnitude. After reviewing the test process, the discrepancy may be attributed to the following: (1) the assumption is ideal, which cannot be fully realized in experiments, i.e., the infinite extension of fractures and the homogeneity of porous media and boundary conditions; (2) density differences between a real sand fill and the measurement before; and (3) winnowing away of some sand particles during testing.

##### 3.2. Deduction of Some Special Known Results

###### 3.2.1. Unfilled Fractured Rock Mass Model

If , then fractures in a rock mass are not filled. Hence, the seepage model reduces to an unfilled fractured rock mass model, and then equation (28) will be reduced to

This result agrees with a well-known result reported earlier by Arbogast and Lehr [34] for the equivalent permeability of a rock mass with unfilled fractures, i.e., , and the difference may come from the choice of different boundary types. In that work, Arbogast chose the BJS boundary condition; this means that velocity is not continuous, and the velocity at the interface of a different medium meets the condition of .

###### 3.2.2. Wholly Filled Fractured Rock Mass Model

If , then the model will be reduced to the wholly filled fractured rock mass model, and equation (28) will be reduced to

###### 3.2.3. The Cubic Law Model

If , and if rock is regarded as a waterproof material, then equation (28) can be rewritten as follows:

This result is in good agreement with the well-known result reported by Lomize [11], which is for the equivalent permeability of a rock mass when considering rock as waterproof. The result is known as the Cubic Law.

##### 3.3. Application

Here, we use the Yuelongmen tunnel of the Chengdu-Lanzhou railway as an applied example of this method. This tunnel is located in the Longmenshan fault zone. The total length of the tunnel is 19.981 km. The main rock types at the tunnel site are limestone, dolomite, and mudstone with occasional granulite, sandstone, and quartzite.

As a result of crossing the Longmenshan fault zone, especially after undergoing the Wenchuan earthquake with Ms 8.0 on May 12th, 2008, the rock mass became relatively fractured, and joints and fissures developed widely. According to a field investigation, the seepage volume at a location 9.4 km from the entrance of the Yuelongmen tunnel amounted to 1682 m^{3}/d, and problems such as serious water gushing, water in-rush, and mud outbursts occasionally occurred when construction was occurring on these fault sections and fissure development zones, which presented a serious obstacle to construction.

To investigate seepage while drilling a tunnel in a zone of fissured rock, a simulation device was designed, as shown in Figure 5, which has been issued by Fu et al. [35, 36]. The size of the model is shown in Figure 6. In the experiment, layered bricks were utilized as a stratified rock mass and sand was used as fill. Water channels on both sides of the model were built as a water source, several water inlet pipes were fixed at the internal side of the channels (see Figure 6) and the water table in the channel was kept at 115 cm while the test was conducted. The tunnel was synthesized using a plastic pipe with a diameter of 10 cm, and massive holes were uniformly distributed along the surface of the pipe as shown in Figure 7. A cylindrical sealing plug was used to plug the pipe before the test. The seepage volume was collected and weighed at 60 s, 120 s, 180 s, 300 s, 480 s, 690 s, 900 s, 1200 s, 1500 s, 1800 s, 2400 s, and 3000 s, and the result is shown in Table 6. The density of water was considered to be 10^{3} kg/m^{3}.

The results from the test and computation based on this method are shown in Figure 8. Figure 8 illustrates that the test values are slightly larger than the calculated values and that the accumulative seepage volume increased as time increased. From Figure 8, it is obvious that the difference between the two values rises with increasing time and that the gap is relatively small before 300 s, is markedly enlarged after 1800 s, and is at the maximum at 3000 s.

#### 4. Conclusions

Rock fractures can greatly influence the hydrological properties of and control the mechanical behaviour of the rock mass. To investigate the seepage characteristics of a rock mass with partly filled fractures, a modified model from the Cubic Law is established, which assumes that a clear fluid in a fracture is governed by the Navier-Stokes equation and that the fluid both inside the porous medium and the rock matrix are subject to the Brinkman-Extended Darcy equation. Through analysis, the analytic solution for the equivalent permeability coefficient () of a rock mass with partially filled fractures is solved, and it can be used to deduce some special known results. To verify the validity of the model, simulation tests were conducted by the authors for flow using the same geometry with a specially designed device.

The results show that the equivalent permeability coefficient () of a rock mass with partially filled fractures is very sensitive to the ratio of fill (), where a smaller value of gives a larger value for the equivalent permeability coefficient. Moreover, comparisons with experimental data previously obtained by the authors for flow using the same geometry are in good agreement. The results determined from the calculated solution for the equivalent permeability coefficient are slightly larger than that from the simulations.

In conclusion, the analytical solution of the equivalent permeability coefficient for a rock mass with partially filled fractures is more comprehensive and reasonable, and it can not only be adopted for calculating the permeability coefficient of a rock mass with filled fractures and unfilled fractures but can also be used for partially filled fractures.

#### Appendix

The flow velocities in infilling sand and in the rock matrix are and , respectively, while denotes the velocity of a clear fluid (see Figure 1).

According to the analysis in Section 2.2, the flow velocities are as follows:

The following three boundary conditions are given: (1)At , according to and , we obtain (2)At , according to and , we obtain Solving equations (A.4) and (A.5) leads to (3)At , according to and , we obtain

Solving equations (A.8) and (A.9), we obtain

Combining equations (A.8) and (A.9), we obtain

As and , so and . According to and , we get . Similarly, we also obtain . Then, substituting them into equations (A.2) and (A.3) leads to

Thus, the average flow velocity () in the direction can be expressed as where

Substituting equations (A.4) and (A.8) into equation (A.15), we obtain where

As and , substituting , , , and into equation (A.19), we obtain

Substituting equation (A.21) into , we obtain

Substituting equation (A.22) into , we obtain

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

#### Authors’ Contributions

The contributions of Prof. Fu included conceptualization; funding acquisition; project administration; supervision; and writing, reviewing, and editing. The contributions of Dr. Ye included conceptualization, data curation, investigation, methodology, validation, and writing the original draft. The contributions of Dr. Duan included conceptualization, data curation, and investigation. The contributions of Dr. Yuan included data curation and investigation. Dr. Duan was our teammate and has graduated from Sichuan University last June 2018. However, when we submitted the original manuscript in August, 2018, we forgot to add his name.

#### Acknowledgments

The authors thank Elsevier’s Webshop for editing and polishing this paper. This work was supported by the International Cooperation Project of Sichuan Province (Grant No. 2018HH0082); the National Nature Science Foundation of China (Grant No. 41772321); the Science and Technology Planning Project of Sichuan Province, China (No. 2017TD0007); and the National Key R&D Program of China (Grant No. 2018YFC1505004).