A Complex Network Model for Analysis of Fractured Rock Permeability

State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou, Jiangsu 221116, China Mechanics and Civil Engineering Institute, China University of Mining and Technology, Xuzhou, Jiangsu 221116, China Fractured Coal Masses Laboratory of Mine Cooling and Coal-Heat Integrated Exploitation, China University of Mining and Technology, Xuzhou, Jiangsu 221116, China


Introduction
e networks of rock fracture are formed by structural deformation and physical diagenesis [1]. On the inner surface of rock stratum, the scale of naturally formed fracture network is expanding and redistributing randomly with different fracture development degrees, which is always difficult to identify. However, researchers usually use the dip angle and azimuth to determine the spatial orientation of fractures. e structural characteristics are consistent with the two directional attributes of structural geology: tendency and trend. In geological objects, the complex trace analysis is used to calculate the obliquity estimation of the three-dimensional data body, i.e., steering cube [2], and then gets the dip angle and azimuth information of each data point. e permeability of fractured reservoirs is very low, and the fracture network controls the fluid flow [3]. Hence, it has an important influence on oil or gas exploitation [4] and geothermal energy extraction [5].
In recent years, researchers around the world have studied permeability characteristics of fracture networks and put forward corresponding models [6][7][8]. Snow [9] in his study established the parallel plate model and obtained a tensor analytical formula of permeability of fracture network. Koudina et al. [10] studied the permeability of fracture network in three-dimensional space by means of numerical simulation. e fracture network was composed of polygonal shapes and the flow of fluid in each fracture satisfied Darcy's law, while comparing it with Snow's model. Xia [11] established the dynamic model of permeability and opening of fracture network under different confining pressures. Van Stappen et al. [12] also connected the seepage model with fracture opening by determining the relationship between fracture permeability and confining pressure. Li et al. [13,14] broke away from the traditional practice of thinking fractured reservoirs as dual media and established a percolation model with equivalent continuum suitable for low-permeability fractured shale reservoirs. De Dreuzy et al. [15] studied the permeability of randomly generated two-dimensional fracture network by numerical and theoretical methods and compared it with natural fractures to verify the accuracy of the model. Klimczak et al. [16] used the parallel plate model to obtain the permeability formula of a single crack under conditions that the fracture length and opening satisfied the power-law relationship and verified the accuracy of the model through numerical simulations. Wei et al. [17] derived a forecasting model of permeability using the electrokinetic relationship between fluid flow and current in microfractures and analyzed the influence of connectivity between fractures on permeability. Li [18] proposed a new model considering fracture connectivity according to the hydraulic fracture morphology of raw coal, the "matchstick" seepage model and the cubic law. However, the above models do not quantitatively relate the permeability of fracture network with porosity, surface density, and microstructural parameters of fracture network, such as fracture connectivity, openness, the dip angle, and azimuth. e randomly distributed fracture networks in rocks have been shown to have statistical self-similarity, which is a basic feature of fractal. Interested readers may consult [19][20][21][22][23][24][25][26] for details. Watanabe and Takahashi [5] used the fractal theory to study the permeability of fracture networks and the extraction of heat in dry hot rocks, but they did not put forward a permeability expression with micro parameters. Yu et al. [27] based on the study of seepage characteristics of porous media in fracture networks by using fractal methods put forward an explicit expression with micro parameters, such as the structure of fracture network and porosity, and then gave the scaling relationship between permeability and the structure of fracture network. Li et al. [28] established the mathematical model of equivalent permeability tensors in fractured reservoirs, based on fracture statistics, the simulation technique of fracture network, and equivalent flow assumption, and then obtained the equivalent permeability tensor of fractured media by using boundary element method. Jafari and Babadagli [29] obtained fractal permeability expressions of random fractures by using multiple regression analysis based on logging data but their empirical relationship contained many empirical constants. Recently, Miao et al. [6] obtained the analytical expression of fracture network permeability according to the basic fractal theory. is model quantitatively connected the fracture length, aperture, the fracture dip angle, and fracture azimuth with permeability of fractured rocks, which did not include any empirical constant.
Most of the above models initiate from the statistical parameters of isolated fractures and macroscopic homogenization. e connectivity of fracture networks is not considered, particularly the influence of the connectivity of a small number of local fractures (maximum degree) on the overall permeability. Starting from the topological structure of fracture network and based on the complex network theory, this paper establishes the network permeability model of fractured rocks and probes into the internal mechanism of the influence of structure parameters of fracture network on permeability, including fracture porosity ∅ M , fracture density D, power index d k , and the maximum node degree k max .

