# A Fractal Model for Characterizing Hydraulic Properties of Fractured Rock Mass under Mining Influence

Authors: Xiaoli Liu; Tao Liang; Sijing Wang; Kumar Nawnit
Journal: Geofluids
 (2019)

---

## Abstract

In this paper, two basic assumptions are introduced: (1) The number and length distribution of fractures in fractured rock mass are in accordance with the fractal law. (2) Fluid seepage in the fractures satisfies the cubic law. Based on these two assumptions, the fractal model of parallel seepage and radial seepage in fractured rock mass is established, and the seepage tensor of fracture network which reflects the geometric characteristics and fractal characteristics of fracture network under two kinds of seepage is derived. The influence of fracture geometry and fractal characteristics on permeability is analyzed, and the validity and accuracy of the model are verified by comparing the calculated results of the theoretical model and physical model test. The results show that the permeability coefficientK of fracture network is a function of the geometric (maximum crack length Lmax, fractured horizontal projection length L0, diameter calculation section porosity Φ, fracture strike α, and fracture angle θ) and fractal characteristics (fracture network fractal dimension Df and seepage flow fractal dimension DT). With the increase of fractal dimension Df, the permeability coefficient increases. With the increase of DT, the permeability coefficient decreases rapidly. And the larger the Df (Df>1.5), the greater the change of permeability coefficient K with DT.


---

## Body

## 1. Introduction

A large number of fractures are distributed in the natural rock mass and artificially disturbed rock mass, which provide channels for fluid seepage. To determine the permeability characteristic of fractured rock is important in the field of geology, geotechnical engineering, petrochemical resource exploitation, groundwater resource development and protection, nuclear waste disposal, and storage of carbon dioxide and others [1–9]. Generally, it is assumed that the fluid flows only in the interconnected fractures, and the permeability characteristics of equivalent fracture network is determined by analyzing the permeability characteristics of the different characteristic fractures [10, 11]. This method is called the discrete fracture network (DFN) method, which has been widely used and developed in recent decades [12–16].The permeability characteristics of the fracture network mainly depend on the fracture characteristics (spatial distribution, density, connectivity, etc.) at the macroscopic level and the fracture characteristics (length, gap width, directivity, roughness, etc.) at the microscopic level. Liu et al. [8] reviewed the current research on the influence of geometrical characteristics of fractured rock mass on the permeability of two-dimensional discrete fracture network and summed up nine parameters which have great influence on the permeability of fracture network. It is the length and distribution of fractures, the width and distribution, fractured surface roughness, dead fracture, fracture cross point, hydraulic gradient, stress condition, anisotropy, and size effect and listed the analytical expressions of the relationship between the relevant fracture parameters and permeability. de Dreuzy and Philippe [17] studied the effect of fracture length and gap width distribution on the infiltration characteristics of two-dimensional random fracture networks. Based on the fractal network statistics and fractal characteristics, a practical method to determine the permeability of fracture network was proposed by Jafari and Babadagli [18, 19]. Sensitivity analysis of the permeability of each parameter was carried out. The results showed that the fracture density and length have the greatest influence on the permeability of fracture network. The permeability of fracture network with exponential distribution was studied by Rossen and Gu [20]. It is found that the single large scale fracture plays a major role in the permeability of the fracture network. The permeability of the study area increases with the increase of the study size. Under normal circumstances, it is very difficult to quantitatively describe the fracture characteristics of the fracture network due to its complexity [21, 22]. Studies have shown that natural fractures have fractal properties that have the same effect on the permeability of fracture networks as well as fracture geometries [18, 23–26].A series of studies have been made on the influence of geometric characteristics and fractal features of fractures on permeability. A two-dimensional random fracture network model considering the characteristics of fracture network (density, length, gap width, directionality, and connectivity) was established. The fracture network statistics and the relationship between fractal characteristics and fracture network permeability were obtained through the multiple regression analysis and artificial neural network processing. An estimation method was proposed to estimate the permeability of multilayer complex fracture network, in order to obtain a good estimate result, which was recommended for comprehensive utilization of drill core data (1D), bedrock exposed data (2D), and drilling test data (3D). Another key factor that affects the permeability of fracture network is the connectivity between fractures. Based on the fractal geometry and percolation theory, the two-dimensional fracture network connectivity was defined, and the relationship between fracture network permeability and fractal dimension and dimensionless fracture density and percolation threshold was given and analyzed. It is pointed out that the estimated value obtained by fractal dimension of fractured network was more accurate than that of fractal dimension of fracture intersections, fracture connectivity rate, and intersection fractal dimension in the X-direction and Y-direction.In recent years, many scholars have established the fractal model of pore and fractured media based on the fractal characteristics of porosity and fractures and have made great achievements in studying fluid seepage, solute transport, and heat conduction [27–32]. Zheng and Yu [33] deduced the seepage characteristics of the gas in the porous media model composed of porous rock and fractal tree fractures. It is shown that the pore fractal dimension, pore bending degree, porosity, the ratio of maximum pore diameter to fractal tree fracture network length, diameter ratio, bifurcation angle, and bifurcation level are the key factors that have significant effects on gas permeability. Yun et al. [34] analyzed the plane radial seepage and parallel seepage of Newtonian fluid in porous media and deduced the permeability coefficient and flow rate and velocity expression in the two seepage cases. The fractal theory and the Monte Carlo method were used to establish the probability model of radial percolation in porous media by Xu et al. [35]. The results showed that the effective radial permeability coefficient decreased drastically with the radial distance increase and the porosity and pore fractal dimension of the radial seepage interface had a significant effect on the effective permeability coefficient. Based on the fractal geometry theory and the laminar cube law of fluid in fractures, a fractal model of seepage in fractured rock mass was deduced by Miao et al. [36]. The theoretical model showed that the permeability coefficient of fractured rock mass was a function of fracture fractal dimension, porosity, fracture density, maximum fracture length, gap width, fracture direction, and inclination angle. Miao et al. [37] assumed that the fractures in the fractured rock mass were randomly distributed, the fluid flow in the fracture was in accordance with the law of cube, the pores in the rock mass connected into tortuous channels, the fracture length and pore diameter distribution had fractal characteristics, and a fractal model of seepage flow in a two-hole media model was established. The Monte Carlo method was used to generate the fractured network, and a fractal flow model was established to reflect the geometrical characteristics of fractured rock mass. It is revealed that the bending degree, the gap width and the random number that reflects fractal regularity of fracture length distribution, had a significant influence on the permeability of fracture network [38]. The seepage fractal model of discrete fracture network was established by using the fractal dimension of fracture geometrical distribution and the fractal dimension of flow line which reflects the surface roughness of fractures. By simulating and calculating the equivalent permeability coefficient of fluid flowing through different geometric characteristics of discrete fractures, it was found that when the fractal dimension of the fractured network is less than 1.5, the seepage of the fracture network is mainly controlled by the small fracture with the length less than the width of the fracture network. With the increase of fractal dimension, the effect of long fractures on fracture network seepage is increasing [7, 8].The permeability of fractured rock mass is influenced by the heterogeneity of fracture direction, inclination, and length distribution, which lead to the fact that the permeability of the fractured rock mass is directional, and many scholars have carried on fruitful research in this area [9, 10, 12, 13, 39, 40]. Based on the previous research, this paper is aimed at establishing a parallel seepage and radial seepage fractal model of fractured rock mass considering the effect of fracture surface roughness in Section 2 and deducing the seepage tensor of the fracture network which reflects the geometrical and fractal characteristics of fracture network in Section 3. The influence of fracture geometry and fractal characteristics on the permeability is analyzed in Section 4. The validity and accuracy of the model are verified by comparing the theoretical model and the physical model water injection test in Section 5.


