Connectivity evaluation for three‐dimensional fracture network in support‐based model: A case study in the Ordos Basin, China

Connectivity is an effective indicator to describe the flow paths and mechanical behavior of reservoir fractured rock masses. The connectivity evaluation of fracture networks plays an important role in the design, assessment, and development of reservoirs in several important engineering applications, such as geothermal recovery and petroleum industry. Based on Allard's definition of a connectivity index, a 3D connectivity evaluation method is proposed in this study by building the connectivity functions of the supports. Meanwhile, different connectivity functions are built for different situations. In addition, a convex polygon is used to simulate the fracture to make the fracture network as real as possible. The accuracy of the generated fracture network was verified by comparing three parameters of fracture, and the results showed that the relative error was less than 5%. Some examples in the Ordos Basin, China, are given to verify the correctness of the program and method. Analysis results show that the 3D method is more realistic than the 2D connectivity evaluation, and the connectivity field of the fracture network is more reliable in the 3D method. Moreover, the proposed method can show the potential preferential flow path in a stereoscopic manner and thus greatly extends the function of connectivity analysis. This method provides important information for the prediction of paths and the evaluation of the connectivity of the fracture network and is important in evaluating the optimal location of geothermal recovery and petroleum industry.


| INTRODUCTION
Shale gas is an important clean and efficient energy source that is widely distributed in the Ordos Basin of China. [1][2][3] The permeability of shale gas reservoir matrix is low, and the production process requires fracturing of the reservoir to generate economic benefits. [4][5][6] In general, the permeability of the reservoir and the fracturing pattern of rock mass are controlled by natural fracture. The natural fractures and rock blocks are connected to each other and control the flow characteristics of the reservoir. In this case, natural fractures have a significant impact on shale gas production capacity and can greatly improve production efficiency. 7,8 Therefore, analyzing the fracture system of shale gas reservoirs is important.
Many researchers have used the discrete fracture network (DFN) model to simulate the fracture system of reservoirs. 9,10 This model can intuitively show the distribution of fractures in underground space. On the basis of the DFN model, the evaluation of the connectivity can be used as an effective method to describe the flow path and mechanical behavior of underground rock mass. [11][12][13] In this model, the hydraulic fracturing and oil and gas development are based on the fracture network. In addition, the connectivity of the DFN will directly affect the permeability of the model. 14 By analyzing the connectivity of fracture networks, optimal drilling locations can be evaluated to improve the efficiency of production, and the potential flow path in rock mass can be predicted. Therefore, evaluating the connectivity of fracture networks is essential to predict the flow behavior of the fracture system in reservoirs.
Most studies are based on 2D space. For example, Knudby and Carrera 15 presented nine indicators of connectivity, which are based on either statistical characteristics or flow and transport parameters. They analyzed these indicators with respect to their ability to explain the difference in flow and transport for sets of fields with low and high connectivity. Hestir 16 and Sahimi 17 used the average number of intersections per fracture, x f , to evaluate the connectivity of fracture network. The results showed that the larger the x f value is, the higher the permeability of the model is. Li et al 18 defined connectivity on the basis of the statistical geometrical parameters of a network of random discontinuities to quantify the hydraulic properties of the network using connectivity. They proposed an analytical method to evaluate connectivity by using the density of degree of freedom (DDOF). The results showed that the connectivity improves and the permeability increases with the increase in DDOF. Tyukhova et al 19 obtained a connected channel network by using information about hydraulic conductivity only and defined static connectivity metrics on the basis of the resistance and width of the individual channels. However, these results are unrealistic because the actual engineering is based on 3D space. Therefore, some researchers have proposed few methods to evaluate the connectivity of fracture networks based on 3D space. Outters et al 20 used the fracture areal intensity parameter P 32 to characterize the connectivity of a 3D model. The result showed that a large value of P 32 indicates high connectivity. Li et al 21 calculated various connectivity parameters, such as the number of intersection lines and the length of intersection lines between fracture surfaces in volume, to analyze the dominant flow direction in 3D fracture networks. Xiong et al 22 established a 3D DFN model for nonlinear flows. The results showed that the hydraulic gradient for the onset of nonlinear flows is low when the fracture surfaces are smooth and the connectivity of the fracture network is high.
In summary, these evaluation methods are based on either 2D space or geometric parameters of fractures, such as the density, size, and direction of fractures. However, the locational relationship between fractures is ignored, and the connecting information between fractures is missing. Facture connections play a decisive role in the connectivity of the DFN model. 23,24 In general, fracture channels are more important than fracture length and density in evaluating the connectivity of fracture networks. Fracture networks with three different geometric characteristics in the 2D plane are shown in Figure 1 for illustration. The fractures shown in Figure 1A are long and have small gaps between them. However, the permeability of the fracture network maybe very low because of the absence of fracture pathways throughout the study area in horizontal and vertical directions. As shown in Figure 1B, the density of fractures is very high. No fracture pathways are observed because of the short length of the fractures. Thus, the connectivity of the fracture network may be low. Figure 1C shows a low-connectivity network with the large average number of intersections per fracture. These examples show that the fracture connections in the fracture network model can directly affect the connectivity. Therefore, the connection between fractures should be considered when evaluating the connectivity of fracture networks.
On the basis of Allard's definition of connectivity index, 25 a connectivity evaluation method a 3D fractured network is proposed to completely preserve the interconnection information between fractures and describe the flow path of fractured rock mass in 3D space. An example is analyzed by combining with the outcropping fracture data around Tongchuan City in the northern region of the Ordos Basin to provide guidance for shale gas exploration and development in the area.