Degree Distribution of Hierarchical Networks.
In order to illustrate the modularity, local clustering, and scale-free topological characteristics of many complex network systems, it is necessary to assume that the modules generate a hierarchical network in some iterative way [30]. Recent studies show that [31] some of the topology modules are well organized hierarchically in the network. e hierarchical network seems to have a very conspicuous feature; that is, the local is similar to the whole in a sense, i.e., the selfsimilarity. Hierarchical network integrates scale-free topology with internal module structure. Song et al. [32] further reveal that self-similarity and degree distribution of scale-free hold true at all coarse-grained stages of the network by adopting renormalization procedure and the degree distribution P(k) of the renormalized network is invariant under renormalization. e power-law relationship can be expressed as follows [32]: where P represents the total number of nodes in the network with degree k, k represents the number of other nodes connected to a node, and c represents the self-similarity index with the range of 1-3, which is transformed by the exponential formula [32].

Basic Features of Fractals.
Most of the trace length of the fracture satisfies the power-law (scale-free) distribution [33,34]. e fractal power-law distribution refers to the fact that the fracture length in nature is random and disordered, showing the characteristics of similarity and fractal. e power-law expression is [35] where D f is the fractal dimension of fracture length, l is the track length of fracture, and N is the total number of fractures. is is the basic expression of fractal scaling law and basis of box counting method.

Power-Law Expression of Fracture Complex Networks.
Covariant analogy is also known as mathematical similarity analogy. Power-law relations (1) and (2) have obviously similar functional relations and equation (1) multiplied by k can be analogous to equation (2). e number of edges of complex network is associated with the number of edges of fracture network. e following expression is obtained: e relationship between c in power-law distribution formula (1) and D f in power-law expression (2) is as follows [32]: where the power exponent d k is 1.5 times of D f [32]. Substitute equation (4) into equation (3) to get a proportional relationship: e parallel plate model is usually used to represent the effective aperture of the fracture and the relationship between crack length and effective aperture has also been studied by a large number of researchers [36,37]. is relationship is given by where β is the proportionality coefficient, which is related to the mechanical properties of the medium around the fracture in the range of 10 − 3 ∼ 10 −1 [16]. a is the effective aperture of the fracture and n is the power exponent. When the power exponent n � 1, the fracture network has the characteristics of self-similarity and fractal [37]. So, for the fracture network with self-similarity [16], equation (6) can be rewritten as Equation (1) can be rewritten as where M represents the number of network edges, and α is the proportionality coefficient. Differentiating equation (8), we can get the number of edges whose node degrees are in the range k to k + dk: wherein the negative sign indicates that the number of edges of a complex network decreases with the increase of node degree, which is in line with the actual situation and −dM(k) > 0. e probability density of an edge with node degrees k is expressed as where M t represents the total number of edges the network has, and f(k) � (α/M t )ck −c is the probability density function of the edge with node degrees k, which satisfies the normalization principle: us, it can be obtained that Evidently, when k min ≪ k max , equation (12) can be expressed as generally, k min /k max ≤ 10 −2 can be taken and complex networks in nature usually meet this requirement. Yu [38,39] studied the power-law relation of fractal distribution of pores in porous media. Likewise, Majumdar and Bhushan gave the cumulative size distribution of islands on the Earth's surface [40]: where N is the total number of islands with the largest area s max greater than s, and D is the fractal dimension for the size distribution of islands. Equation (14a) indicates that there is the largest island on the Earth's surface; in addition, Majumdar and Bhushan [40] used this power-law formula to describe the contact points on engineering surfaces, where the maximum point area s max � gλ 2 max , a point area s � gλ 2 with λ being the diameter of a point and g being a geometry coefficient.
Since self-similarity is one of the basic characteristics of fractal, the self-similarity of porous media with fractures needs to satisfy a certain power-law relationship [41]. Hence, equation (14a) is used to describe islands on the Earth's surface and points on the engineering surface can be extended to describe the size distribution of nodes on the surface of a fracture network. In the complex network theory, the characteristic size of a single node includes outdegree and in-degree [30]: where k omax k imax represents the maximum node size with k omax and k imax , respectively, being the maximum out-degree and maximum in-degree, k o k i is a node size with the outdegree and in-degree being k o and k i , respectively. When the direction of the degree is ignored, equation (14b) can be simplified as From equation (14c), the cumulative number of nodes whose degrees are greater than k can be expressed as Advances in Civil Engineering 3 where N is the cumulative number of nodes in a fracture network. From equation (14d), the total number of nodes in a complex fracture network is obtained: Because the contribution of one edge to degree is 2, the average degree of complex networks is Inserting equation (15) into equation (16), Inserting equation (17) into equation (13) to get the proportionality coefficient, en, we insert equation (18) into equation (9) which gives Equation (19) is an important power-law distribution relation of edges with certain node degree in complex networks. Furthermore, by the same logic, the average degree of complex networks can be obtained as 2.4. Surface Porosity of Fracture Networks. Self-similarity is closely related to the fractal. Yu and Li [42] deduced the relationship between porosity and fractal dimension in porous media based on the fractal theory: where λ min and λ max are, respectively, the minimum pore diameter and the maximum pore diameter. D p is the fractal dimension of pores. d E is the Euclidean dimension: in two dimensions, d E � 2; in three dimensions, d E � 3. Equation (21) is appropriate not only for precise fractal geometry but also for statistical fractal geometry. As long as the pores of porous media fall within the self-similar range of λ min ∼ λ max , forming a fractal set, equation (21) holds accurate regardless of the shape of the pores. erefore, in a hierarchical complex network, it is embedded into the matrix as a fracture, forming a network model with fracture properties. e edges of the complex network, that is, the fractures with a node, satisfy the above equation (21) within the self-similar range of k min ∼ k max and are independent of the shape of the node. It can be rewritten as where ∅ M is the effective porosity of fractures in the rock and k min and k max are the minimum and maximum of the nodes, respectively. On the cross section of the representative elementary volume, the surface porosity of the fracture network is defined as [6] where A M represents the cross-sectional area of the representative elementary volume in which the fracture network is located, and A PM represents the total area of fracture pores on this area. According to equations (5), (7), and (19), we can get the total cross section area of the fracture [6]: Inserting equation (22) into equation (24), where the porosity ∅ M is used in the two-dimensional space of equation (22), i.e., d E � 2.