## 2. Fractal Characteristics of Fractured Network of Rock Mass

Studies have shown that the cumulative surface area distribution of the Earth’s surface islands is subject to power distribution [36, 41], i.e.,
(1)NA>a∝a−D/2,where N is the total number of islands with area A greater than constant a and D is the fractal dimension representing the area distribution of the island. On this basis, amax is used to represent the area of the largest island [42]. Equation (2) is given:
(2)NA>a=amaxaD/2.The relationship between the equivalent gap width and the crack length can be expressed [36, 43–45]:
(3)e=βln,where e is the equivalent gap width; β is the proportional coefficient, which is related to the mechanical properties of the surrounding rock and in the range of 0.001-0.1 [37, 44]; l is the fracture length; and n is a constant that reflects the fracture characteristics, which ranges from 0.5 to 2.0 [37, 44, 45]. When n=1, it indicates that the equivalent gap width and the length of the fracture are linearly distributed. The fractured network has self-similarity and fractal characteristics [37, 38, 43, 45]. Equation (3) can be written as follows:
(4)e=βl.Many studies have shown that the fracture length distribution satisfies the fractal law [15, 36, 45–49]. Therefore, equation (2), which describes the distribution law of island area, is used to describe the distribution of fractures area in fractured media:
(5)NA≥a=emaxlmaxelDf/2,where emax and lmax are the maximum gap width and maximum length, respectively, and e and l are the fracture width and length, respectively. Equation (4) is substituted into equation (5), and equation (6) is given as follows:
(6)NL≥l=lmaxlDf,where Df is the fractal dimension of the fracture length, for the two-dimensional problem 0<Df<2 and for the three-dimensional problem 0<Df<3. By substituting l in equation (6) with lmin, the total number of fractures in the fracture network can be expressed as follows:
(7)NtL>lmin=lmaxlminDf.In general, the number of fractures in the fracture network is large, so equation (6) can be approximated as a continuous differential equation and equation (8) can be obtained by solving differential l in equation (6):
(8)−dNl=DflmaxDfl−Df+1dl.Equation (8) indicates the number of fractures in [l,l+dl], and the negative sign indicates that the number of fractures decreases with increasing fracture length. Equation (8) divided by equation (7) is equal to the following equation:
(9)−dNNt=DflminDfl−Df+1dl=fldl,where fl=DflminDfl−Df+1 is the probability density function, which is satisfied by equation (10) by the definition of the probability function:
(10)∫−∞+∞fldl=∫lminlmaxfldl=1−lminlmaxDf≡1.Therefore, equation (11) is given as follows:
(11)lminlmaxDf≅0.In general,lmin<<lmax in equation (11) is a necessary condition for the fractal network to show fractal characteristics. Many researchers [22, 38] used lmin/lmax≤0.001 as a threshold for the fractal model considering the two-dimensional fracture network seepage. This criterion is also used in this paper.The relationship between porosity and fractal dimension can be expressed [37, 50, 51]:
(12)Df=dE+ln∅lnlmax/lmin,where Df is the fractal dimension of the fracture length; Φ is the porosity of the fracture network; lmin and lmax are the minimum and maximum values of the fracture length; dE is the European dimension, in which for the two-dimensional problem dE=2 and for the three-dimensional problem dE=3.According to the definition of porosity, the following equation is available:(13)ϕ=AfA0,where Φ is the porosity of the fracture network, A0 is the area of the cross section, and Af is the pore area.
(14)Af=−∫lminlmaxeldNl=βDflmax22−Df1−lminlmax2−Df.By solving simultaneously equations (12) and (14), the following equation becomes available:
(15)Af=βDflmax22−Df1−ϕ.By substituting equation (13) into equation (15), the following equation becomes available:
(16)A0=βDflmax22−Df1−ϕϕ.


## 3. Fractal Model of Permeability of Fractured Rock Mass

### 3.1. Formula Derivation