NET WORK
Fractures are widely distributed in the interior of the rock mass and play a major control role in the physical behavior of the rock mass. In the 3D space, the rock mass is not transparent. This condition limits the information about the real fracture distribution and geometric characteristics. Most studies have shown that fractures are distributed uniformly, logarithmically, and normally in 3D space. 26 However, the distribution must be represented by certain geometric features, such as the size, shape, orientation, and location of the fractures, regardless of its type. Moreover, the density of fractures and other factors should be considered in regions where fractures are sufficiently developed. This study constructs the fracture network mathematical model using MATLAB software.

| Fracture shape
In the 3D space, a fracture is presented as a surface, which is different from that in 2D space where it is considered a line. Due to the fractures have a certain shape, many studies have adopted the disk model to simulate fractures. 27,28 However, the disk fractures cannot simulate rock anisotropy effectively, as shown in Figure 2. Figure 2A shows that a disk fracture is present in the rock mass, and the mechanical properties and permeability are consistent in the horizontal and vertical directions. However, when the shape of the fracture is polygonal, the anisotropy of the fractured rock mass can be displayed evidently. In fact, the disk is the limit case of an n-sided shape, that is, the fracture is a disk when n approaches infinity. Shi et al 28 proved that disk and polygonal fractures can be used to describe the geometric characteristics of the structural plane. Polygonal fracture can simulate the heterogeneity and anisotropy of fractured rock mass effectively, as shown in the dotted line in Figure 2B. Thus, polygons can be utilized instead of disks in describing fractures.
Polygons can be divided into a series of triangles, and the efficiency of numerical calculation is optimized. Some research has shown that the triangular fracture plane has poor ductility, whereas the quadrilateral one has good ductility and vertices that are rarely second only to those of the triangular one. Thus, the latter can simulate the topological relations of fracture networks with simple data effectively. Therefore, this study uses quadrangles to simulate fractures of rock mass in 3D space. Figure 3 shows the algorithm to achieve quadrilateral fractures. The z-coordinate is ignored to show how to construct quadrilateral fractures. First, two arbitrary points (a, b) and (c, d) in space are taken to construct rectangle e 1 e 2 e 3 e 4 , as shown in Figure 3. Four points A, B, C, and D are randomly selected on the four sides of the rectangle as the four vertices of the quadrilateral fracture. rand is a random number between 0 and 1. The convex quadrilateral fractures of any shape can be constructed by transforming the coordinates of the four vertices. This method can avoid the generation of concave quadrilateral. The calculation of the centroid and area of concave quadrilateral in 3D space is complex. The intersection line between two concave quadrilaterals in 3D space is also difficult to determine, as shown on the left of Figure 4. The aforementioned parameters are important in studying the geometrical topological characteristics of fracture networks. When the centroid is located on the outside, the area of the fracture can be calculated as S ⌷1234 = S ∆125 + S ∆235 + S ∆541 − S ∆534 . When the shape of the fracture is convex quadrilateral, the area of the fracture is S ⌷1234 = S ∆ 125 + S ∆235 + S ∆345 + S ∆541 . However, the symbols for the concave quadrilateral should be determined. The symbols in the formula are consistent, and this way greatly reduces the complexity of calculating the fracture area and improves the efficiency of the numerical calculation. Figure 4 shows the intersection line of two fractures with different shapes. The intersection lines between concave quadrilateral fractures are two line segments T 1 S 1 and T 2 S 2 , and that for the convex quadrilateral fracture is line segment T 1 T 2 . In terms of modeling calculation, the intersection of concave quadrilateral fractures is more complicated than that of convex ones. Therefore, convex quadrilateral is adopted to simulate fractures of rock mass in this study.