Relationship between Surface Density and the Self-Similarity Index
According to equation (5), the total length of fractures on the cross section of the representative elementary volume is as follows: Inserting equation (22) into equation (26),

Advances in Civil Engineering
In the two-dimensional fracture network, the surface density refers to the density of the cross section in a unit cell of fractures (not a single fracture), which is defined by [43] where D is the surface density of fractures, and L total is the total length of all fissures on the cross section of the representative elementary volume body, which is related to the complex network model. Equations (23), (25), and (27) are inserted into equation (28) to get the surface density of the fracture, i.e., Equation (29) shows that the surface density in the twodimensional complex fracture network is a function of the fracture porosity ∅ M , self-similarity index c, power index d k related to fractional dimension, proportion coefficient β, and the degree k max of the largest node.
In order to study the relationship between the surface density of these four fracture networks and the self-similarity index, the prediction results of the surface density for the complex fracture network are compared with the four random fracture networks generated by Zhang and Sanderson [43] through the numerical method of selfavoiding walking. In their simulation, the critical fractal dimensions lie in a narrow range from 1.22 to 1.38 (average 1.30) for those critical clusters with variations in the lower limit of length from 0.005 to 1.5 m, in the dispersion angle of fracture direction from 0 to 50°, and in exponents from 1.2 to 1.8. erefore, during calculation, the minimum degree of the node is taken to be 1 and the maximum degree is taken to be 300 according to the same ratio coefficient. Meanwhile, the average power index is 2 and the average self-similarity index is 1.65 through equation (4), and the average porosity ∅ M is calculated through equation (22). From Figure 1, it can be observed that the predicted results are in good agreement with numerical simulations. Meanwhile, Figure 1 shows that the surface density of fracture network increases with the increase of fracture self-similarity index. Figure 2 shows the relationship between fracture surface density and porosity when the maximum degree of node k max is 300 and β is 0.006. It can be observed from Figure 2 that the surface density of fracture network increases with the increase of fracture porosity. is is because the larger is the porosity, the greater will be the pore area of the fracture network. Under certain conditions of β, the longer is the total length of the fracture, the stronger the connectivity will be as mentioned. is result is consistent with the simulations by Miao et al. [6]. It can be explained that the change of numerical values will not affect the general trend between them.