In a three-dimensional space, the direction of the fracture is determined by the direction and inclination, as shown in Figure1. The x-axis and the y-axis are the two coordinate directions of the horizontal plane, which is in the same direction as the geodetic coordinate; that is, the x-axis is N and the y-axis is E; α is the angle between the trend line and the y-axis, θ is the angle between the fracture plane and the horizontal plane (inclination), l is the fracture length, e is the fracture width, Lft is the seepage length of the fluid in the fracture considering the fracture roughness, Lf is the linear length in the direction of fracture seepage, Jf is the tangential component of the hydraulic gradient in the fracture, and Jn is the normal component of the hydraulic gradient in the fracture. L0 is the linear length of Lf projected to the x-axis direction.Figure 1
Three-dimensional fracture seepage diagram.Usually, there are a large number of cracks in the fractured rock mass, and it is impossible to determine the direction and inclination of the fractures one by one. Studies have shown that the directionality of many fractures in a given area is not exactly the same but can usually show a tendency [36, 52]. For example, the results of 1878 fractures in the study showed an average inclination of 70° and an average orientation of N-S [53]. Furthermore, the fracture network studied at present is only caused by mining. Compared with the natural fracture networks affected by complex geological factors, the fracture network studied here is more regular. Therefore, we take the trend and dip angle in the calculation model for a certain range of statistical average.It is assumed that the flow of fluid in the fracture can be described by the cubic law [39, 44, 54]:
(17)ql=e3l12μΔPLft,where ql is a single fracture flow, μ is the dynamic viscosity coefficient, e is the fracture width, l is the fracture length, Lft is the length of the seepage of the fluid in the fracture considering the fracture roughness, and ΔP is the pressure difference across the fracture.Because of the rough surface of the fracture, the fluid flow path in the fracture is a curve, which leads to the extension of the flow path and the decrease of the effective flow capacity, as shown in Figure1. The relationship between Lft and Lf can be expressed as follows:
(18)Lft=e1−DTLfDT,where DT is the fractal dimension of the flow line of the seepage flow. DT reflects the nonlinearity of the streamline, and the streamline is a straight line when DT=1, Lft=Lf.By substituting equations (4) and (18) into equation (17), the following equation becomes available:
(19)ql=β2+DTl3+DT12μΔPLfDT=β2+DTl3+DT12μDTLfDT−1dPdLf.In a three-dimensional space, the hydraulic gradient is divided into normal and tangential components along the fracture, as shown in Figure1. Only the tangential hydraulic gradient produces the seepage flow, that is,
(20)ql=β2+DTl3+DT12μDTLfDT−1Jf.By substitutingL0=Lfsina into equation (20), the following equation becomes available:
(21)ql=β2+DTl3+DT12μDTL0sinα1−DTJf.It is assumed that the angle between the normal direction of the fracture surface and the coordinate axes isα1, α2, and α3, respectively. The normal direction vector of the fracture surface can be expressed as
(22)nn=cosα1i+cosα2j+cosα3k.At the same time, the hydraulic gradient vector is decomposed along the fracture plane normal and tangential components [55]:
(23)J=Jf+Jn.Taking equation (24) into consideration,
(24)Jn=J∙nnnn.The hydraulic gradient is decomposed in the direction of the coordinate axis; equation (25) is available.
(25)J=Jx+Jy+Jz=Jxi+Jyj+Jzk.According to theorem of vector operation, equation (26) is available.
(26)Jf=J−Jn.Equations (22), (23), (24), (25), and (26) are simultaneously solved; equation (27) is available.
(27)Jf=Jx1−cosα1cosα1−Jycosα2cosα1−Jzcosα3cosα1i+−Jxcosα1cosα2+Jy1−cosα2cosα2−Jzcosα3cosα2j+−Jxcosα1cosα3−Jycosα2cosα3+Jz1−cosα3cosα3.Equation (27) is written in a matrix form:
(28)Jf=1−cos2α1−cosα2cosα1−cosα3cosα1−cosα1cosα21−cos2α2−cosα3cosα2−cosα1cosα3−cosα2cosα31−cos2α3JxJyJz.The total seepage flow can be obtained by integrating seepage flow of single fracture along the fracture length on the calculated cross section.(29)Q=−∫lminlmaxqldNl.By substituting equations (8), (21), and (28) into equation (29), equation (30) becomes available:
(30)Q=∫lminlmaxβ2+DTl3+DT12μDTL0sinα1−DTJfDflmaxDfl−Df+1dl=β2+DTlmax3+DT12μDTL0sinα1−DTDf3+DT−Df1−lminlmax3+DT−Df1−cos2α1−cosα2cosα1−cosα3cosα1−cosα1cosα21−cos2α2−cosα3cosα2−cosα1cosα3−cosα2cosα31−cos2α3JxJyJz.Consideringlmin<<lmax, 1<DT<2, 1<Df<2, lmin/lmax3+DT−Df<<1, equation (30) is simplified and rewritten as a matrix:
(31)qxqyqz=β2+DTlmax3+DT12μDTL0sinα1−DTDf3+DT−Df1−cos2α1−cosα2cosα1−cosα3cosα1−cosα1cosα21−cos2α2−cosα3cosα2−cosα1cosα3−cosα2cosα31−cos2α3JxJyJz.Figure1 shows that the relationship between the direction cosine of fracture normal plane and the direction and inclination of the geodetic coordinate is as follows:
(32)cosα1=−sinθcosα,cosα2=−sinθsinα,cosα3=cosθ.By substituting equation (32) into equation (31), equation (33) can be obtained as follows:
(33)qxqyqz=β2+DTlmax3+DT12μDTL0sinα1−DTDf3+DT−Df1−sin2θcos2α−sin2θsinαcosαsinθcosθcosα−sin2θsinαcosα1−sin2θsin2αsinθcosθsinαsinθcosθcosαsinθcosθsinαsin2θJxJyJz.The permeability coefficient can be expressed as(34)K=μQA0J.By substituting equations (16) and (33) into equation (34), equation (35) is obtained as follows:
(35)K=β1+DTlmax1+DT12DTL0sinα1−DT2−Df3+DT−Df∅1−∅1−sin2θcos2α−sin2θsinαcosαsinθcosθcosα−sin2θsinαcosα1−sin2θsin2αsinθcosθsinαsinθcosθcosαsinθcosθsinαsin2θ.


### 3.2. Radial Seepage

In a three-dimensional space, the radial flow diagram is shown in Figure2. Map coordinates and parameter definitions are consistent with Figure 1. The x- and y-axes are two coordinate directions, which are consistent with the direction of geodetic coordinates. It is assumed the x-axis is N, the y-axis is E. α is the angle between the strike line and the y-axis, and θ is the angle between the fracture plane and the horizontal plane (inclination). l is the fracture length, e is the fracture width, ra is the seepage length of the fluid in the fracture considering the fracture roughness, r is taken as the straight line length along fracture seepage direction, rw is radius of the hole, Jf is the tangential component of the hydraulic gradient in the fracture, and Jn is the normal component of the hydraulic gradient in the fracture.Figure 2
Radial flow diagram: (a) elevation view and (b) section view.
(a)
(b)It is assumed that the flow of fluid in the fracture can be described by the cubic law [39, 44, 54]. Equation (17) can be rewritten as follows under the condition of radial flow:
(36)ql=e3l12μdpdra,where ql is a single fracture flow, μ is the dynamic viscous coefficient, e is the fracture width, l is the fracture length, and p is the water pressure. Due to the rugged surface of the fracture, the fluid seepage path in the fracture is a curve, resulting in a decrease in the flow path and an effective overcurrent capacity, as shown in Figure 2.The relation betweenra and r can be rewritten according to equation (18) as follows:
(37)ra=e1−DTrDT.Equation (38) can be obtained by substituting equations (4) and (37) into equation (36) as follows:
(38)ql=β2+DTl3+DT12μDTrDT−1dpdr.Similarly, the hydraulic gradient vector will be decomposed along the fracture plane in normal and tangential components [55]. Only the hydraulic gradient along the tangential direction of the fracture will produce seepage flow [12], i.e.,
(39)ql=β2+DTl3+DT12μDTrDT−1Jf.On the inner side of the radial hole wall, the seepage flow of single fracture is integrated along the fracture length, and the total seepage flow can be obtained.(40)Q=2πrwhA0∫lminlmaxqldN=πhrwβ2+DTlmax3+DT6μDTA0rDT−1Df3+DT−Df1−sin2θcos2α−sin2θsinαcosαsinθcosθcosα−sin2θsinαcosα1−sin2θsin2αsinθcosθsinαsinθcosθcosαsinθcosθsinαsin2θJxJyJz.Equation (16) is substituted into equation (40) and is rewritten in a matrix form:
(41)qxqyqz=πhrwβ1+DTlmax1+DT6μDTrDT−12−Df3+DT−Df∅1−∅1−sin2θcos2α−sin2θsinαcosαsinθcosθcosα−sin2θsinαcosα1−sin2θsin2αsinθcosθsinαsinθcosθcosαsinθcosθsinαsin2θJxJyJz.Combining equations (41) and (34) and A0 being replaced by Ar=2πrh, equation (42) is obtained as follows:
(42)Kr=rwβ1+DTlmax1+DT12DTrDT2−Df3+DT−DfΦ1−Φ1−sin2θcos2α−sin2θsinαcosαsinθcosθcosα−sin2θsinαcosα1−sin2θsin2αsinθcosθsinαsinθcosθcosαsinθcosθsinαsin2θ.The permeability coefficient can be expressed as follows at the inner side of the radial hole wall:(43)Kw=β1+DTlmax1+DT12DTrwDT−12−Df3+DT−DfΦ1−Φ1−sin2θcos2α−sin2θsinαcosαsinθcosθcosα−sin2θsinαcosα1−sin2θsin2αsinθcosθsinαsinθcosθcosαsinθcosθsinαsin2θ.From equations (42) and (43), the dimensionless permeability coefficient is obtained as follows:
(44)K+=KrKw=rwrDT.Equation (44) shows that the dimensionless permeability coefficient is closely related to DT and r/rw. It decreases with the increase of DT and decreases with the increase of r/rw, as shown in Figure 3.Figure 3
The relationship betweenK+ and r/rw.


