Impact of Fracture Geometry and Topology on the Connectivity and Flow Properties of Stochastic Fracture Networks

In a low permeability formation, connectivity of natural and induced fractures determines overall hydraulic diffusivity in fluid flow through this formation and defines effective rock permeability. Efficient evaluation of fracture connectivity is a nontrivial task. Here, we utilize a topological concept of global efficiency to evaluate this connectivity. We address the impact of key geometrical properties of stochastic fracture networks (fracture lengths, orientations, apertures, and positions of fracture centers) on the macro‐scale flow properties of a shale‐like formation. Six thousand different realizations have been generated to characterize these properties for each fracture network. We find that a reduced graph of a fracture network, which consists of the shortest paths from the inlet nodes (fractures) to all outlet nodes, contributes most to fluid flow. Three‐dimensional (3D) fracture networks usually have higher global efficiency than two‐dimensional (2D) ones, because they have better connectivity. All geometrical properties of fractures influence quality of their connectivity. Aperture distribution impacts strongly global efficiency of a fracture network and its influence is more significant when the system is dominated by large fractures. Fracture clustering lowers global efficiency in both 2D and 3D fracture networks. Global efficiency of 2D and 3D fracture networks also decreases with the increasing exponent of the power‐law distribution of fracture lengths, which means that the connectivity of the system decreases with an increasing number of small fractures. Realistic fracture networks, composed of several sets of fractures with constrained preferred orientations, share all the characteristics of the stochastic fracture networks we have investigated.