Complex Network Model of Permeability of Fractured Rocks
Generally, the production of low-permeability reservoirs often depends on the seepage system of fracture network. When there are differences of temperature and pressure in the system, there will be fluid flow or heat transfer between the fracture networks. In these processes, however, the laws of mass, momentum, and energy transfer among fluids are very complicated. Moreover, the geometric aspects of fractures cannot be determined relatively, including density and surface roughness. For fractured reservoirs, we can proceed from macroheterogeneity, because the degree distribution, length, aperture, and orientation of the fracture are often random and disordered. Complex network can provide an effective method for representing irregular objects. e topological model of complex networks normally considers the positional relationship between nodes but not their shape and size. erefore, network space determined by the azimuth and dip angle has an important influence on the seepage characteristics of fracture network. Nevertheless, the spatial orientation of fractures is usually random and the number of fractures in space is so large that it is almost impossible to express the orientation of each fracture precisely [44]. Generally, the statistical method in the field of engineering has been adopted to show the location of fracture network, which is to take the average value of fracture dip angle and fracture azimuth [45], and this is often used in petroleum engineering, shale gas exploitation, and geothermal energy extraction. erefore, in this paper, we assume that the average dip angle of the complex fracture network is θ and the average azimuth of the fracture is α, as shown in Figure 3. e cubic law of single fracture is based on the model of parallel plate, which becomes the basic theory of network seepage of fractured rocks and this is usually considered simple and effective. e flow rate along the flow direction through a fracture can be described by the famous cubic law [46,47]: where L 0 denotes the length of the representative elementary volume, a denotes the fracture aperture, l denotes fracture trace length, Δp denotes the pressure drop across a fracture along flow direction, and μ denotes dynamic viscosity coefficient of the fluid. If the spatial orientation of fracture is considered, the flow rate of a single fracture can be expressed as [6,48] e total flow rate of fluid through a set of complex fracture networks can be obtained by integrating equation (31) from minimum degree to maximum degree in a unit cross section; i.e.,

Advances in Civil Engineering
In general, k min ≪ k max . According to equation (4) and [35], since 1 < c < 2.3 in the two-dimensional plane, and (k min /k max ) 4/d k (− c+1) ≪ 1, consequently, equation (32) can be simplified as It can be seen from equation (33) that the total flow rate of fluid in the complex fracture network is related to the index of self-similarity c, fractional-dimension related power index d k , fracture azimuth α, and fracture dip angle θ and the flow rate is very sensitive to the maximum degree k max of nodes.
Darcy's law for Newtonian fluid flow in porous media is given by [6] e permeability of complex fracture network can be obtained by inserting equation (33) into equation (34):

Fracture
Flow direction x y Horizontal plane α θ Figure 3: e average orientation of fractures in the three-dimensional space, the plane of the coordinate axis is the horizontal plane and the direction of water flow is along the x-axis. e included angle between the fracture direction and y-axis is α, that is, the azimuth of the fracture. e θ angle between the fracture plane and the horizontal plane is the dip angle of the fracture.

Advances in Civil Engineering
By inserting equation (28) and equation (29) into equation (35), the permeability of fracture network can be expressed by the surface density of fracture: Equation (36) suggests that permeability is a function of self-similarity index c, power index d k , structural parameters (maximum degree of a node k max , fracture surface density D, fracture azimuth α, and fracture dip angle θ), and fracture porosity ∅ M in a medium formed by a complex fracture network. Equation (36) further reveals that permeability is strongly dependent on the maximum degree k max of the node. e higher is the node degree, the stronger is the connectivity of fracture network. e fluid capacity increases with increase of the flow path, leading to higher permeability. erefore, this model has more advantages than the traditional model and can better explain the influence of node failure on fluid flow in the fracture network.