### 3.3. Parameter Impact Analysis

This section mainly analyzes the influence of the parameters of equations (35) and (42) on the permeability coefficient K, and the values of each parameter are shown in Table 1.Table 1
List of calculation parameters.
Parameter
Value
Parameter description
Df
1.1-1.8
Fractal dimension of fracture
DT
1.0-1.7
Fractal dimension of seepage streamline
β
0.001-0.1
Ratio of fracture width and length
α
0-π
Trend
θ
0-π/2
Dip angle
lmax
3-30 m
Maximum fracture length
lmin/lmax
0.001
Minimum fracture length/maximum fracture length
L0
3-180 m
Computational element lengthFigure4 shows the relationship between the fractal dimension Df and the permeability coefficient Kxx in the case of different DT, and the other parameters are β=0.01, θ=π/6, and α=π/2. As shown in Figure 4, the permeability coefficient increases with the increase of fractal dimension Df. And with the increase of DT, the permeability coefficient Kxx decreases. The results show that the increase of the fracture surface roughness leads to the increase of the flow path and the decrease of flow velocity.Figure 4
The relationship betweenKxx and Df (β=0.01, θ=π/6, and α=π/2).Figure5 shows the relationship between DT and permeability coefficient Kxx in the case of different Df, and the other parameters are β=0.01, θ=π/6, and α=π/2. As shown in Figure 5, with the increase of DT, the permeability coefficient decreases rapidly. The overall permeability coefficient increases with the increase of Df, and the larger the Df (Df>1.5), the greater the change rate of Kxx with DT, which means that when Df>1.5, the weakening effect of the fracture surface roughness on the permeability of the fracture network increases rapidly.Figure 5
The relationship betweenKxx and DT (β=0.01, θ=π/6, and α=π/2).As can be seen from Figures4 and 5, Df and DT have a significant effect on the permeability of the fracture network, but they have different physical meanings and influences. The larger the value of Df, the more complex the fractured network, the higher the fracture density [18], and the greater the permeability coefficient. DT reflects the flexural degree of the fluid seepage path in the fracture due to fracture surface roughness. With the increase of DT, the seepage path is lengthened and the effective overcurrent capacity of the fracture is reduced, which leads to the decrease of the permeability coefficient of the fracture network.Figure6 shows the relationship between porosity Φ and permeability coefficient Kxx in the case of different DT, and the other parameters are β=0.01, θ=π/6, and α=π/2. As shown in the figure, the permeability coefficient increases with the increase of porosity Φ. In general, it can be seen from Equation (12) and Figure 7 that, assuming lmin/lmax=0.001, Φ is positively correlated with Df, and thus the increase of porosity Φ and Df increase have identity, and all reflect the increase of seepage cross section. It also shows in Figure 6 that the overall permeability coefficient decreases with the increase of DT, and the larger the DT (DT>1.2), the smaller the change rate of Kxx with Φ, which indicates that when DT>1.2, the roughness of the fracture surface will greatly offset the effect of increased permeability of the fracture network due to the increase of the seepage area.Figure 6
The relationship betweenKxx and Φ (β=0.01, θ=π/6, and α=π/2).Figure 7
The relationship betweenΦ and Df.Figure8 shows the relationship between fractal dimension Df and permeability coefficient Kxx of the fracture network under the condition of the change of the fracture dip angle θ, and the other parameters are DT=1.1, β=0.01, and α=π/4. As shown in the figure, the permeability coefficient increases with the increase of porosity Df, which is the same as that of Figure 4. At the same time, as the fracture dip angle θ increases, the permeability coefficient Kxx decreases.Figure 8
The relationship betweenKxx and Df (DT=1.1, β=0.01, and α=π/4).Figure9 shows the relationship between the fracture dip angle θ and the permeability coefficient Kxx, and the other parameters are Df=1.4, DT=1.1, β=0.01, and α=π/4. As shown in Figure 9, the permeability coefficient Kxx decreases with the increasing of fracture dip angle θ from 0 to π/2. Figure 10 shows the relationship between the fracture strike α and the permeability coefficient Kxx, and the other parameters are Df=1.4, DT=1.1, β=0.01, and θ=π/4. As shown in Figure 10, the permeability coefficient Kxx increases as the dip angle α of fracture increases from π/18 to π/2. The variation of the permeability coefficient component with the fracture strike and dip angle reflects the nonuniformity and seepage directivity of fractured rock mass.Figure 9
The relationship betweenKxx and θ (Df=1.4, DT=1.1, β=0.01, and α=π/4).Figure 10
The relationship betweenKxx and α (Df=1.4, DT=1.1, β=0.01, and α=π/4).Figure11 shows the relationship between the maximum fracture length Lmax and permeability coefficient Kxx under condition of different Df. The other parameters are DT=1.1, β=0.01, θ=π/6, and α=π/2. As shown in Figure 11, the permeability coefficient increases with the increase of Lmax. At the same time, with the increase of Df, the overall permeability coefficient will increase, and with the increase of Df (Df>1.5), the change rate of Kxx with Lmax get sharply greater, indicating that when Df>1.5, with the increase of fractal dimension, the influence of long fractures on fracture network seepage is increasing, which is consistent with the results of Li et al. [22].Figure 11
The relationship betweenLmax and Kxx (DT=1.1, β=0.01, θ=π/6, and α=π/2).Figure12 shows the relationship between β and the permeability coefficient Kxx, and the other parameters are Df=1.4, DT=1.01, θ=π/6, and α=π/2. As shown in Figure 12, the permeability coefficient increases with the increase of β. Considering equation (4), this also reflects the relationship between the fracture width and the permeability coefficient to a certain extent.Figure 12
The relationship betweenβ and Kxx (Df=1.4, DT=1.01, θ=π/6, and α=π/2).Figure13 shows the relationship between L0 and permeability coefficient Kxx, and the other parameters are Df=1.4, DT=1.1, β=0.01, θ=π/4, and α=π/2. As shown in Figure 13, the permeability coefficient decreases gradually with the increase of L0. This is because the increase in L0 is the growth of the seepage path, which means that the seepage resistance increases and the hydraulic gradient decreases.Figure 13
The relationship betweenL0 and Kxx (Df=1.4, DT=1.1, β=0.01, θ=π/4, and α=π/2).


## 3.1. Formula Derivation