| Fracture location and size
The location of any quadrilateral fracture in 3D space can be represented by quadrilateral centroid, which is mainly determined by four vertex coordinates. Figure 5 shows that the centroid coordinate (x c , y c , z c ) can be expressed as follows according to Helen's formula: where (x c1 , y c1 , z c1 ) and (x c2 , y c2 , z c2 ) are the centroid coordinates of ∆ABC and ∆BCD, respectively; and S 1 and S 2 are the areas of ΔABC and ΔBCD, respectively.
where a, b, c, d, and e are the side lengths of ∆ABC and ∆BCD p 1 and p 2 are half the circumference of ∆ABC and ∆BCD, respectively. The size of the fracture can be expressed by the area of the quadrilateral, that is, S = S 1 + S 2 .

| Fracture orientation
In 3D space, the normal vector of the fracture is obtained from the cross product of two adjacent edge vectors to describe its direction when the coordinates of four vertices of the quadrilateral fracture are known. As shown in Figure 5, the AB vector is (x 2 − x 1 , y 2 − y 1 , z 2 − z 1 ), and the AB vector is (x 4 − x 1 , y 4 − y 1 , z 4 − z 1 ). Then, the normal vector of fracture surface can be obtained by

| Generation of random fracture network
When the z-axis coordinates are considered, convex quadrilateral fractures at any angle can be constructed in 3D space by rotating the fracture surface. Figure 6A shows the spatial transformation process of convex quadrilateral fractures in 3D space. Figure 6B,C shows the randomly generated 3D fracture network. The coordinates of the four vertices of the rectangle in 3D space, e 1 , e 2 , e 3 , and e 4 , can be expressed as The fracture surface is rotated around the centroid by establishing the transformation matrix corresponding to the rotation angle. At this point, the normal vector of fracture surface can be described by the inclination angle, inclination, strike, and other parameters. 27 The relationship between the normal vector of fracture, dip angle α, and dip direction β can be expressed as

| Distribution of fracture geometric parameters
In general, each geometric parameter is assumed to obey a certain distribution function in 3D DFNs. 29,30 For example, the location (the centroid) of fractures in 3D space obeys Poisson distribution. Distribution functions, such as Bingham distribution, Arnold hemispheric normal distribution, Fisher distribution, and bivariate normal distribution, are commonly assumed for fracture orientations. In particular, Fisher distribution has a good fit for fracture orientations by analyzing a large number of outcrop fracture data. The existing technology cannot measure the exact area of the fracture surface from the field but can obtain the measurement data of fracture trace length. If the fracture trace length is assumed to be circular or quadrilateral, then the area of the fracture can be estimated using the crack trace length. The relationship between the two can be established, as shown in the following expression: where S is the average fracture area and S′ is the fracture area calculated in accordance with the length of the fracture trace. Figure 7 shows the flowchart of the 3D fracture network modeling.