Results and Discussion
Jafari and Babadagli [49] analyzed 22 different fracture networks in nature. e digitized fracture patterns were exported to commercial fracture modeling software (FRACA) to calculate their equivalent fracture network permeability. A 3D model with a grid block size of 100 m × 100 m × 10 m was constructed. Each digitized 2D fracture pattern (i.e., the digitized mapped fracture traces from outcrops) was imported into the 3D model in such a way that all fractures were considered to be vertically touching the top and the bottom of the layer, wherein the maximum fracture length is 2 m and the dip angle of the fracture is 0°. erefore, in the calculation, the minimum degree of the node is 1, and the maximum degree is 6. Furthermore, since the model of parallel plate mainly depends on the effective aperture of a single fracture, the actual tortuosity of the fracture is not considered by using this simplified model. Via equation (4) and equation (20), the average power index and the average degree of the node are calculated. All the structural parameters used in theoretical calculations are listed in Table 1. Figure 4 shows that the predicted values of our model are in good agreement with the results of numerical simulations.
We discuss the influence of model parameters on permeability. From equation (36), it is observed that the parameters that play a decisive role mainly include fracture porosity ∅ M , fracture dip angle θ, fracture surface density D, power index d k , and maximum node degree k max . Figure 5 shows the relationship between permeability and fracture porosity of the complex network model at different dip angles. In the calculation, the maximum degree k max � 398 of the fracture node (at β � 0.006) is taken. It can be seen from Figure 5 that the permeability of fracture network increases with the increase of fracture porosity. In addition, with the same porosity, the larger is the fracture dip angle, the smaller is the permeability of fracture network. is is because the flow resistance of fluid increases with the increase of the fracture dip angle. Figure 6 shows the relationship between permeability and fracture surface density in the complex network model. In this calculation, the maximum degree k max � 398 of fracture node, fracture dip angle θ � 45°, fracture azimuth α � 0°, and β � 0.006 are taken. It can be seen from Figure 6 that the permeability of fracture network increases with the increase of fracture surface density. is is because as the density of fracture surface increases, the porosity of fracture also increases. Hence, the permeability of fracture network increases. Figure 7 shows the relationship between permeability of complex network model and the power index. In the calculation, the minimum and maximum degrees of fracture   Advances in Civil Engineering 7 network nodes as 1 and 6 are taken, respectively. Equations (4) and (20) are used to calculate the average self-similarity index c av � 1.67 and the average degree of node k av � 4 and we take the dip angle of fracture θ � 0°and β � 0.001. It can be seen from Figure 7 that the permeability of fracture network decreases slowly with the increase of power index. Miao et al. [6] have verified that the permeability of the fractal fracture network model increases slowly with the increase of fractal dimension. By considering equation (4), that is, the internal correlation between scale-free property of complex networks and fractal scaling law, it can be concluded that there will be a competitive relationship between the inhibition of seepage flow by power index and the promotion of seepage flow by fractal dimension. Henceforth, it leads to the discontinuous phase that is not always occurring. Figure 8 shows the relationship between the permeability of complex network model and the maximum node degree. When a node in a network has multiple edges connected to it, the number of edges is the degree of the node, regardless of its direction. In the calculation, we take the dip angle (from 0°to 180°) of fracture θ � 45°, azimuth of the fracture α � 0°, β � 0.006 with the range of 0.001∼0.1 [16], and average surface density D � 10 (m/m 2 ). It can be seen from Figure 8 that permeability of fracture network increases sharply with the increase of the maximum degree of nodes. Since connectivity of the entire fracture network is strongly dependent on maximum degree of a node, it is equivalent to the connection hub of entire complex network. When a small number of edges are removed from the network, the overall connectivity of the network will not be greatly affected. us, the complex network has a high robustness to the node destruction. At the same time, if a node with the maximum degree is deliberately attacked, the entire network will paralyze quickly and the fluid can only flow through a few paths. is is also the vulnerability of complex network to deliberate attacks on nodes.

Conclusion
is paper applies the complex network theory and topological model to fractured rocks, while describing the   Figure 8: e relationship between fracture permeability and maximum node degree. 8 Advances in Civil Engineering fracture network as a hierarchical network with self-similarity. Meanwhile, the fracture network model of surface density is obtained based on the power-law distribution relation of network edges. en, the permeability model of fractured rocks is deduced in accordance with the famous cubic law, Darcy's law, and complex network theory. Compared with the existing numerical simulations, the predicted results show that the above models are accurate. Besides, the effect of structural parameters on the permeability of fractured media is also discussed. e permeability of fracture networks increases with the increases of porosity and surface density. e permeability of fracture networks increases exponentially with the increase of the maximum node degree and its power exponent is 3/d k .
Data Availability e data (numerical simulation) used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding the publication of this paper.