developed open source software that generates 2D and 3D fracture networks with MAT-LAB. As a high-level programming language, MATLAB is insufficient to deal with hundreds of thousands or millions of fractures, especially in 3D. For Alghalandis (2017)'s software, it takes 23 s for a 2D network and 6 min for a 3D network with 100,000 fractures. In this research, we use the in-house developed C++ code. By applying a fast cluster check algorithm by Newman and Ziff (2001), our code can generate a large number ( 6 (10 )  ) of fractures in both 2D and 3D. Figure 1 shows the computational time of generating fractures and checking clusters.
In a low permeability formation, connectivity of natural and induced fractures determines the overall hydraulic diffusivity, and defines the effective permeability of fluid flow through the formation, (Bour & Davy, 1997;Maillot et al., 2016;Manzocchi, 2002;Renshaw, 1999). However, efficient evaluation of this connectivity is a nontrivial task.
Percolation theory is widely used to analyze connectivity of fracture networks (Berkowitz, 1995;Berkowitz et al., 2000;Bour & Davy, 1997, 1998de Dreuzy et al., 2001;De Dreuzy et al., 2002;Ji et al., 2011;Masihi et al., 2007;Odling, 1997;Robinson, 1983). In percolation theory, clusters of fractures are identified and labeled. In shales, connected fractures are commonly considered as the only permeable pathways for fluid flow. A spanning cluster of fractures is the largest cluster that connects boundaries of a domain under consideration, and creates global connectivity of the fracture network. However, a parameter that measures percolation is required before predicting the formation of the spanning cluster. When the spanning cluster forms in the system, the corresponding percolation parameter reaches the percolation threshold. For simple fracture networks, where fractures are distributed uniformly and have constant lengths, quantities like the total excluded area, total self-determined area and number of intersections per fracture are valid percolation parameters (Balberg et al., 1984;Bour & Davy, 1997, 1998Park et al., 2006;Robinson, 1983;Zhu et al., 2018). For fracture networks with variable lengths and non-uniform distributions of fracture centers, we found that these quantities are not percolation parameters (Zhu et al., 2018). It means that (a) these quantities cannot predict the formation of the spanning cluster, and (b) they cannot quantify fracture connectivity. Allard (1993) proposed a connectivity function to describe connections between two points, and further developed a cell-based method to evaluate the connectivity of the domain that has been discretized into pixels (cells) (Renard & Allard, 2013   approach to describe connections between two cells through fractures, and proposed a connectivity index/ field to evaluate the connectivity of fracture networks. Their binary connectivity function uses 1 and 0 to indicate whether two cells are connected through fractures or not. This function is better than a percolation parameter, such as intersection intensity, because it incorporates the density, intersection, and clustering information together. However, this function does not incorporate the apertures and variable lengths of the fractures, and it cannot evaluate how the different aspects of fracture connectivity impact the overall hydraulic conductance of the system. Intersection relationships of fractures can be classified as the fracture intersections (X-nodes), abutments or splays (Y-nodes), and the isolated fracture tips (I-nodes). Barton and Hsieh (1989) introduced a ternary diagram to characterize connectivity, in which the relative frequencies of the three intersection (node) types present in a system are plotted as points. Barton and Hsieh (1989)'s method captures topology of the intersections, and is adopted by many researchers to characterize fracture networks (Igbokwe et al., 2020;Manzocchi, 2002;Sanderson & Nixon, 2018;Sanderson et al., 2019). However, it neglects to capture how the various geometrical properties of a fracture network, such as fracture lengths and apertures, impact the hydraulic connectivity of the network.
In this study, we adopt a concept from graph theory, called global efficiency, (Hope et al., 2015;Valentini et al., 2007), to quantify the global connectivity of a fracture network and investigate how the different aspects of fracture geometry and topology impact this connectivity (Andresen et al., 2013;Huseby et al., 1997;Latora & Marchiori, 2001). In Table 1, we compare the global efficiency using all three commonly used approaches to demonstrate its advantages. "Yes" for each parameter (length, aperture or center position) means that a given connectivity index can quantify the indicated property of a fracture network. The global efficiency index is the only quantity we are aware of that captures the impacts of all geometrical properties of arbitrary (within confines of this study) fracture networks.
To apply a topological analysis, one needs to convert a fracture network to its graph representation. Flow calculations on a fracture graph, (Hyman et al., 2017;O'Malley et al., 2018;Viswanathan et al., 2018;Yang et al., 1995), are far more efficient than those using the discrete fracture-matrix models of flow in fractured porous media (Berrone et al., 2013;Botros et al., 2008;Karra et al., 2015;Roubinet et al., 2010).
The graph-based flow calculation neglects the impact of the matrix. For formations with low permeability, this simplification is valid, since fractures contribute the majority of the flow (de Dreuzy et al., 2001;De Dreuzy et al., 2002). To quantify the importance of matrix flow, we have performed a full-scale, embedded discrete fracture model simulation with UNCONG software (Li et al., 2015). We have calculated flow rates under the same boundary conditions, but with different ratios (K m /K f ) of matrix and fracture permeability in seven 2D and 3D benchmark cases. The resulting flow rates in each scenario are shown in Figure 2 and Table 2. Table 2 indicates that if the permeability ratio between the matrix and fractures is smaller than 10 −5 , the fracture network contributes more than 90% of the total flow. Fractures commonly have permeability of darcies or tens of darcies, while matrix permeability varies with rock type and depositional environment. For shale reservoirs, the matrix permeability is measured in tens of nanodarcies (Cho et al., 2013).
The purpose of this study is to investigate how the three key geometrical properties of fractures (lengths, apertures, and center positions), and the corresponding topological structures impact connectivity of a ZHU ET AL.

Materials and Methods
In this section, we discuss how to generate fracture networks in 2D and 3D, and how to convert them to the corresponding graphs. We also calculate flow through fracture networks and the measure of their global efficiency.

Generation of 2D and 3D Fracture Networks
The multiscale subsurface fracture networks are largely unknowable in detail, and it is difficult to quantify their properties exactly. To reduce modeling complexity, line segments are used to represent fractures in 2D, and random convex polygons with four vertices are adopted in 3D. The four geometrical parameters that describe each fracture are length, orientation, apertures, and position of the fracture center. Each of these parameters can be characterized by a different statistical distribution.
The fracture lengths are characterized by a power-law distribution (Bour & Davy, 1997) where n(l)dl is the number of fractures with lengths in the interval [l, l + dl], α is the coefficient of proportionality and a is the power-law exponent. Bonnet et al. (2001) and Zhu et al. (2018) showed that this exponent ranges approximately from 2 to 3 in most cases. The probability of generating very long fractures decreases sharply as a increases.
The fracture orientations follow uniform distributions, because subsurface rocks several have many different sets of fractures formed during their geologic history (Barton et al., 1995).
The positions of fracture centers are sampled from a uniform or fractal distribution. The fractal spatial density distribution (Darcel et al., 2003;Meakin, 1991) introduces clustering effects in the fracture network. Darcel et al. (2003) and Zhu et al. (2018) showed that in real fracture networks fracture centers cluster.
The fracture apertures are assumed to be constant, proportional to the fracture lengths, or follow a uniform or lognormal distribution.
ZHU ET AL.  Examples of different fracture networks in 2D and 3D are shown in Figures 3 and 4. Each network has a spanning cluster that connects all four sides of the domain boundary in 2D networks and the six boundary faces in 3D networks. The spanning cluster is the only pathway of fluid flow across the entire system and determines the effective permeability of the system.

Graph Representation of Fracture Networks
Connectivity is an essential aspect of fracture networks that determines a flow rate under an imposed pressure gradient. A topological concept, global efficiency, can be used to evaluate the connectivity of a given graph. Therefore, the conversion of a fracture network to its graph representation is necessary. Furthermore, it is more efficient to calculate flow using the node-link formalism (Fadakar-A et al., 2013;Hyman et al., 2017;Patzek, 2001), rather than to calculate it directly with the finite difference or finite element ZHU ET AL.
10.1029/2020WR028652 5 of 15  methods. Graph representation preserves all (in 2D) or most (in 3D) of the topological structure of a fracture network.

Graph Representation of 2D Fracture Networks
For a 2D fracture network, the graph representation is based on the percolation backbone of the fracture network obtained after removing all dead-end fractures (Bour & Davy, 1997). For coding practice, it is worth noting that new dead-end fractures may be generated once the old ones have been removed. Therefore, the backbone-check algorithm needs to be iterated until no new dead-end fractures appear. The backbone networks of each network in Figure  The graph representations of 2D fracture networks are exact, and capture all of the geometrical properties of these networks, including the lengths, orientations, apertures, and positions of the fracture centers.

Graph Representation of 3D Fracture Networks
For 3D fracture networks, a simplified graph is incomplete and lacks some geometrical properties of the fractures. 3D fracture shapes cannot be mapped onto a graph without loss of information, simply because the arbitrarily oriented fracture planes and fracture intersections cannot be projected onto a 2D plane, while preserving their 3D features. However, a graph representation of a 3D fracture network preserves important ZHU ET AL.
10.1029/2020WR028652 6 of 15 aspects of its connectivity, enough to determine most of the flow properties of this network. Each polygon is reduced to its centroid, and the intersection segment of two polygons is reduced to its central intersection point. If two fractures intersect, their centroids are connected via the intersection point. As a result, the nodes in the graph representation include centroids of all fractures and their intersection points. After converting to the graph representation, dead-end nodes and links are removed, because they do not contribute to flow. The graph representations of 3D fracture networks are shown in Figure 6.

Flow Rate Calculation
To calculate the single-phase flow rate in a fracture network, Darcy's law is applied. If a fracture in 2D is approximated by two parallel plates, the fracture permeability and aperture are related by the cubic law (Polibarinova-Kochina & De Wiest, 1962;Snow, 1965;Witherspoon et al., 1979Witherspoon et al., , 1980  2 12 f a k ( 2) where k is the fracture permeability and a f is its aperture.
For 3D fractures, we also assume flow between two parallel plates. After a transformation, Darcy's law can be restated as (Priest, 2012) where Q is the volumetric flow rate (m 3 /s); a f is the fracture aperture (m); b denotes the third dimension of the fracture (m), set here to 1 for 2D and 3D networks; g is the acceleration of gravity (m/s 2 ); ν corresponds to the kinematic fluid viscosity (m 2 /s), set to 1 × 10 −6 m 2 /s; L is the fracture length (m); ΔH is the hydraulic head drop (m); C is the lumped hydraulic conductance (m 2 /s). The correlation between hydraulic pressure drop ΔP and hydraulic head drop ΔH is The mass conservation must hold at all graph nodes. Using Darcy's equation and mass conservation law, the hydraulic head, H j , of node j and the conductance, C ij , between nodes i and j can be related as ZHU ET AL.
10.1029/2020WR028652 7 of 15 where n is the coordination number of node j. We prescribe a constant pressure boundary condition for all networks. To make the pressure gradient physically realizable and constrain the Reynolds number to a realistic range,

3
(10 )  , the hydraulic head at the inflow boundary (the left hand side of the 2D networks and the left face of the 3D networks) is set to 20 m, and all the other boundaries (the remaining three sides in 2D networks and the five faces in 3D networks) have zero heads. This head difference yields a macroscopic pressure gradient of 2 kPa/m (0.0884 psi/ft). The flow calculation in this paper is calibrated with a flow simulator developed by Li et al. (2015). For a 2D fracture network, the flow calculation based on its graph representation is exact. However, the flow calculation for a 3D fracture network based on its graph representation is only approximate, because the graph representation cannot capture all geometrical properties of this network.

Global Efficiency Calculation
Different distributions of fracture geometries change the connectivity property of the resulting fracture network. In topology, global efficiency, E, can be used to measure the connectivity of a network (Latora & Marchiori, 2001). The global efficiency is defined as where d ij is the shortest distance between nodes i and j, N is the number of nodes. d ij is either the physical distance or a weighted distance. In this work, the weights are the reciprocals of the hydraulic conductances (resistances) of the connecting links.

Results
After choosing the distributions that describe the geometrical properties of fracture networks, a set of corresponding stochastic fracture networks (DFN realizations) can be generated. Examples of stochastic fracture networks in 2D and 3D are shown in Figures 3 and 4. Subsequently, we convert the stochastic fracture networks to their graph representations (Figures 5 and 6) and calculate their global efficiency that describes the connectivity of the system. After imposing boundary conditions, we calculate the flow rate according to Equations 2 and 3. For each considered case, the results are averaged over 6000 random realizations.

Reference Cases for 2D and 3D Fracture Networks
To investigate the impact of different probability distributions of fracture geometry on the connectivity and flow rate of fracture networks, a reference (base) case is needed. For 2D fracture networks, the reference case is the average of 6000 stochastic fracture networks, where 1. Fracture lengths follow a power-law distribution with the exponent of 2.5 2. Orientations follow a uniform distribution between 0 and π 3. Positions of the fracture centers follow a uniform distribution inside a 100 × 100 m square 4. Apertures are set to 500 micrometers.
The 3D reference case has the same distributions of lengths, orientations (both strikes and dips), positions of the fracture centers, and apertures. The system domain is a 100 × 100 × 100 m cube. The boundary condition was defined in Section 3 and the minimum length of the fractures is 10 m. The hydraulic head distributions in the reference 2D and 3D networks are shown in Figures 7a and 7b.
In the generated fracture networks, many fractures do not contribute to the macroscopic pressure drop. These fractures significantly slow down the calculation of global efficiency, because the algorithm has the complexity of  (Cormen et al., 2009) is applied to find the shortest (least resistant) paths from the inflow nodes to all outflow nodes, and produce a reduced graph representation of the fracture network. Similar concepts of reduced (pruned) graphs have been adopted by many other researchers to simplify complex fracture networks and reduce computational time Viswanathan et al., 2018). The hydraulic head distribution of the reduced graph representation is shown in Figures 7c and 7d).
For the reference 2D case, the number of nodes in the original graph and reduced graph are 90.3 and 60.3 on average, and the corresponding flow rates are 1.65 × 10 −4 m 3 /s and 1.57 × 10 −4 m 3 /s, respectively, with the same boundary condition. This calculation indicates that on average, 67% of the nodes contribute 95% of the flow. For the 3D reference case, the respective numbers of nodes are 87.8 and 68.0, and the flow rates are 2.42 × 10 −4 m 3 /s and 2.36 × 10 −4 m 3 /s. Therefore, 77% of the nodes contribute 98% of flow on average. This phenomenon is more significant when the system size increases, meaning that a reduced network can capture the majority of the flow in both 2D and 3D networks.
In the global efficiency calculation, the weight of each link is the reciprocal of its conductance (resistance) normalized by the minimum weight of all links in the network. Therefore, the global efficiency calculation ZHU ET AL. depends only on the topological structure of the graph and can be compared with other graphs as well. The average global efficiencies of the original and reduced graph representations are 0.0132 and 0.0164 for the 2D reference networks. For the 3D reference networks, these values are 0.018 and 0.020, respectively. The reduced networks have global efficiency higher than the original networks, because of the removal of the least-conductive links. The 3D fracture networks have higher efficiency than the 2D ones, consistent with the fact that it is easier to percolate in 3D, rather than in 2D.