| Principles of connectivity
On the basis of the definition of connectivity index, 25 we assume that the connectivity of two points only depends on whether a fracture connects two points in the fracture network. Figure 4 shows two different connectivity measuring methods. The lattice (cell)-based connectivity of two points is defined only by the relative locations of the grids in Figure 8 (left) (eg, overlapping and adjoining). In other words, the information of intersection between fractures cannot be preserved, and it mainly includes two types of 2D connections (eg, 4 or 8 connectivity scenarios for lattices) and three types of 3D connections (eg, 4 or 8 or 26 connectivity scenarios for lattices). Figure 8 (right) shows that the connectivity measuring methods we assumed are defined only by the connections between fractures in support (a support is any defined subspace of the study region varying in size from a point to the entire region)-based methods. These methods can retain any interconnection information between fractures and involve two different sampling methods, namely systematic and random sampling. The figure shows that the result is the same regardless of the sampling method. In random sampling, even if the two sample supports are not adjacent, they can be considered connected because of the presence of a pathway between them via fractures. However, when no pathways are present between two sample supports, they are considered not connected even if they overlap each other. Similarly, the connection between supports is independent of the relative locations of (5) n = ( sin sin , sin cos , cos ) Centroid location of fracture the supports in system sampling. In other words, the connectivity of any two supports depends only on the connections between fractures and not on the relative locations of the two supports. Figure 8 shows fundamental differences between the two methods. Supports are not physical entities and require some specific grids to describe them. If a grid is used to describe the support, then it is not considered a lattice because of the definition of the support. The supportbased method is used in this study because it partially solves the uncertainty errors caused by lack of data or data errors and can evaluate the connectivity of the research area effectively. Figure 9 shows three different spatial positional relations between fracture surfaces in 3D space. Figure 10 shows the connections between fractures in 3D space.
F I G U R E 6 3D fracture network F I G U R E 7 Flowchart of 3D fracture network modeling using fracture data from outcrops

| Connectivity function
The connectivity function C(x, y) is based on the definition of Allard. 31 Two supports x and y (eg, a volume in 3D or an area in 2D) in the region R. If they are connected, then they are denoted by x ↔ y; that is, a fracture pathway exists between the supports. At this time, the connectivity function value is 1; otherwise, it is 0, as shown below: The function quantifies the connectivity between any two supports in a space based on the connections between fractures. The connectivity function can determine whether where N is the number of MC simulations. The location of support x 0 remains unchanged, and the location of support y 0 is changed to ensure it covers the entire study region by using the system sampling scheme. As a result, the connectivity function is expressed as where j is the number of supports satisfying y j ∈ R. The global connectivity (GC) (ie, covering the entire study region) of support x 0 can be evaluated by constantly changing the location of support y 0 in Equation (9). The connectivity value for each support y 0 is then calculated. The right side of the expression is the estimated formula for the discrete case. When the supports are located through systematic sampling on a coarse 3D grid (3 × 3 × 3 supports) for the specific fracture network, we assume that each support y j has a number; that is, y j = y αβγ . α, β, and γ are the location markers of the subspace in the research area, as shown in Figure 11. In summary, the procedure for evaluating the GC in three dimensions begins by finding the extent (the research area R) of the fracture network, which is then divided into cubic shapes (ie, the 3D support v) in the resolution (the side length of cubic shapes satisfies the calculation requirements). A 2D example is calculated, as shown in Figure 12, to better describe the function of Equation (9). We assume that the scale of the research area is 10 m × 10 m, and fracture locations follow random Poisson distribution. Using the system sampling, the entire 2D research area can be arranged into 100 supports. Supports x 0 are positioned at x 55 , x 56 , x 57 , x 65 , x 66 , x 67 , x 75 , x 76 , and x 77 , and the location of support y 0 is changed to ensure it covers the entire study region by using the system sampling scheme. The connection relations between support x 0 and other supports in the research area can be obtained using Equation (7). The region inside the black matrix represents supports x 0 . The dark gray area represents the area connected to supports x 0 ; that is, the connectivity function C(x 0 , y j ) = 1, and the white region represents the area that is not connected to supports x 0 ; that is, the connection function C(x 0 , y j ) = 0. The connections between supports are independent of the relative locations of the supports in the process of system sampling. In other words, the connectivity of any two supports depends only on the connections between fractures and not on the relative locations of two supports. This way avoids the problem that the relationships between fractures (eg, intersection and closeness) cannot be (8) C N x 0 ,y 0 = 1 N ∑ N C i x 0 ,y 0 (9) C R x 0 ,y = ∫ R C x 0 ,y j ,x 0 ,y j ∈ R ≈ ∑ j C x 0 ,y j ,x 0 ,y j ∈ R F I G U R E 1 0 Support-based connection in 3D space F I G U R E 1 1 Numbering of supports during system sampling fully preserved due to the discretization of fracture network. The proposed method can preserve the relationships between fractures, and this capability is important in predicting the potential flow pathways and determining optimal drilling locations for reservoirs. Similarly, the connection relation of any two supports in 3D space is only dependent on the connections between fractures, and this condition is consistent with that in 2D space. Supports x 0 are positioned at x 222 , x 333 , and x 444 , and the connections between supports are shown in Figure 13. When multiple fracture networks in the research area are generated by MC simulation, the connectivity function between supports x 0 and y j can be calculated as where N is the number of MC simulations. For a specific fracture model, the GC of support x 0 in multiple fracture networks can be evaluated through the above mentioned Equation. A numerical example is given to demonstrate the connectivity function (10) under the assumptions that the scale of the research area is 100 m × 100 m × 100 m and fracture locations follow random Poisson distribution. The number of fracture (10)