In a three-dimensional space, the direction of the fracture is determined by the direction and inclination, as shown in Figure1. The x-axis and the y-axis are the two coordinate directions of the horizontal plane, which is in the same direction as the geodetic coordinate; that is, the x-axis is N and the y-axis is E; α is the angle between the trend line and the y-axis, θ is the angle between the fracture plane and the horizontal plane (inclination), l is the fracture length, e is the fracture width, Lft is the seepage length of the fluid in the fracture considering the fracture roughness, Lf is the linear length in the direction of fracture seepage, Jf is the tangential component of the hydraulic gradient in the fracture, and Jn is the normal component of the hydraulic gradient in the fracture. L0 is the linear length of Lf projected to the x-axis direction.Figure 1
Three-dimensional fracture seepage diagram.Usually, there are a large number of cracks in the fractured rock mass, and it is impossible to determine the direction and inclination of the fractures one by one. Studies have shown that the directionality of many fractures in a given area is not exactly the same but can usually show a tendency [36, 52]. For example, the results of 1878 fractures in the study showed an average inclination of 70° and an average orientation of N-S [53]. Furthermore, the fracture network studied at present is only caused by mining. Compared with the natural fracture networks affected by complex geological factors, the fracture network studied here is more regular. Therefore, we take the trend and dip angle in the calculation model for a certain range of statistical average.It is assumed that the flow of fluid in the fracture can be described by the cubic law [39, 44, 54]:
(17)ql=e3l12μΔPLft,where ql is a single fracture flow, μ is the dynamic viscosity coefficient, e is the fracture width, l is the fracture length, Lft is the length of the seepage of the fluid in the fracture considering the fracture roughness, and ΔP is the pressure difference across the fracture.Because of the rough surface of the fracture, the fluid flow path in the fracture is a curve, which leads to the extension of the flow path and the decrease of the effective flow capacity, as shown in Figure1. The relationship between Lft and Lf can be expressed as follows:
(18)Lft=e1−DTLfDT,where DT is the fractal dimension of the flow line of the seepage flow. DT reflects the nonlinearity of the streamline, and the streamline is a straight line when DT=1, Lft=Lf.By substituting equations (4) and (18) into equation (17), the following equation becomes available:
(19)ql=β2+DTl3+DT12μΔPLfDT=β2+DTl3+DT12μDTLfDT−1dPdLf.In a three-dimensional space, the hydraulic gradient is divided into normal and tangential components along the fracture, as shown in Figure1. Only the tangential hydraulic gradient produces the seepage flow, that is,
(20)ql=β2+DTl3+DT12μDTLfDT−1Jf.By substitutingL0=Lfsina into equation (20), the following equation becomes available:
(21)ql=β2+DTl3+DT12μDTL0sinα1−DTJf.It is assumed that the angle between the normal direction of the fracture surface and the coordinate axes isα1, α2, and α3, respectively. The normal direction vector of the fracture surface can be expressed as
(22)nn=cosα1i+cosα2j+cosα3k.At the same time, the hydraulic gradient vector is decomposed along the fracture plane normal and tangential components [55]:
(23)J=Jf+Jn.Taking equation (24) into consideration,
(24)Jn=J∙nnnn.The hydraulic gradient is decomposed in the direction of the coordinate axis; equation (25) is available.
(25)J=Jx+Jy+Jz=Jxi+Jyj+Jzk.According to theorem of vector operation, equation (26) is available.
(26)Jf=J−Jn.Equations (22), (23), (24), (25), and (26) are simultaneously solved; equation (27) is available.
(27)Jf=Jx1−cosα1cosα1−Jycosα2cosα1−Jzcosα3cosα1i+−Jxcosα1cosα2+Jy1−cosα2cosα2−Jzcosα3cosα2j+−Jxcosα1cosα3−Jycosα2cosα3+Jz1−cosα3cosα3.Equation (27) is written in a matrix form:
(28)Jf=1−cos2α1−cosα2cosα1−cosα3cosα1−cosα1cosα21−cos2α2−cosα3cosα2−cosα1cosα3−cosα2cosα31−cos2α3JxJyJz.The total seepage flow can be obtained by integrating seepage flow of single fracture along the fracture length on the calculated cross section.(29)Q=−∫lminlmaxqldNl.By substituting equations (8), (21), and (28) into equation (29), equation (30) becomes available:
(30)Q=∫lminlmaxβ2+DTl3+DT12μDTL0sinα1−DTJfDflmaxDfl−Df+1dl=β2+DTlmax3+DT12μDTL0sinα1−DTDf3+DT−Df1−lminlmax3+DT−Df1−cos2α1−cosα2cosα1−cosα3cosα1−cosα1cosα21−cos2α2−cosα3cosα2−cosα1cosα3−cosα2cosα31−cos2α3JxJyJz.Consideringlmin<<lmax, 1<DT<2, 1<Df<2, lmin/lmax3+DT−Df<<1, equation (30) is simplified and rewritten as a matrix:
(31)qxqyqz=β2+DTlmax3+DT12μDTL0sinα1−DTDf3+DT−Df1−cos2α1−cosα2cosα1−cosα3cosα1−cosα1cosα21−cos2α2−cosα3cosα2−cosα1cosα3−cosα2cosα31−cos2α3JxJyJz.Figure1 shows that the relationship between the direction cosine of fracture normal plane and the direction and inclination of the geodetic coordinate is as follows:
(32)cosα1=−sinθcosα,cosα2=−sinθsinα,cosα3=cosθ.By substituting equation (32) into equation (31), equation (33) can be obtained as follows:
(33)qxqyqz=β2+DTlmax3+DT12μDTL0sinα1−DTDf3+DT−Df1−sin2θcos2α−sin2θsinαcosαsinθcosθcosα−sin2θsinαcosα1−sin2θsin2αsinθcosθsinαsinθcosθcosαsinθcosθsinαsin2θJxJyJz.The permeability coefficient can be expressed as(34)K=μQA0J.By substituting equations (16) and (33) into equation (34), equation (35) is obtained as follows:
(35)K=β1+DTlmax1+DT12DTL0sinα1−DT2−Df3+DT−Df∅1−∅1−sin2θcos2α−sin2θsinαcosαsinθcosθcosα−sin2θsinαcosα1−sin2θsin2αsinθcosθsinαsinθcosθcosαsinθcosθsinαsin2θ.


## 3.2. Radial Seepage