Impact of Aperture Variation
For the reference cases, all apertures are considered constant (500 micrometers), while in real fracture networks apertures vary. To investigate the impact of variable apertures on flow and connectivity, more nodes are added along each link. The number of added nodes is defined as where add i N is the number of additional nodes of link i, l i is the length of link i, and round is the rounding function that forces  2. Aperture of each sub-link is log-normally distributed over (0.8, 1.2) × a mean , and a mean is the aperture of the link calculated in scenario 2.
The method to generate a random number in a given range following a lognormal distribution has been discussed in (Baghbanan & Jing, 2007). Different aperture distributions are shown in Figure 8.
The results for 2D and 3D networks with variable apertures are listed in Table 3.
Since long fractures usually contribute more to flow and are endowed with larger apertures in Scenarios 2-4, flow rates are higher there. However, in general, the flow rate is not a good indicator of the connectivity of the whole network because it strongly depends on boundary conditions, as well as on the number and permeability of inlet fractures. Instead, global efficiency is a better measure, which relies only on the topology of the network. The constant-aperture scenario has the highest global efficiency, while scenarios 3 and 4 have the lowest efficiency in both 2D and 3D. Scenario 2 has an intermediate global efficiency, since apertures of sub-links are the same, while only apertures of different links vary based on their lengths. Aperture variations can reduce global efficiency by a factor of two. Efficiency reduction in 3D fracture networks is more significant, which indicates that these networks are more sensitive to aperture variations.