F I G U R E 1 3 Connectivity algorithm
in 3D space surfaces is 150, and fracture direction obeys random distribution. Support x 0 is positioned at the center of the research area. The MATLAB program written by the authors is used to eliminate the isolated or closed fractures in the network, as shown in Figure 14B. The distribution of fracture clusters in the study area can be observed from the figure. The isolated fracture surface is gray, and the fracture surface belonging to the same fracture cluster is of the same color. A few studies have shown that whether the seepage reaches the critical point of percolation in the 3D fracture network is related to not only the density of fracture surface but also the maximum fracture cluster. Therefore, the approximate flow region in the fracture network can be obtained by separating the isolated fractures and determining the maximum fracture cluster. Figure 14C shows the GC of support x 0 in 3D space, and its contour maps are shown in Figure 14D,E,F. These maps are generated using 5, 50, and 500 times of MC simulations. The dark blue region in the figure is the disconnected region, and its connectivity function value is close to 0. The dark red region is the connected region, and the value is close to 1. The GC presents isotropy and more realistic simulation results as simulation times increase. The results show that the probability of connecting between supports is only related to the distance between them when the fracture orientation is random distribution, and the correctness of the program is proven. Figure 14G,H,I shows the connected region obtained according to GC. When the fracture orientations in the study area are along a certain direction, the fracture network model will show anisotropy. We assume that the fracture plane is approximately perpendicular to z 0 y plane in the simulation process to show the anisotropy of fracture network intuitively. The inclination angles of fracture surface are also set to 0°, 45°, and 90°, as shown in Figure 15. Equation (10) is used to analyze the GC of support x 0 in three different fracture networks, as shown in F I G U R E 1 4 Numerical example in 3D space Figure 16, for studying the influence of fracture orientations on the GC. Figure 16A,B,C shows the elevation, top view, and side view of the generated GC, while the inclination of fracture surface is approximately 0° in the 3D fracture network. The other graphs follow the same pattern. Figure 16 shows that fracture orientations greatly influence the GC. The connectivity is relatively high in the area along the orientation of fracture surface and relatively low in the area perpendicular to the orientation of fracture surface. This finding indicates that interconnected regions are easy to form on the orientation of fracture. The connected regions also extend along the orientation of fracture and present anisotropy. The results show that the contour lines will form an elliptic cluster centered on support x 0 with the increase in MC simulation times.
On the basis of the aforementioned simulation results, Figure 17 shows the connectivity distribution in different paths. The paths are parallel to the x-axis, y-axis, and z-axis and cross the center point (support x 0 ) of the simulation region, as shown in Figure 17A. When the fracture orientation is 0°, paths 1 and 2 have roughly the same connectivity distribution. The connectivity at both ends of path 3 shows a rapid decrease. This condition indicates that the probability of the fracture plane intersecting in the z-axis direction is small; that is, the probability of forming the fracture channel is low. When the fracture orientation is 90°, the connectivity distribution of paths 1 and path 3 is roughly the same, and both ends of path 2 show a rapid decrease. Figure 17C shows that the three curves slightly differ when the fracture orientation is 45°. All these phenomena indicate that the orientation of fracture significantly influences the connectivity distribution.
The figures indicate that Equations (8)-(10) are used to evaluate the connectivity of a fixed support x 0 . In evaluating the connectivity of the entire research area, the location of support x 0 in the aforementioned Equations is changed to ensure it covers the entire study region by systematic sampling distribution. The connectivity function of supports covering the entire study region is where k is the number of supports satisfying x k ∈ R and j is the number of supports satisfying y j ∈ R. The three-dimensional connectivity field (3DCF) can be evaluated by constantly changing the locations of support x k and support y j in Equation (11). The 2D example in Figure 12 is calculated to describe the function of Equation (11) simply and intuitively. Figure 18A,B shows the 2D fracture network and 2DCF calculated by Equation (11). The dark red region indicates that the connection of the region is good, while the dark blue region indicates that the connection of the region is poor. The connection relation only depends on the intersection relation between the fractures. Figure 18C shows the contour map of 2DCF. The map can directly display the connectivity distribution in the study area and predict the potential flow pathways. As shown in  Figure 18D, the area inside the blue dotted circle is considered to be very likely to change its own seepage state; that is, the region is the fracture expansion area. This information is useful in identifying the potential connections in fracture networks, such as fracture extensions in geothermal reservoir applications. Even though the fractures in the area within the dotted circle are not connected to each other, they play as potential paths in the flow. The GC can predict the potential flow path of a DFN, and this capability is important in practical engineering. Figure 19A-K shows the 3DCF of a random example. The scale of the study area is 100 m × 100 m × 100 m, and the number of fractures is 150. The centroid of the fracture follows random Poisson distribution, and the size of the fracture is distributed according to the power law with the value range of 1-20 m 2 . The fracture orientations follow Fisher distribution with the Fisher constant k which is 0. Figure 19 shows the resulting 3DCF slice and volumetric maps. The simulated 3D fracture network is shown in Figure 19A. Figure 19B,C shows the 3DCF of the study region, and the area with large connectivity, which can be directly observed. Figure 19D,H shows the section and outer surface of the contour map of 3DCF in the study area. Figure 19E-G and I-K is three views of them. In the figures, the red region has the best connectivity, while the blue region has the worst connectivity. The 3DCF in the study region presents anisotropy. Figure 20 shows the connectivity distribution curve on each path, which exhibits great randomness.
We assume that the geometric parameters of fracture surface obey a certain distribution through analysis of outcrop fracture data. 32 Therefore, the Monte Carlo method can be used to calculate the 3DCF in the study area. It can be obtained from Equation (11): where N is the number of MC simulations, and other parameters are consistent with those in Equation (11).
In summary, this study proposes a 3D connectivity assessment method based on the 2D connectivity assessment method proposed by Xu 25 and Allard. 31 The results for verification are shown in Figures 12, 13, 18, and 19. The support-based connectivity evaluation can preserve the spatial relation information between fractures. Therefore, the proposed method can provide highly reliable calculation results because it directly evaluates the connection of fractures, and this feature directly affects the seepage flow in the study area. The definition of CF and its extensions are generic and can thus be applied to 2D and 3D cases. The method has also been verified in previous studies of 2D fracture network calculation examples. 25,31 (12)