In a three-dimensional space, the radial flow diagram is shown in Figure2. Map coordinates and parameter definitions are consistent with Figure 1. The x- and y-axes are two coordinate directions, which are consistent with the direction of geodetic coordinates. It is assumed the x-axis is N, the y-axis is E. α is the angle between the strike line and the y-axis, and θ is the angle between the fracture plane and the horizontal plane (inclination). l is the fracture length, e is the fracture width, ra is the seepage length of the fluid in the fracture considering the fracture roughness, r is taken as the straight line length along fracture seepage direction, rw is radius of the hole, Jf is the tangential component of the hydraulic gradient in the fracture, and Jn is the normal component of the hydraulic gradient in the fracture.Figure 2
Radial flow diagram: (a) elevation view and (b) section view.
(a)
(b)It is assumed that the flow of fluid in the fracture can be described by the cubic law [39, 44, 54]. Equation (17) can be rewritten as follows under the condition of radial flow:
(36)ql=e3l12μdpdra,where ql is a single fracture flow, μ is the dynamic viscous coefficient, e is the fracture width, l is the fracture length, and p is the water pressure. Due to the rugged surface of the fracture, the fluid seepage path in the fracture is a curve, resulting in a decrease in the flow path and an effective overcurrent capacity, as shown in Figure 2.The relation betweenra and r can be rewritten according to equation (18) as follows:
(37)ra=e1−DTrDT.Equation (38) can be obtained by substituting equations (4) and (37) into equation (36) as follows:
(38)ql=β2+DTl3+DT12μDTrDT−1dpdr.Similarly, the hydraulic gradient vector will be decomposed along the fracture plane in normal and tangential components [55]. Only the hydraulic gradient along the tangential direction of the fracture will produce seepage flow [12], i.e.,
(39)ql=β2+DTl3+DT12μDTrDT−1Jf.On the inner side of the radial hole wall, the seepage flow of single fracture is integrated along the fracture length, and the total seepage flow can be obtained.(40)Q=2πrwhA0∫lminlmaxqldN=πhrwβ2+DTlmax3+DT6μDTA0rDT−1Df3+DT−Df1−sin2θcos2α−sin2θsinαcosαsinθcosθcosα−sin2θsinαcosα1−sin2θsin2αsinθcosθsinαsinθcosθcosαsinθcosθsinαsin2θJxJyJz.Equation (16) is substituted into equation (40) and is rewritten in a matrix form:
(41)qxqyqz=πhrwβ1+DTlmax1+DT6μDTrDT−12−Df3+DT−Df∅1−∅1−sin2θcos2α−sin2θsinαcosαsinθcosθcosα−sin2θsinαcosα1−sin2θsin2αsinθcosθsinαsinθcosθcosαsinθcosθsinαsin2θJxJyJz.Combining equations (41) and (34) and A0 being replaced by Ar=2πrh, equation (42) is obtained as follows:
(42)Kr=rwβ1+DTlmax1+DT12DTrDT2−Df3+DT−DfΦ1−Φ1−sin2θcos2α−sin2θsinαcosαsinθcosθcosα−sin2θsinαcosα1−sin2θsin2αsinθcosθsinαsinθcosθcosαsinθcosθsinαsin2θ.The permeability coefficient can be expressed as follows at the inner side of the radial hole wall:(43)Kw=β1+DTlmax1+DT12DTrwDT−12−Df3+DT−DfΦ1−Φ1−sin2θcos2α−sin2θsinαcosαsinθcosθcosα−sin2θsinαcosα1−sin2θsin2αsinθcosθsinαsinθcosθcosαsinθcosθsinαsin2θ.From equations (42) and (43), the dimensionless permeability coefficient is obtained as follows:
(44)K+=KrKw=rwrDT.Equation (44) shows that the dimensionless permeability coefficient is closely related to DT and r/rw. It decreases with the increase of DT and decreases with the increase of r/rw, as shown in Figure 3.Figure 3
The relationship betweenK+ and r/rw.


## 3.3. Parameter Impact Analysis

This section mainly analyzes the influence of the parameters of equations (35) and (42) on the permeability coefficient K, and the values of each parameter are shown in Table 1.Table 1
List of calculation parameters.
Parameter
Value
Parameter description
Df
1.1-1.8
Fractal dimension of fracture
DT
1.0-1.7
Fractal dimension of seepage streamline
β
0.001-0.1
Ratio of fracture width and length
α
0-π
Trend
θ
0-π/2
Dip angle
lmax
3-30 m
Maximum fracture length
lmin/lmax
0.001
Minimum fracture length/maximum fracture length
L0
3-180 m
Computational element lengthFigure4 shows the relationship between the fractal dimension Df and the permeability coefficient Kxx in the case of different DT, and the other parameters are β=0.01, θ=π/6, and α=π/2. As shown in Figure 4, the permeability coefficient increases with the increase of fractal dimension Df. And with the increase of DT, the permeability coefficient Kxx decreases. The results show that the increase of the fracture surface roughness leads to the increase of the flow path and the decrease of flow velocity.Figure 4
The relationship betweenKxx and Df (β=0.01, θ=π/6, and α=π/2).Figure5 shows the relationship between DT and permeability coefficient Kxx in the case of different Df, and the other parameters are β=0.01, θ=π/6, and α=π/2. As shown in Figure 5, with the increase of DT, the permeability coefficient decreases rapidly. The overall permeability coefficient increases with the increase of Df, and the larger the Df (Df>1.5), the greater the change rate of Kxx with DT, which means that when Df>1.5, the weakening effect of the fracture surface roughness on the permeability of the fracture network increases rapidly.Figure 5
The relationship betweenKxx and DT (β=0.01, θ=π/6, and α=π/2).As can be seen from Figures4 and 5, Df and DT have a significant effect on the permeability of the fracture network, but they have different physical meanings and influences. The larger the value of Df, the more complex the fractured network, the higher the fracture density [18], and the greater the permeability coefficient. DT reflects the flexural degree of the fluid seepage path in the fracture due to fracture surface roughness. With the increase of DT, the seepage path is lengthened and the effective overcurrent capacity of the fracture is reduced, which leads to the decrease of the permeability coefficient of the fracture network.Figure6 shows the relationship between porosity Φ and permeability coefficient Kxx in the case of different DT, and the other parameters are β=0.01, θ=π/6, and α=π/2. As shown in the figure, the permeability coefficient increases with the increase of porosity Φ. In general, it can be seen from Equation (12) and Figure 7 that, assuming lmin/lmax=0.001, Φ is positively correlated with Df, and thus the increase of porosity Φ and Df increase have identity, and all reflect the increase of seepage cross section. It also shows in Figure 6 that the overall permeability coefficient decreases with the increase of DT, and the larger the DT (DT>1.2), the smaller the change rate of Kxx with Φ, which indicates that when DT>1.2, the roughness of the fracture surface will greatly offset the effect of increased permeability of the fracture network due to the increase of the seepage area.Figure 6
The relationship betweenKxx and Φ (β=0.01, θ=π/6, and α=π/2).Figure 7
The relationship betweenΦ and Df.Figure8 shows the relationship between fractal dimension Df and permeability coefficient Kxx of the fracture network under the condition of the change of the fracture dip angle θ, and the other parameters are DT=1.1, β=0.01, and α=π/4. As shown in the figure, the permeability coefficient increases with the increase of porosity Df, which is the same as that of Figure 4. At the same time, as the fracture dip angle θ increases, the permeability coefficient Kxx decreases.Figure 8
The relationship betweenKxx and Df (DT=1.1, β=0.01, and α=π/4).Figure9 shows the relationship between the fracture dip angle θ and the permeability coefficient Kxx, and the other parameters are Df=1.4, DT=1.1, β=0.01, and α=π/4. As shown in Figure 9, the permeability coefficient Kxx decreases with the increasing of fracture dip angle θ from 0 to π/2. Figure 10 shows the relationship between the fracture strike α and the permeability coefficient Kxx, and the other parameters are Df=1.4, DT=1.1, β=0.01, and θ=π/4. As shown in Figure 10, the permeability coefficient Kxx increases as the dip angle α of fracture increases from π/18 to π/2. The variation of the permeability coefficient component with the fracture strike and dip angle reflects the nonuniformity and seepage directivity of fractured rock mass.Figure 9
The relationship betweenKxx and θ (Df=1.4, DT=1.1, β=0.01, and α=π/4).Figure 10
The relationship betweenKxx and α (Df=1.4, DT=1.1, β=0.01, and α=π/4).Figure11 shows the relationship between the maximum fracture length Lmax and permeability coefficient Kxx under condition of different Df. The other parameters are DT=1.1, β=0.01, θ=π/6, and α=π/2. As shown in Figure 11, the permeability coefficient increases with the increase of Lmax. At the same time, with the increase of Df, the overall permeability coefficient will increase, and with the increase of Df (Df>1.5), the change rate of Kxx with Lmax get sharply greater, indicating that when Df>1.5, with the increase of fractal dimension, the influence of long fractures on fracture network seepage is increasing, which is consistent with the results of Li et al. [22].Figure 11
The relationship betweenLmax and Kxx (DT=1.1, β=0.01, θ=π/6, and α=π/2).Figure12 shows the relationship between β and the permeability coefficient Kxx, and the other parameters are Df=1.4, DT=1.01, θ=π/6, and α=π/2. As shown in Figure 12, the permeability coefficient increases with the increase of β. Considering equation (4), this also reflects the relationship between the fracture width and the permeability coefficient to a certain extent.Figure 12
The relationship betweenβ and Kxx (Df=1.4, DT=1.01, θ=π/6, and α=π/2).Figure13 shows the relationship between L0 and permeability coefficient Kxx, and the other parameters are Df=1.4, DT=1.1, β=0.01, θ=π/4, and α=π/2. As shown in Figure 13, the permeability coefficient decreases gradually with the increase of L0. This is because the increase in L0 is the growth of the seepage path, which means that the seepage resistance increases and the hydraulic gradient decreases.Figure 13
The relationship betweenL0 and Kxx (Df=1.4, DT=1.1, β=0.01, θ=π/4, and α=π/2).