Impact of Fractal Spatial Density Distribution
The fractal spatial density distribution causes fractures to cluster, as commonly observed in nature. In 2D, if the fractal dimension is 2, the fractal spatial density distribution reduces to a uniform distribution. If the dimension is smaller than 2, there will be fracture clustering. Similarly, in 3D, the corresponding limiting dimension for uniform fracture distribution is 3. In this study, we choose 1.5 and 2.5 for the fractal dimensions of the 2D and 3D fracture networks. The method to generate a fractal spatial density distribution is discussed in (Darcel et al., 2003;Meakin, 1991). The comparisons between networks with the fractal spatial density distributions and the reference networks are listed in Table 4.
With the same constant head boundary conditions, 56.6% of nodes contribute 85% of the flow in the 2D fractal cases, and 70.3% of nodes contribute 95% of the flow in the 3D fractal cases. The global efficiency of the fractal case is 55% lower than the efficiency of the reference case in 2D networks, but it is only 15% lower in 3D networks, which indicates that the 2D networks are more sensitive to clustering effects.
ZHU ET AL.   Figure 9. flow-rate ratio versus fracture length exponent in 2D and 3D fracture networks. 2D fracture networks with positions following uniform or fractal spatial density distributions have the legends, Uniform, 2D, and Fractal, 2D. 3D fracture networks have the similar legends, Uniform, 3D, and Fractal, 3D.