| CASE CALCULATION
The fracture data used in this study are based on the outcrop survey around Tongchuan City in the northern part of Ordos Basin. 21 Ordos Basin is located in the western region of China and has the characteristics of large area, resource potential, large reserves, and large scale. The prospect of oil exploration and development is immeasurable, the terrain in this area is relatively gentle, and the fracture distribution in this area is very representative. Figure 21 shows the topographic satellite map around Tongchuan City. It can be seen from figure that there are three major geological structural features of rock mass, such as uplift, fractures, and unconformity. The fracture data are shown in Table 1.
In this research, a 10 m × 10 m × 10 m modeling volume based on the initial geometric parameters of fractures in Table 1 is generated for connectivity analysis. To verify the validity of the fracture model, three parameters are extracted from the simulated fracture network model for comparison: (a) fracture areal intensity parameter P32, which is termed conductive intensity and used in some studies to quantized fracture model 20 ; (b) fracture trace length; and (c) fracture orientation. The simulated fracture parameters are obtained using Monte Carlo method to simulate 10 000 times, and the relative error between the simulation and measured results is shown in Table 2. The relative errors observed in the table are less than 5%. Therefore, the model can be considered to meet the calculation requirements. Figure 22A shows a randomly generated fracture network model. The blue fracture is fracture set 1, and the red fracture is fracture set 2. Figure 22B,C,D is their three views. The figures show the distribution and connection of fracture plane in 3D space.
The supports are distributed according to system sampling in 3D space, and their size meets the calculation requirements, that is, 1/25 of the boundary length. The connectivity field of the 3D fracture network can be calculated using the connectivity function, as shown in Figure 23. Figure 23A-G shows the resulting 3DCF slice and volumetric maps. These maps clearly demonstrate several well-connected areas in the 3D fracture network (areas shown in red). These areas can F I G U R E 2 2 3D fracture network model also be considered the domain that controls the flow characteristics. The 3DCF can be used to determine the existing and possible flow areas in the 3D fracture network. Figure 23I shows the flow region in the fracture network determined on the basis of the 3DCF. The 3DCF can be used to determine the optimal drilling locations for reservoirs, such as establishing a new well in the region with the largest connectivity value to maximize the productivity of reservoirs. The section of the flow areas determined based on the 3DCF is presented in Figure 24, which clearly shows the dominant flow path in the study area. The results show that the 3DCF can depict the connectivity of the research area in an all-round way. Fracture density is an important parameter in the evaluation of DFN. Therefore, it is necessary to compare connectivity with fracture density. Figure 25 shows the relationship between connectivity and fracture density, and the two parameters are dimensionless. Note that the fracture density refers to the total area of the fractures in each support. It can be found that the two parameters are closely related, and the curve changes in the same direction. The difference is caused due to the existence of isolated fractures. The results show that compared with fracture density, the connectivity is better to eliminate the influence of isolated fractures. Further studies need to be proved due to the limitation of simulation times.