## 4. Experimental Verification

### 4.1. Physical Model Test

Figure14 shows a large physical model (partial) profile. According to the generalized stratigraphic parameters shown in Table 2, the similar materials arranged according to a similarity ratio are poured up and down, and the bottom pouring height of no. 2 coal seam and no. 5 coal seam was 2.25 m and 0.5 m, respectively. The thickness of the layer is 0.05 meters, the net length of working face is 2.1 meters, and the reserved pillar width is 0.3 meters. The displacement and stress changes of the overlying strata and the expansion and evolution of the fractures are studied under the conditions of different burial depth and mining speed. Figure 15 shows the distribution of the fractures of the overlying strata after the mining of no. 5 coal seam.Figure 14
Physical model (partial) profile.Table 2
Generalized stratigraphic parameters.
No.
Stratum
Thickness (m)
σ (MPa)
ρ (kg/m3)
1
Aeolian sand and loess
30
2000
2
Mudstone
15
34.86
2304
3
Fine sandstone and siltstone (aquifer no.1)
10
50.00
2550
4
Mudstone
20
25.22
2204
5
Fine sandstone and siltstone (aquifer no.2)
10
5.26
2133
6
Mudstone
15
25.22
2204
7
Coal seam no. 2
5
20.07
1233
8
Mudstone
15
34.86
2304
9
Fine sandstone and siltstone
20
5.26
2133
10
Mudstone
15
34.86
2304
11
Medium-thick grained sandstone (aquifer no.3)
20
50.00
2550
12
Mudstone
20
25.22
2204
13
Fine sandstone and siltstone
10
5.26
2133
14
Mudstone
10
25.22
2204
15
Medium-thick grained sandstone (aquifer no.4)
20
50.00
2550
16
Mudstone
20
25.22
2204
17
Medium-thick grained sandstone
20
5.26
2133
18
Coal seam no. 5
5
20.07
1233
19
Mudstone
15
34.86
2304
20
Medium-thick grained sandstone
15
50.00
2550
21
Mudstone
10
34.86
2304
22
Medium-thick grained sandstone
10
50.00
2550Figure 15
Distribution map of mining-induced fractures of no. 5 coal seam mining.


### 4.2. Water Injection Test in the Physical Model

In order to obtain the permeability parameters of the rock mass in different fracture zones and verify the fractal model of the fractured rock mass, the water injection test in the physical model was carried out after the physical model excavation was completed. In order to compare with the evolution of fracture expansion, the working face is divided into seven zones along the excavation direction, and the width of each zone is 30 cm. In order to reduce the interaction of the boundary effect and the water injection holes, the water injection holes were arranged in plum blossom in the middle of the partition, as shown in Figure16.Figure 16
Plane layout of water injection test points (units: mm).The location of water injection test sites is carried out at the top of the model in accordance with Figure16 when constructing, and then the hydraulic drill is set up to conduct drilling construction, as shown in Figures 17 and 18. The water injection test is carried out by the top-down stratified water injection method to obtain the permeability parameters of different plane positions and different elevation rocks in the mining face. Meanwhile, the corresponding fracture distribution images are obtained by drilling imaging an analyzer, and the fracture images are processed digitally to obtain the fracture expansion images, as shown in Figure 19. The fracture expansion images are analyzed to extract the fracture geometric characteristic parameters, and thus the fractal characteristic parameters, DT and Df, are calculated by the box counting method [51, 56–59].Figure 17
Water injection test profile.Figure 18
Bench drill and drilling imaging analyzer.Figure 19
Inner view of different holes and their fractures developed pattern in different depth.
(a)
Inner view in ZS3-6 (60-75 cm)
(b)
Fractures developed pattern in ZS3-6 (55-65 cm)
(c)
Fractures developed pattern in ZS3-6 (65-75 cm)
(d)
Inner view in ZS3-6 (175-191 cm)
(e)
Fractures developed pattern in ZS3-6 (175-191 cm)
(f)
Inner view in ZS3-6 (225-240 cm)
(g)
Fractures developed pattern in ZS3-6 (225-235 cm)
(h)
Fractures developed pattern in ZS3-6 (235-245 cm)
(i)
Inner view in ZS3-4 (170-185 cm)
(j)
Fractures developed pattern in ZS3-4 (175-185 cm)


## 4.1. Physical Model Test

Figure14 shows a large physical model (partial) profile. According to the generalized stratigraphic parameters shown in Table 2, the similar materials arranged according to a similarity ratio are poured up and down, and the bottom pouring height of no. 2 coal seam and no. 5 coal seam was 2.25 m and 0.5 m, respectively. The thickness of the layer is 0.05 meters, the net length of working face is 2.1 meters, and the reserved pillar width is 0.3 meters. The displacement and stress changes of the overlying strata and the expansion and evolution of the fractures are studied under the conditions of different burial depth and mining speed. Figure 15 shows the distribution of the fractures of the overlying strata after the mining of no. 5 coal seam.Figure 14
Physical model (partial) profile.Table 2
Generalized stratigraphic parameters.
No.
Stratum
Thickness (m)
σ (MPa)
ρ (kg/m3)
1
Aeolian sand and loess
30
2000
2
Mudstone
15
34.86
2304
3
Fine sandstone and siltstone (aquifer no.1)
10
50.00
2550
4
Mudstone
20
25.22
2204
5
Fine sandstone and siltstone (aquifer no.2)
10
5.26
2133
6
Mudstone
15
25.22
2204
7
Coal seam no. 2
5
20.07
1233
8
Mudstone
15
34.86
2304
9
Fine sandstone and siltstone
20
5.26
2133
10
Mudstone
15
34.86
2304
11
Medium-thick grained sandstone (aquifer no.3)
20
50.00
2550
12
Mudstone
20
25.22
2204
13
Fine sandstone and siltstone
10
5.26
2133
14
Mudstone
10
25.22
2204
15
Medium-thick grained sandstone (aquifer no.4)
20
50.00
2550
16
Mudstone
20
25.22
2204
17
Medium-thick grained sandstone
20
5.26
2133
18
Coal seam no. 5
5
20.07
1233
19
Mudstone
15
34.86
2304
20
Medium-thick grained sandstone
15
50.00
2550
21
Mudstone
10
34.86
2304
22
Medium-thick grained sandstone
10
50.00
2550Figure 15
Distribution map of mining-induced fractures of no. 5 coal seam mining.