Impact of Power-Law Distribution
Field observations and analog experiments suggest that the fracture lengths follow a power-law distribution. Bonnet et al. (2001) and Zhu et al. (2018) showed that the power-law exponent varies between 2 and 3 for most cases. The impact of the power-law exponent on a flow-rate ratio (R flow ) is depicted in Figure 9.
The flow-rate ratio decreases with the increase of the exponent a for all systems. 3D networks always have higher flow-rate ratios than 2D fracture networks. The 2D fractal cases have the most significant decrease reaching 20% when the exponent a is 2.9. It means that for 2D fracture networks with clustering effects and dominated by small fractures, it is not appropriate to replace the original fracture network with its reduced version.
Global efficiency variations of the original network and four scenarios with various apertures for four systems are shown in Figure 10.
Different distributions of fracture lengths, positions, and apertures all have significant impacts on the global efficiency of resulting networks. A larger number of small fractures (a higher exponent in the power-law distribution) can lower this efficiency severalfold. The impact of aperture variations is more important when fracture lengths are large. Clustering effects always reduce global efficiency. However, the difference between uniform and fractal 3D cases is smaller, indicating that clustering affects 3D networks less.
ZHU ET AL.
10.1029/2020WR028652 12 of 15 Figure 10. Global efficiency versus the fracture length exponent in 2D and 3D fracture networks. Solid markers refer to the fracture networks with centers following a uniform distribution; Open markers denotes the fracture networks with centers obeying a fractal spatial density distribution. The fracture property scenarios are listed in Table 3.