| CONCLUSIONS
The connectivity evaluation is a highly effective method, especially in assessing the characteristics of reservoirs (eg, flow characteristics) on the basis of fracture systems. Most of previous studies are based on 2D space. Therefore, on the basis of Allard's definition of connectivity index, a 3D connectivity evaluation method is proposed by building the connectivity functions of the supports. This method depends only on the connections between fractures and is more realistic in assessing connectivity than the latticebased method that is defined only by the relative locations connectivity and fracture density of the supports. Meanwhile, some connectivity functions are built for different situations and used as deterministic methods for a fixed fracture network in 3D space and as probabilistic methods for a stochastic case. In addition, a convex polygon is used to simulate the fracture to make the fracture network as real as possible, and some examples are given to verify the correctness of the program and method. An example is calculated to compare connectivity with fracture density. The results show that the 3D connectivity evaluation method is feasible and more realistic and comprehensive than the 2D one. Thus, this method greatly improves the reliability of connectivity evaluation and plays an important role in the placement of pollution sources and the prediction of potential flow pathways, such as the placement of nuclear waste in areas with less connectivity value to reduce nuclear contamination. The 3DCF can be used to determine optimal drilling locations for reservoirs, such as establishing a new well in the region with the largest connectivity to maximize the productivity of reservoirs. This method also has certain theoretical reference importance to the flow problem of fractured rock mass encountered widely in hydraulic engineering.