## 4.2. Water Injection Test in the Physical Model

In order to obtain the permeability parameters of the rock mass in different fracture zones and verify the fractal model of the fractured rock mass, the water injection test in the physical model was carried out after the physical model excavation was completed. In order to compare with the evolution of fracture expansion, the working face is divided into seven zones along the excavation direction, and the width of each zone is 30 cm. In order to reduce the interaction of the boundary effect and the water injection holes, the water injection holes were arranged in plum blossom in the middle of the partition, as shown in Figure16.Figure 16
Plane layout of water injection test points (units: mm).The location of water injection test sites is carried out at the top of the model in accordance with Figure16 when constructing, and then the hydraulic drill is set up to conduct drilling construction, as shown in Figures 17 and 18. The water injection test is carried out by the top-down stratified water injection method to obtain the permeability parameters of different plane positions and different elevation rocks in the mining face. Meanwhile, the corresponding fracture distribution images are obtained by drilling imaging an analyzer, and the fracture images are processed digitally to obtain the fracture expansion images, as shown in Figure 19. The fracture expansion images are analyzed to extract the fracture geometric characteristic parameters, and thus the fractal characteristic parameters, DT and Df, are calculated by the box counting method [51, 56–59].Figure 17
Water injection test profile.Figure 18
Bench drill and drilling imaging analyzer.Figure 19
Inner view of different holes and their fractures developed pattern in different depth.
(a)
Inner view in ZS3-6 (60-75 cm)
(b)
Fractures developed pattern in ZS3-6 (55-65 cm)
(c)
Fractures developed pattern in ZS3-6 (65-75 cm)
(d)
Inner view in ZS3-6 (175-191 cm)
(e)
Fractures developed pattern in ZS3-6 (175-191 cm)
(f)
Inner view in ZS3-6 (225-240 cm)
(g)
Fractures developed pattern in ZS3-6 (225-235 cm)
(h)
Fractures developed pattern in ZS3-6 (235-245 cm)
(i)
Inner view in ZS3-4 (170-185 cm)
(j)
Fractures developed pattern in ZS3-4 (175-185 cm)


## 5. Analysis and Discussion

According to the field test conditions, the continuous head water injection test is used to determine the seepage per unit time in the case of constant water head and then calculate the permeability coefficient of the simulated strata according to equation (44).
(45)K=0.423QlHlg2Hr,where K is the permeability coefficient of the soil layer (cm/s), Q is the stable water flow (cm3/s), H is the test water head (cm), l is the length of the water injection section (cm), and r is the water injection test hole radius (cm).Fracture development images are analyzed in the previous section; then the fractal characteristic parameters are extracted, and the relevant parameters are substituted into equation (42). The calculated results are compared with the Kxx values calculated by equation (43), as shown in Table 3. And the fractal characteristics are calculated, as shown in Figure 20. The relevant parameters are substituted into equation (41); then the calculated results are compared with the measured Q values for the water injection test, as shown in Figure 21. The calculated values of Kxx and Q are compared with the measured values of water injection, as shown in Figure 22. From Figures 20–22, it shows that the theoretical calculation is in good agreement with the experimental results, which proves the validity of the theoretical model.Table 3
List of theoretical models and water injection test.
Df
Kxx calculated value (cm/s)
Kxx test value (cm/s)
Kxx calculated value/Kxx test value
Q calculated value (L/min)
Q test value (L/min)
Q calculated value/Q test value
1.266
0.16
0.15
107.26%
13.16
13.68
96.19%
1.324
0.36
0.36
101.06%
19.20
21.46
89.48%
1.306
0.10
0.10
98.36%
21.90
20.00
109.51%
1.313
0.15
0.16
96.41%
32.37
29.53
109.61%
1.401
1.26
1.18
106.76%
175.53
157.95
111.13%
1.451
2.31
2.22
103.95%
322.91
298.43
108.21%Figure 20
Comparison of theoretical model and water injection testKxx.Figure 21
Comparison of theoretical model and water injection testQ.Figure 22
Theoretical calculated value/water injection test measured histogram ofKxx and Q.


## 6. Conclusions

In this paper, a fractal model for characterizing hydraulic properties of fractured rock mass under mining influence considering fracture surface roughness is established, and the infiltration tensor of fracture network is deduced under two kinds of seepage conditions. The effects of fractal geometry and fractal characteristics on the permeability are analyzed. The validity and accuracy of the model are verified by comparing the theoretical model with the physical model water injection test. The results show the following:(1)
The permeability coefficientK of the fracture network is a function of the geometric of fracture network. In general, the maximum fracture length Lmax, calculation section porosity Φ and fracture network fractal dimension Df reflect the related indicators of the calculation section, such as the fracture length, discharge area, and density, which play a role in promoting permeability. The horizontal projection length L0 and the fractal dimension DT of the seepage flow line reflect the obstacle of the surface roughness and the increase of the seepage path to the seepage fluid, which inhibit the permeability
(2)
With the increase of fractal dimensionDf, the permeability coefficient increases. While with the increase of DT, the permeability coefficient decreases rapidly, and the larger the Df (Df>1.5), the greater the change of permeability coefficient K with DT, which indicates that when Df>1.5, the weakening effect of fracture surface roughness on the permeability of fracture network increases rapidly
(3)
Df and DT have a greater impact on the permeability of the fracture network, but their physical meanings and influences are different. The larger the Df value, the more complex the fractured network in the cross section, and the higher the fracture density, thus the higher the permeability coefficient. DT reflects the rupture degree of the fluid seepage path in the fracture due to the rough surface of the fracture. With the increase of DT, the seepage path is lengthened and the effective overcurrent capacity of the fracture is reduced, which leads to the decrease of the permeability coefficient of the fracture network
(4)
With the increase of porosityΦ, the permeability coefficient increases. While with the increase of DT, the permeability coefficient decreases. The higher the DT (DT>1.2), the smaller the change of permeability coefficient K with the change of Φ, which shows that when DT>1.2, the fracture surface roughness will be greatly offset due to the increase of water seepage area leading to increasing effect on permeability coefficient
(5)
With the increase ofLmax, the permeability coefficient increases. Simultaneously, with the increase of Df, the overall permeability coefficient will increase, and with the increase of Df (Df>1.5), the rate of permeability coefficient K with Lmax change sharply increases, which shows that when Df>1.5, with the increase of fractal dimension, the influence of long fractures on fracture network seepage is increasing



---

*Source: 8391803-2019-12-20.xml*