Global Efficiency and Flow Rate in Realistic Fracture Networks
Real fracture networks seldom have all of their fractures follow the same distributions of lengths, orientations, and center positions. These networks may be composed of several sets of fractures (Holland et al., 2009;Zhu et al., 2019). Each set of the fractures has its own preferential orientation because of the stress history. An example of a fracture system at Ras Al Khaimah in the United Arab Emirates is shown in Figure 11a.
Assuming that the maximum principal earth stress, S 1 , is oriented in the east-west direction, realistic fracture networks with four sets of fractures are constructed as in Figures 11b and 11c. In realistic fracture networks, the first set of fractures consists of a few large tensile fractures that span the domain and have similar strike directions (north-south), which may be an analogy to the tensile joints in a fold structure caused by bending. The other three sets of fractures include conjugate shear joints, small tensile joints, and random shear joints. Conjugate shear joints have strike angles ranging around N60°E or N60°W, deviating by about 30° from S 1 . Small tensile joints have strike angles similar to the first set of fractures. Random shear joints have random strikes and dips because of the local stress anisotropy. Except for the random shear joints, all other sets of fractures have dips around 90°. Clustering effects are also present in this system. The detailed information of each set of fractures is listed in Table 5. The global efficiencies of realistic networks are listed in Table. 6. A reduced network can still capture the majority of fluid flow and has a higher global efficiency in realistic fracture networks. 82.2% of original nodes contribute 96.6% of the flow in 2D, and 73.8% nodes contribute 96.9% of the flow in 3D. The 3D networks have higher efficiency than the 2D ones. Variable apertures still have a strong impact on global efficiency and reduce it significantly. In general, realistic networks share the characteristics found in the stochastic networks analyzed in this work.

Conclusions
We have analyzed the impact of fracture geometry and topology on the flow properties and global efficiency of stochastic fracture networks in 2D and 3D. The key conclusions from our research are: 1. A reduced graph composed of the shortest (least resistant) path from the inlets to outlets contributes the majority of the flow rate. In most cases, a reduced graph can be used to replace the original graph to reduce computational time. However, for fracture networks with a fractal spatial density distribution dominated by small fractures, this replacement is inappropriate. 2. 3D fracture networks have higher global efficiency than 2D fracture networks, which implies that the connectivity of 3D networks is better. 3. Aperture distribution decreases strongly global efficiency and impacts flow distribution. Variable apertures can reduce this efficiency severalfold, and their influence is more significant in networks dominated by large fractures. 4. Fracture clustering reduces global efficiency both in 2D and 3D; however, 3D fracture networks are less sensitive to clustering effects. 5. The global efficiency of 2D and 3D fracture networks decreases with the increase of the exponent of the power-law distribution of fracture lengths. This means that a larger number of small fractures lowers the connectivity of the network.
ZHU ET AL. e L max denotes the maximum fracture length while the minimum fracture length is set to 0.1 L; a is the exponent of the power-law distribution.