Permeability of Three‐Dimensional Numerically Grown Geomechanical Discrete Fracture Networks With Evolving Geometry and Mechanical Apertures

Fracture networks significantly alter the mechanical and hydraulic properties of subsurface rocks. The mechanics of fracture propagation and interaction control network development. However, mechanical processes are not routinely incorporated into discrete fracture network (DFN) models. A finite element, linear elastic fracture mechanics‐based method is applied to the generation of three‐dimensional geomechanical discrete fracture networks (GDFNs). These networks grow quasi‐statically from a set of initial flaws in response to a remote uniaxial tensile stress. Fracture growth is handled using a stress intensity factor‐based approach, where extension is determined by the local variations in the three stress intensity factor modes along fracture tips. Mechanical interaction between fractures modifies growth patterns, resulting in nonuniform and nonplanar growth in dense networks. When fractures are close, stress concentration results in the reactivation of fractures that were initially inactive. Therefore, GDFNs provide realistic representations of subsurface networks that honor the physical process of concurrent fracture growth. Hydraulic properties of the grown networks are quantified by computing their equivalent permeability tensors at each growth step. Compared to two sets of stochastic DFNs, GDFNs with uniform fracture apertures are strongly anisotropic and have relatively higher permeabilities at high fracture intensities. In GDFN models, where fracture apertures are based on mechanical principles, fluid flow becomes strongly channeled along distinct flow paths. Fracture orientations and interactions significantly modify apertures, and in turn, the hydraulic properties of the network. GDFNs provide a new way of understanding subsurface networks, where fracture mechanics is the primary influence on their geometric and hydraulic properties.

The methods by which mechanical and hydraulic properties are calculated on DFNs vary widely between implementations, including boundary element (Tuckwell et al., 2003), discrete element (Lisjak & Grasselli, 2014), finite difference (Karimi-Fard et al., 2004), finite element (Kim & Deo, 2000;Lang et al., 2018), finite volume (Koudina et al., 1998), and virtual element methods (Benedetto et al., 2017). Combined methods incorporate both continuum and discontinuum approaches to handle solution of forces and displacements, and handling of dynamics and contact, respectively (e.g., Elmo & Stead, 2010;Lisjak et al., 2016;Marschall et al., 2017). DFN studies are instrumental in computing the permeability tensors of fractured domains, because connectivity and channeling in domains are explicitly defined by the intersections between discrete fractures. Effective or equivalent permeability tensors can then be used to upscale the influence of fracture networks on large-scale models, where representing all the fractures which are relevant to flow is not feasible (Bisdom et al., 2015;Durlofsky, 1991Durlofsky, , 2005. Upscaling requires the definition of a representative elementary volume which is sufficiently large and intensely fractured. This poses a challenge in sparsely fractured rocks where properties are highly scale-dependent (Azizmohammadi & Matthäi, 2017).

10.1029/2019JB018899
The complexity of real fracture networks presents a challenge for DFN modeling. Previous studies have used stochastic (Poisson) methods to create networks, where fracture properties, including location, orientation, aperture, and length, are selected randomly (Ebigbo et al., 2016;Leung & Zimmerman, 2012;Mourzenko et al., 2011;Saevik et al., 2013), or generated using fractal theory (e.g., Long et al., 1992). Subsurface fracture networks have many features which are not well reproduced by such methods, including the following: • (i) Anisotropic fracture orientations: Subsurface fracture networks form as a result of discrete deformational phases associated with tectonic events (Gillespie et al., 1993) or local disturbances such as tunnel excavation (Baek et al., 2017). These states of stress are almost never isotropic. This results in sets of fractures with similar orientations. Variations in offset within a set are the result of mechanical interaction between fractures. Stochastic models can contain sets of fractures and vary orientations by small amounts within these sets, in order to more closely resemble anisotropically fractured rocks. • (ii) Shape: Nonplanar geometries frequently arise in nature during mixed mode growth, particularly as a result of interaction (Olson & Pollard, 1989, 1991Pollard et al., 1982), whereas fractures in stochastic DFNs are typically represented as planar disks. Furthermore, real fractures rarely have uniform shapes (Frash et al., 2019). • (iii) Spacing: Pre-existing fractures influence how additional fractures may nucleate. For example, in a tensile domain surrounding a fracture, stress is reduced along its strike and increased at its tips (Kachanov, 1987), which has been shown to have significant influence on the density of fractures in networks (Davy et al., 2010). In layered rocks, fracture spacing is also controlled by layer thickness (Olson, 2004). • (iv) Aperture: Fractures have complex walls and can be filled with different materials, resulting in extremely heterogeneous apertures and flow distributions across large fractures (Black et al., 2016;Snow, 1969;Tsang & Neretnieks, 1998), particularly at intersections (Leckenby et al., 2005). Apertures vary in response to stresses (Lang et al., 2018) but also due to chemical effects including dissolution and precipitation (Lang et al., 2015), and thermal effects (Salimzadeh et al., 2018). In situ stress and fracture interaction may have a more significant influence on apertures than fracture length (Olson, 2003;Pollard & Segall, 1987;Renshaw & Park, 1997). • (v) Length: The distribution of fracture lengths in a population is usually described using power laws (Bonnet et al., 2001;Davy et al., 2010;Gillespie et al., 1993). Fracture length also correlates with fracture aperture (Olson, 2003;Renshaw & Park, 1997). This means larger fractures generally have much more impact on flow than do smaller ones. The exponent of the power law governing a network's length distribution determines the extent to which flow is divided between large and small fractures (De Dreuzy et al., 2001). • (vi) Intersections: Stochastic networks have fundamentally different topological properties compared to real networks (Sanderson & Nixon, 2015). In Poisson DFNs, intersecting fractures will always cross one another, whereas in realistic networks, both crossing and abutting are possible and are determined by the mechanics of the system. Rules governing whether fractures cross in DFNs have been shown to be a major control on final domain properties, particularly fracture intensity (P 32 , fracture surface area per unit volume) (Davy et al., 2010). • (vii) Three-dimensional effects: In order to respect the full physics of the subsurface, DFNs must be three dimensional. Two-dimensional networks consistently underestimate the connectivity and permeability of three-dimensional networks (Lang et al., 2014). Mechanical interaction between fractures can increase or decrease stress along different parts of a fracture's tip, depending on their relative location and the orientation of stress, with significant implications for three-dimensional growth patterns (Thomas et al., 2017(Thomas et al., , 2020b. In literature, the term "discrete fracture network" (DFN) has been applied to describe fracture networks which have been stochastically generated, mapped from outcrops, or numerically simulated (Lei et al., 2017). In this work, the term "geomechanical" DFN is reserved for describing numerically generated DFNs that take into account the physical deformation of the rock matrix to dictate the growth of multiple fractures within the domain. Geomechanical networks are distinguished here from pseudo-geomechanical networks, which instead use heuristical approaches to emulate the formation of a network. Mechanical influence can be incorporated into Poisson DFNs by using power law distributions of fracture lengths (Bonneau et al., 2016;Thovert et al., 2017), alternative shapes for fractures (Cvetkovic et al., 2004), and correlations between fracture length and aperture (Baghbanan & Jing, 2008;Ebigbo et al., 2016). More advanced approaches use pseudo-mechanical rules to control different aspects of fracture insertion. Aspects of network development 10.1029/2019JB018899 which are dynamically evaluated during DFN generation include fracture nucleation, extension, termination, and mechanical interaction (Bonneau et al., 2016;Davy et al., 2013;Josnin et al., 2002;Maillot et al., 2016). These approaches do not directly simulate the mechanics of rock deformation or fracture growth and can therefore rapidly generate domains containing many fractures.
Examples of two-dimensional fracture growth have demonstrated the potential for grown fracture networks to reproduce realistic network properties (Olson, 2004;Olson & Pollard, 1989;Philip et al., 2005;Renshaw & Pollard, 1994). Paluszny and Matthäi (2010) introduced a finite element method-based approach where automatic nonuniform meshing of DFNs allows growth in any 2-D direction. The method separates the fracture geometry from the domain discretization, automatically generating the latter to conform to the evolving fracture surfaces. This method also enables stress-dependent fracture aperture distributions to be computed, showing that assuming uniform apertures on fractures in DFNs may cause an overprediction of domain permeability.
The present work makes significant progress toward realistic fracture network generation through the use of a three-dimensional, finite element-based fracture growth simulator (Paluszny & Zimmerman, 2011;Paluszny et al., 2013). Thomas et al. (2020b) validated fracture growth under tension and compression, and describe the advanced computational framework enabling the generation of geomechanical DFNs (GDFNs). In contrast to Poisson or pseudo-mechanical DFNs, these networks are generated by continually evaluating deformation of the volumetric domain and tracking the resulting changes in the network geometry. The present work describes the growth of six dense geomechanical fracture networks and examines their geometric and hydraulic properties. Uniaxial tension drives the extension of multiple sets of initially small fractures, which grow to make the domain three orders of magnitude more intensely fractured. Dense networks arise as a result of the stress state and fracture interaction in the volume. The equivalent permeability tensors of these GDFNs are compared against two types of stochastically generated DFNs with equivalent fracture intensities, demonstrating the various differences in network hydraulic properties when growth is geomechanical.

Methodology
A detailed overview and validation of fracture growth using the methodology presented herein, implemented in C++ in the Imperial College Geomechanics Toolkit, is given in Thomas et al. (2020b). The method is based on linear elastic fracture mechanics theory and grows fractures quasi-statically in distinct growth steps, starting from an initial distribution of flaws in a three-dimensional volumetric domain (Growth Step 0). In each step, the growth process is split into five stages: discretization (section 2.1), deformation (section 2.2), stress intensity factor computation (section 2.3), evaluation of growth criteria (section 2.4), and fracture extension (section 2.5). Flow and equivalent permeability analysis at each growth step is performed during postprocessing after growth (section 2.6). The present work focuses on linear elastic tensile deformation and Darcy flow, and in other work the ICGT has also been applied to compressive stresses (Thomas et al., 2020b), friction and contact (Nejati et al., 2016), and coupled thermohydromechanics (Salimzadeh et al., 2018).

Geometry and Discretization
The simulation domain can have any arbitrary shape, and in the present work, the domain is a cube with 20 m sides. Two representations of the domain are handled during fracture growth simulations (Paluszny & Zimmerman, 2011). The first of these is the geometry of the domain, composed of simple geometric objects (e.g., disks, triangles, or squares) and sorted into groups (one group for each boundary and one group per fracture). Fracture surfaces are described as any number of connected two-dimensional surfaces. Disks provide a convenient initial fracture shape, because fractures which grow quasi-statically will converge toward circular shapes if unaffected by stress perturbations (Lazarus, 2003), aiding the analysis of final fracture patterns. When fractures grow, new surfaces (triangles) are added at the tip of a fracture, determined by the length and orientation of the growth vectors along the tip (see section 2.4).
In order to solve for forces and displacement throughout the domain, an octree mesher is used to automatically generate a mesh at each growth step. Two types of elements are used in the domain: isoparametric quadratic tetrahedra in the rock volume, and isoparametric quadratic triangles on surfaces (fractures and boundaries). The meshing process uses an optimization-based approach, dynamically sizing volumetric and surface elements to place additional discretization in fracture tip regions. This reduces the total number of degrees of freedom for the finite element computation, while maintaining accurate stress and displacement fields near the fracture tip. A typical mesh for step 0 of the GDFNs in this work contains up to 320,000 nodes and 250,000 elements.
Each fracture tip is discretized into a number of tip nodes, connected by tip segments. In this work, fractures have 20 tip nodes, meaning each initial tip length is the fracture's circumference divided by 20. The minimum 10.1029/2019JB018899 size for finite elements along the fracture tip is defined by the length of one tip segment. In regions further from the fracture tip, elements gradually increase in size. Fracture surfaces are meshed as two coincident surfaces which connect at the tip, shown in Figure 1a. During deformation, the two opposite sides of the fracture deform separately. Therefore, fractures are simulated as a discontinuity in the volumetric domain. A discussion of the tip size effect and quantification of its effect on the error in the stress intensity factor (SIF) is given in Thomas et al. (2017). Growth of a tip causes the new tip nodes to be further apart, increasing tip lengths and thus local element sizes as the simulation progresses. Simulations can also maintain tip lengths and thus keep refinement relatively constant.

Linear Elastic Deformation
The rock matrix is treated as a linear-elastic, homogeneous, isotropic material. The stress and deformation throughout the rock is solved using the finite element method. Stress and strain are governed by where is the Cauchy stress tensor, is the infinitesimal strain tensor, 0 and 0 are the initial stress and strain, and D is a linear elastic stiffness matrix (Cook, 2007). For a static system, where is the differential matrix operator, and F represents body forces exerted on the object. In this work, the rock matrix is a 20 × 20 × 20 m cube, and the remote boundary condition driving growth is a fixed stress applied over the top boundary (+Z). The bottom boundary (−Z) is fixed (no displacement) and the other boundaries are traction-free. Gravitational body forces are not imposed on the domain. Tensile stresses are reckoned positive in this work.

SIF Computation
SIFs are defined along fracture tips and quantify the energy that deforms the tip in the three propagation modes: opening (mode I), in-plane shear (mode II), and out-of-plane shear (mode III) (Anderson & Anderson, 2005;Lawn, 1993). Each of the modes has a corresponding SIF (K I , K II , and K III ). SIFs are a function of the stress and displacement close to the tip and are influenced by fracture geometry, including inclination relative to the load direction and fracture interaction (Thomas et al., 2017). Locally evaluating SIFs around a tip enables growth criteria to be evaluated, which provide the means of extending fractures (see section 2.4) (Paluszny & Zimmerman, 2011).
The present work adopts the displacement correlation (DC) method for computing the SIFs at each tip node along a fracture's tip, by sampling two points on the fracture's surface (Nejati et al., 2015). A plane strain condition is known to prevail locally around a fracture, such that three-dimensional deformation fields approach two-dimensional plane strain fields (Nakamura & Parks, 1988). Williams (1957) expressed the linear elastic stress fields for a two-dimensional cracked plate subjected to an arbitrary load, now known as the Williams series expansions (Anderson & Anderson, 2005). Near the tip, the first terms of the Williams series are dominant, providing an expression for the singular stress field near the tip. Thomas et al. (2020b) validated DC for fractures under tension and compression for all three SIF modes within the present implementation of the method.
Large variations in the SIFs along a fracture's tip in complex domains, such as the present GDFNs, can result in jagged fronts which do not converge to circular shapes. This effect is noticeable when large parts of a fracture's front are handled with relatively few fracture tips. A NURBS-based (Non-Uniform Rational B-Spline) interpolation method is employed to smooth small variations in K I , K II , and K III along the fronts of single fractures. A two-dimensional second-degree rational B-spline curve approximates each modal stress intensity factor, whereby each point of the spline is defined as (x, K i ), where x is the one-dimensional distance along the tip curve, and K i is the modal stress intensity value. The result is a curve which is smooth and is defined to smoothly approximate all the original SIF values, for each modality in an independent manner. The curve is defined to have a low number of control points, so as to reduce variability in the SIF prediction. The resulting distribution of SIFs is similar to the unsmoothed SIF distribution and ensures that high values are spread out over multiple tip nodes. Similar approaches are used to improve growth in other SIF-based growth schemes, for example, moving averages (Gupta et al., 2017).

Growth Criteria
Following SIF computation, three criteria are evaluated to determine the local behavior of the fracture tip. A mixed mode failure criterion is adopted, such that a fracture tip will grow if the comparative stress intensity factor (K v ) is the same as, or greater than, the critical SIF of the material (K IC ) (Richard et al., 2014). The comparative SIF is a combination of all three SIFs, and is calculated as follows: where 1 and 2 are the ratios of the mode I fracture toughness to the modes II and III fracture toughness To match the properties of brittle low-permeability rocks, K IC = 1.8 × 10 6 Pa m 1/2 , K IC ∕K IIC = 1∕1.3268 and K IC ∕K IIIC = 1∕0.6297 (Aliha et al., 2015;Chang et al., 2002;Wei et al., 2016) are chosen in the present work.
For each growing fracture tip, a vector is generated to extend the fracture tip. This requires two further values: an extension length, and an extension angle ( ). This angle is measured between the fracture plane and normal to the fracture plane, such that 0 o growth is along the plane, and 90 o is normal to the plane (Thomas et al., 2020b). Analysis of fracture growth experiments identify the following expression as being accurate for many types of three-dimensional, mixed mode loads: where A = 140ř,B= −70 • . Angle is negative when K II > 0 and positive for K II < 0 and K I ≥ 0 (Richard et al., 2004(Richard et al., , 2013. The twisting angle is not incorporated to avoid simulating mode III breakdown (e.g., Martel & Boger, 1998), which may cause complex tip geometries to emerge which are more challenging to discretize dynamically than smooth surfaces.
Fracture extension lengths are assigned based on the value of the strain energy release rate () at each tip, which is a combination of the three SIFs: where subscript n refers to one particular tip node, E ′ is the Young's Modulus for plain strain conditions (i.e., E ′ = E∕(1 − 2 ), where is the Poisson's ratio (Nakamura & Parks, 1988), and is the shear modulus of the material (Rice, 1968). A modified two-step Paris-type law is implemented to handle concurrent fracture growth in three dimensions (Thomas et al., 2020b). The approach is based on the Paris law for fatigue fracture growth and is modified to improve its ability to model the simultaneous growth of multiple fractures and fracture tips. In each growth step, Δa v is the maximum extension length for all tips in the volume (2 m in the present work). For tip node n on fracture f , an extension length Δa n is selected by applying two power laws with different exponents, providing a length less than or equal to Δa v . First, a maximum extension length for each fracture (Δa f ) is found: where  is the largest value of  n on the fracture f ,  v is the largest value of  n in the volume, and f is the fracture growth rate exponent. Subsequently, the local tip extension is assigned, scaled by this value: where n is the tip growth rate exponent.
The use of the modified two-step Paris-type law is instrumental for growing geomechanical networks, for several reasons. In many brittle geological materials (e.g., in outcrops, reservoirs, and deep crystalline rock), fractures are found to have grown together in populations under the same stress conditions (Segall, 1984;Lamarche et al., 2017). A standard Paris-type law in three dimensions focuses growth on only the most stressed tips, leading to the activation of relatively few fractures in the network, unless very low values of the growth exponent are used. The two-step approach models wider activation of fractures, enabling more fractures to grow competitively and concurrently in each growth step. Concurrent growth is also favorable for reducing the computational time required to reach high fracture intensities. The values of f and n used to generated the present networks are 0.35 and 2, respectively, providing conditions where fractures activate widely in the domain, and mechanical interaction affects fracture growth extents. The choice of these parameters is examined further in Thomas et al. (2020b).

Fracture Extension
At each node which meets the mixed mode failure criterion (equation (3)), a growth vector p is defined with length Δa n . The vector's direction is along the direction of growth in the previous step, modified by angle (equation (4)). An example growth vector shown in the context of the finite elements at a fracture's tip is shown in Figure 1a. The extended fracture geometry is updated by creating a new fracture tip loop (constructed from connected tip nodes) along the points at the end of each vector, or through the original tip node if no growth occurs. This is constructed as a strip of triangles, with two triangles tessellating the space between each vector (Thomas et al., 2020b). Figure 1b shows how growth vectors provide the sides and corners of the new triangular surfaces. The new fracture geometry is brought forward to the next step, and subsequently the domain is remeshed. The new mesh is constrained in the same way as the original, that is, the mesh is finer around the fracture tip and gradually coarsens with distance, and finite elements do not conform exactly to the extended fracture surface.

Fracture Intersection
Intersecting fractures can form either "T" or "X" intersections, depending on whether fractures terminate against one another or cross, respectively (Sanderson & Nixon, 2015). Empirical criteria attempt to model fracture behavior at the intersections; for example, Renshaw and Pollard (1995) for general frictional surfaces in linear elastic materials, and Cheng et al. (2014) for hydraulic fractures. These criteria require additional sampling of local stress fields at each intersected crossing fracture at every step, adding significant computational cost to each growth step. In three dimensions, further complexity is introduced, where a fracture may cross at one point but terminate in another. Furthermore, fracture intersections between nonplanar fractures cannot be characterized by a single angle for use in existing crossing criteria. Therefore, a straightforward approach to handling fracture intersections is adopted here, with the intention that future investigations will examine the influence of both "T" and "X" intersections in geomechanical networks. This approach is consistent with the use of a tensile boundary condition, as no stress can be transmitted across the plane of a fracture whose surfaces are not in contact.
In the present work, tips that intersect boundaries and fractures are deactivated so that they do not propagate further. Intersections are detected by checking whether growth vectors intersect surfaces in their local neighborhood. An example of intersection during growth is shown in Figure 1b, where two tip nodes terminate on another surface. If some tip nodes are intersected, the other tip nodes of the fracture can continue to grow.
If a tip is close to the boundary, in order to avoid performing calculations at volume boundaries where values become unknown, tips are "snapped" onto the boundary along the most direct path. To handle the special case where two growth vectors intersect different boundaries, a single triangular surface is created which fills the corner. Therefore, when new surfaces are added to the fracture, modifications to the growth vectors ensure that fractures do not go out of bounds and that "T" intersections are created when fractures first intersect. "X" intersections can arise in the present implementation of this methodology in two ways. First, the initial flaws can crosscut one another due to random placement. This "X" intersection will then remain in the simulation from the start. Second, if both growth vectors on one fracture tip do not intersect a fracture, but the line between the ends of the growth vectors do, the new surfaces are created using these vectors, and the new tip is deactivated. This case will result in a small part of the fracture crossing.

Simulation End Condition
In this work, GDFNs are grown until no additional new fracture surfaces are added in each step. This cessation of growth occurs when all tips have terminated against other fractures or boundaries, or if the mixed-mode failure criterion is not met at any tips (equation (3)). This implies that the network has reached the maximum possible fracture intensity for the applied boundary condition. In future work, changing the boundary condition (i.e., increasing the magnitude or changing the orientation of the stress or displacement) will enable additional fracture growth in the domain to occur.

Equivalent Permeability Computation
To investigate their hydraulic properties, the equivalent permeability tensor of the GDFNs are calculated for each growth step. The method described and validated by Lang et al. (2014) is implemented in the ICGT. The method begins by solving three separate single-phase Darcy flow problems using the finite element method. Each experiment uses opposite boundaries of the DFN domain as the source and sink for flow (−X to +X, −Y to +Y, and +Z to −Z), with source pressure 100 Pa and sink pressure 10 Pa, and all other boundaries are no-flow boundaries. The Darcy flow problem models the flow of an incompressible, single-phase fluid, where ▿p is the pressure gradient, and k e is the element permeability. In volumetric (matrix) elements, k e has a fixed value of 10 −15 m 2 . On fracture surface elements, a permeability is assigned based on the local aperture of the fracture, h (Zimmerman & Bodvarsson, 1996): For each flow experiment, element-volume weighted averages of the pressure gradients and fluxes (⟨u i ⟩) along the three coordinate axes (i = x, y, z) are computed over a subdomain for each flow problem (Durlofsky, 2005): ⟨ p where u is the flow velocity, V is the volume of the domain, e is a finite element, u e i is the flow velocity along direction i in element e, V e is the volume of element e, p is the pressure, and p e x i is the pressure gradient. Surface elements on fractures have zero volume, so a local thickness is assigned from the aperture of each fracture's element.
The quantities in equations (10) and 11 are related to permeability by Darcy's law for porous media, where k is the equivalent permeability tensor, with eigenvalues k 1 , k 2 , and k 3 , and is the viscosity (1 × 10 −3 Pa·s throughout this work) (Durlofsky, 2005). Equation (8) provides the terms for a linear equation of the form Ax = b, where A is a matrix of the average pressure gradients, x is the equivalent permeability tensor k (in vertical form as a 1 × 9 vector), and b is a vector of the average flow velocities. Tensor k is formulated as follows: where k xx is the permeability of the domain in the x direction, and so on. The off-diagonal terms describe the flow rate in one direction due to the pressure gradient in another. The terms in the matrices A and b are provided in full by Lang et al. (2014). This overconstrained linear equation is solved for x using LU decomposition, thus providing the equivalent permeability tensor of the fractured domain.
In order to avoid boundary effects, averages for equations (10) and (11) are computed only on elements with barycenters located within the central 90% of the domain. This approach has been referred to as undersampling or flow jacketing (Lang et al., 2014) and produces results which better reflect the overall connectivity of the domain (Lang et al., 2018). Use of Darcy's law assumes that the medium can be replaced by an equivalent homogeneous block of a small size with permeability k, known as a representative elementary volume. For fracture networks this requirement may not be met, as relatively small representative elementary volumes relative to fracture lengths do not have ellipsoid-shaped permeabilities (Long et al., 1982). In this work, we Young's modulus (E) 2 0 × 10 9 Pa (Krech et al., 1974) Poisson's ratio ( ) 0.2 (Krech et al., 1974) Critical stress intensity factor (K IC ) 1× 10 6 Pa m 1/2 (Wei et al., 2016) K IC ∕K IIC 0.7540 (Chang et al., 2002) K IC ∕K IIIC 1.5881 (Aliha et al., 2015) Matrix permeability (k m ) 1 0 −15 m 2 Fracture growth exponent ( f ) 0.35 (Renshaw & Pollard, 1994) Fracture growth exponent (  do not explore the sensitivity of our results to the scale of observation, and highlight that the terms of computed k tensors may change with the sampling scale. Furthermore, k should be considered the equivalent, rather than effective, permeability tensor for a particular domain, as effective would imply that a sufficiently large region has been sampled to encompass all the scales of variation in the permeability field (Durlofsky, 1991;Renard & de Marsily, 1997).

Geometry of GDFNs
This section describes the growth and geometric properties of six GDFNs. The six networks share the same material properties and stress boundary condition driving growth. They contain different numbers and distributions of initial fractures. Fracture intensity (also referred to as fracture density in literature) in the networks is quantified using P 32 (Sanderson & Nixon, 2015): where i is the fracture number, n f is the number of fractures in the domain, f ai is the surface area of fracture i, and V is the volume of the domain (8,000 m 3 ). P 32 is particularly suited to networks containing fractures with irregular surface areas. The domains are visualized with fractures represented as opaque planes. This rendering process simulates a light source such that some fractures will cast shadows onto others and different parts of fracture surfaces are darker or lighter depending on the viewing angle, to improve visualization of three-dimensional domains.  Table 1. P 32 values (m −1 ) quantify the fracture intensity (or density) of the domain (see equation (14)). Table 1 details the numerical properties, material properties, and boundary conditions used to grow GDFNs and analyze their hydraulic properties. Material properties for hard rocks (e.g., granite, granodiorite) were chosen to make results particularly relevant to radioactive waste disposal and other deep rock applications. Each of the six GDFNs is named with a letter (A-F), and the geometry of the domains in their initial and final growth step are shown in Figure 2. GDFNs A and B have 30 initial fractures, C and D have 60, and E and F have 90 fractures. Initial fractures are randomly placed in the central 7 × 7 × 7 m region of the domain, with no two fracture centers of the same inclination being less than 1 m apart. Domains with the same number of fractures have different random distributions. Once discretized, the number of elements in each mesh correlates with the number of fractures. GDFN E has the most nodes (321,041) and elements (247,187), and GDFN B had the least (121,677 and 93,359 respectively).

Initial Domains and Material Properties
The initial radii and orientations of the fractures have a significant effect on their growth. Five in every six of the initial fractures in the GDFNs have a radius of 0.5 m and are oriented with a vector normal to their surface of (0, 0, 1); one in six fractures has a radius of 2.5 m and normal vector (4, 5, 1). These vectors are relative to the coordinate axes shown in Figure 2. In each growth step, a tensile stress of = 2 × 10 6 Pa is applied, using the uniaxial tension boundary condition described in section 2.2. This stress state results in fracture growth on the initially horizontal fractures from step 0, as they immediately meet the failure criterion (equation (3)). Initially, the inclined fractures do not grow because their inclination reduces the SIFs at their tip. As fractures in the domain grow closer together, mechanical interaction leads to stress concentration near fracture tips, activating the inclined fractures along critically stressed parts of their tips. Contact and friction can be handled within the present methodology for shear and compression boundary conditions (Nejati et al., 2016;Thomas et al., 2020b). However, because domains are deformed in tension, fracture surfaces separate and tracking the contact and friction between surfaces is not required.

Growth of Geomechanical Networks
In step zero (Figure 2), all six GDFNs contain a random cluster of fractures. In the final steps, the fractures have grown across most of the volume, with many horizontal fractures intersecting the side boundaries. Fractures which intersect along some of their surfaces can become elongated, as growth continues on the rest of the fracture while the intersected tips remain deactivated and stationary. Most fracture growth in the GDFNs is horizontal (perpendicular to the stress direction), including the inclined fractures. Some surfaces are curved, particularly in later growth stages and near boundaries. In all GDFNs, overlapping fractures reorient their surfaces, propagating such that they overlap as little as possible (i.e., preferentially growing toward the boundaries, or curving toward less densely fractured regions of the domain). However, a significant amount of overlap still arises despite their preference to avoiding overlapping. Fracture tips near the traction-free boundaries experience interface stress concentration (e.g., Al-Ostaz et al., 1998;Pollard & Holzhausen, 1979), altering their SIFs, thereby promoting additional curvature by changing growth angles (equation (4)). The domain is relatively more displaced at the top of the domain, leading to larger fractures being more common near the top boundary.
GDFNs A and B contain 30 fractures and grew for a relatively large number of steps compared to the initially denser GDFNs (E and F). GDFNs A and B's relatively low initial fracture intensity results in more space being available for fracture growth. GDFN B reaches a much higher density than A, because its arrangement of initial fractures provides greater space for fractures to grow. GDFNs seeded with the most fractures tend to grow for fewer steps, but reach fracture densities similar to other GDFNs with fewer initial fractures. Domain densities could be increased by continuing the growth simulations, or increasing the magnitude or orientation of stress boundary condition.  Figure 3a shows the simulation step against two measures of fracture intensity (total fracture surface area and P 32 ) for GDFN B. Both values follow a curve with a reducing gradient during growth. In the early steps, the increase in fracture surface area is rapid, with approximately 80% of the surface area being present by the 20th growth step. In later steps, the growth slows down as fracture tips stop growing due to intersection and deactivation.
Changes in total fracture surface area added per step can be examined in Figure 3b, which plots the surface area of each fracture over simulation time. Fractures are colored by their initial inclination, and at Step 0, the difference in the initial surface areas of each set can be seen. The horizontal fractures grow quickly in the first 10 steps. The inclined fractures grow at a more steady rate throughout the simulation. Fracture deactivation occurs when no additional surface area is added, and the growth step at which this occurs varies significantly between different fractures. These fractures may be mechanically active in other ways (i.e., slip, or change aperture) even though they no longer propagate. Figures 3d and 3e show the surface area of the domain and individual fractures respectively for GDFN C. Both GDFNs B and C have similar trends of fracture surface area during growth, with fractures being generally smaller and deactivating earlier in GDFN C due to the higher number of initial fractures (60 rather than 30).
Fracture populations are often analyzed using histograms of fracture length (Armand et al., 2014;Fang et al., 2017;Hardebol et al., 2015). Due to their nonuniform shapes, surface areas of fractures in the final growth step in GDFNs B and C are used to create histograms, shown in Figures 3c and 3f. The present fracture surface area distributions are similar to a normal or log-normal distribution, unlike natural fracture networks which are characterized by power laws (e.g., Davy et al., 2010;Gonzalez-Garcia et al., 2000). Renshaw and Pollard (1994) also observe a log-normal distribution of fracture lengths in networks generated with two-dimensional fracture growth model. If relatively small fractures, or a wider range of fracture sizes, were included in the present GDFNs, selective fracture activation would cause most of them to remain inactive, due to smaller fractures acquiring lower SIFs (Kassir & Sih, 1975). This mechanism of concentrating strain toward a small number of large fractures has been observed in two-dimensional networks (Renshaw To generate the ellipsoids, eigenvectors (shown as red, green, and blue lines) are multiplied by their respective eigenvalues, and a surface is passed through the vectors. In the simulations that generated the equivalent permeability tensors, fractures have uniform apertures of 0.001 m. & Pollard, 1994). Including more small fractures would improve the realism of the network, in terms of the distribution of fracture sizes, at the significant cost of additional mesh refinement.

Equivalent Permeability of GDFNs
The equivalent k calculation method introduced in section 2.6 was used to quantify the hydraulic properties of the GDFNs with uniform apertures (section 4.1) and mechanically induced apertures (section 4.2). These apertures replace those which emerge during fracture growth. Both uniform and mechanically induced apertures are analyzed to distinguish which flow and equivalent permeability properties arise in each case, with the former condition being common in other DFN studies (e.g., Vu et al., 2018).

Uniform Apertures
In uniform aperture flow simulations for GDFNs and SDFNs, the matrix was assigned a permeability of 10 −15 m 2 , and all fractures are assigned constant apertures of 0.001 m, leading to a fracture permeability of 8.33 × 10 −11 m 2 . These parameters create a high contrast between flow in the fractures and in the matrix, meaning permeability is significantly influenced by the fracture network and its connectivity. Higher refinement meshes are used in the flow simulations to ensure that fracture intersections are discretized accurately. Figure 4 visualizes the geometry of GDFN E at four different growth steps, and the equivalent k tensor as an ellipsoid at each step. The ellipsoid axes correspond to the tensor eigenvectors and its axes half-lengths correspond to its eigenvalues. The ellipsoid serves to visualize the tensor in three dimensions with respect to the geometry of the domain. The two early growth steps shown (0 and 1) demonstrate that the tensor is initially relatively isotropic, with all three eigenvectors having similar magnitudes. Percolation (the emergence of a flow path between source and sink boundaries along fractures) occurs in GDFN E in growth Step 4, increasing the horizontal permeabilities (k xx and k yy ) from a magnitude of 10 −14 to 10 −12 m 2 . This is demonstrated by the ellipsoid in growth Step 5 (Figure 4c) which becomes almost flat in the x-y plane. The tensor becomes relatively more strongly anisotropic in the final growth step (19).
The hydraulic properties of GDFNs are compared to two sets of stochastic DFNs (SDFNs) to quantify their changes in permeability as fracture intensity increases. SDFNs are generated by randomly placing circular fractures into a 20 × 20 × 20 m domain, and are oriented such that the fracture normal vectors have a uniform distribution on a sphere. Different fracture intensities are achieved by varying the number of fractures 10.1029/2019JB018899 in each domain, and parts of those fractures which protrude outside the domain do not contribute to fracture intensity or flow. Due to the simple method used to generate the SDFNs, they exclusively contain X intersections.

Comparison to Isotropic Stochastic Networks (SDFN Set 1)
In SDFN Set 1, fractures have a fixed radius of 4.3 m. This radius is selected as it matches the approximate mean radius of fractures in GDFN C at its final step, assuming fractures are circular with radius r = √ a∕ , where a is each fracture's surface area. Fractures are randomly oriented with their centers placed in a region extending ±11 m from the center of the domain, allowing some fractures to intersect the boundaries. For clarity, SDFN and GDFN permeabilities are compared by plotting diagonal values of k, which correspond to the permeabilility of the network in the x, y, and z directions. Figure 5 compares the directional permeabilities of each growth step of the six GDFNs to 81 SDFNs from Set 1. Figures 5a-5f plot k xx , k yy , and k zz for each domain and growth step against their P 32 values, with the same results from SDFN Set 1 shown as colored dots in each graph. Examples SDFNs in Set 1 are shown in Figures 5g-5j with increasing fracture intensity. The GDFN directional permeabilities are shown as crosses connected by lines, and have several common features that occur in all six GDFNs. Initially, the GDFNs have isotropic permeabilities close to k m . The largest change in permeability occurs at Step 4 or 5 in all GDFNs, when k xx and k yy increase significantly, indicating horizontal percolation. The stress boundary condition promotes growth in the x-y plane, causing fractures to intersect the vertical boundaries and each other, creating percolating flow paths at around Step 4 or 5. This demonstrates the strong dependence of GDFN permeability on fracture growth.
Vertical percolation does not occur in any of the GDFNs, except B where a larger k zz arises in the final steps, and in A for one step. In these cases, fractures intersect the top boundary, but not the bottom, providing most of the required flow path, meaning the jump in permeability is not as significant as those in k xx and k yy . In the other domains, the vertical permeability k zz generally fluctuates because the surfaces added in each growth step do not provide pathways which promote vertical flow. This suggests that the permeability of each axis also depends on the permeability of the other axes. The P 32 at which percolation occurs varies significantly with the number of fractures in the GDFN. Domain properties such as the total number of simulation steps (i.e., the accumulated growth over time), initial fracture orientations, and the maximum tip extension per step are important in determining the percolation threshold of the present GDFNs.
Networks in SDFN Set 1 percolate at lower fracture intensities than the GDFNs. This percolation occurs over a range of P 32 values between ∼0.075 and 0.16 m −1 , which is a similar range to the GDFNs (0.15 to 0.25 m −1 ). In the final steps of growth, the GDFNs have higher k xx and k yy values than Set 1 SDFNs, due to the emerging anisotropy on the x-y plane.

Comparison to Anisotropic Stochastic Networks (SDFN Set 2)
In SDFN Set 2, fractures have a fixed radius of 4.3 m and have two orientations: five in every six have a surface normal vector of (0, 0, 1), and one in six have a surface normal vector of (4, 5, 1). Fracture centers are placed in a region extending ±11 m in the x and y axes, and ±9 m in the z direction, centered on the midpoint of the 20 × 20 × 20 m domain. SDFN Set 2 therefore more closely resembles the GDFNs due to the fact that majority of fractures are normal to the z axis and a small number of vertical fractures providing connectivity between the horizontal fractures in the domain.
Figures 6a-6f compares the directional permeabilities of the six GDFNs to 114 SDFNs from Set 2, and demonstrates the geometry of four domains at different fracture intensities (G-J). Set 2 SDFNs percolate at higher P 32 values than SDFN Set 1 (Figure 5), more closely matching the densities at which the GDFNs percolate. Vertical percolation is possible if the inclined fractures intersect both the bottom and top boundaries, and becomes more likely as intensity increases, with all networks that have P 32 > 0.4 percolating in the z direction. The lack of intersections between horizontal fractures lead to the lowest permeability magnitudes at high fracture intensities among SDFN Sets 1, 2, and the GDFNs. This suggests that the intersections between almost coplanar fractures, arising in GDFNs due to fracture reorientation, are important features which are missing from stochastic networks. Such intersections can be seen on the cut planes through the geomechanical network in Thomas et al. (2020b). (a-f) Graphs comparing P 32 against directional permeability (k xx , k yy , and k zz ) in GDFNs A-F (crosses connected by lines) and Set 1 SDFNs (dots). Dots, crosses, and lines are colored by the directional permeability they correspond to, with the legend shown in (a). (g-j) Geometry of example SDFNs from Set 1 with four different fracture intensities. Figure 6. (a-f) Graphs comparing P 32 against directional permeability (k xx , k yy , and k zz ) in GDFNs A-F (crosses connected by lines) and Set 2 SDFNs (dots). Dots, crosses, and lines are colored by the directional permeability they correspond to, with the legend shown in (a). (g-j) Geometry of example SDFNs from Set 2 with four different fracture intensities. The directional permeability results in Figures 5 and 6 highlight several key hydraulic properties arising in GDFNs. The hydraulics of the present GDFNs can be traced to the material properties of the host rock, the initial fracture distribution, and the stress condition driving growth, which result in the domain becoming strongly anisotropic and connected at the final growth step. In GDFNs containing the same number of initial fractures but different distributions (i.e., A and B, C and D, or E and F), the hydraulic properties follow a similar trend, but still have different percolation intensities and final maximum permeabilities. This highlights the variation which arises due to the random initial fracture distribution. It must be reiterated that GDFN hydraulic properties depend on many factors, including the initial fractures (radii, orientations, and locations), material properties, and boundary conditions. In particular, the inclusion of vertical fractures in the present GDFNs enables rapid percolation due to the additional connectivity they provide to the fracture network. Analyses of fracture networks grown according to geomechanical considerations provide opportunities to investigate how initial conditions, material properties, and stress conditions driving growth interact to form networks with different hydraulic properties.
SDFN set two demonstrates that stochastic networks can be tailored to match the anisotropy which emerges in GDFNs. This enables them to better reproduce the anisotropy and percolation properties of a mechanically grown network compared to an isotropic network (i.e., SDFN Set 1). However, the lack of mechanical influences (e.g., fracture interaction, curvature, and selective growth based on geometry) leads to different  Figure 7) and uniform apertures in the right column (b, d, and f). Flow is from the −x boundary to the +x boundary. Fractures are partially transparent surfaces, and the flow velocity at each element is placed at its center. Vectors are scaled in size and color by their magnitude; hence, flow occurs in the matrix but at a much lower velocity than in the fractures. permeability estimates at different fracture intensity ranges, and high variations in vertical permeability for different realizations.

Mechanically Derived Apertures
In this section, flow and k results are presented in GDFNs where fracture apertures are induced by a tensile stress, rather than being uniform. This adds further realism to the analysis of the generated geomechanical networks, where the influence of the stress boundary condition, accumulated growth, fracture interaction, and connectivity can be studied together. In the discretized domain, fractures are represented by triangular surface elements arranged as two coincident surfaces which connect at the tip (see Figure 1). Once the domain is deformed, the discontinuity between two initially coincident nodes can be used to find the slip (displacement along the fracture surface) and aperture (displacement normal to the fracture surface) locally on a fracture's surface. This aperture is assigned to each element as the average of its three corner nodes, and a local permeability is calculated using equation (9). Figure 7 shows the spatial distribution of the apertures on fracture surfaces in GDFN E at growth Step 18. A tensile stress boundary condition of 1 × 10 7 Pa at the top boundary drives the opening of fractures. No fluid pressures are considered during this deformation. The highest apertures are found at the center of flat, large fractures, toward the sides of the domain. At the domain edges, fractures are less densely packed, so less stress shielding takes place here, allowing fractures to acquire large apertures. The traction-free condition on the side boundaries also contributes to higher apertures at the domain edges. There are large variations in aperture magnitude, with most surfaces having apertures in the range of 10 −5 to 10 −3 m. Inclined fractures have very low apertures because they are parallel to the direction of tension, a factor which strongly influences the permeability of the mechanically stressed GDFN during its growth.
For three different tensile stress magnitudes (1, 0.5, and 0.25 × 10 7 Pa), flow velocity and k were calculated in GDFN E for multiple different growth steps. Figure 8 compares the fluid velocity in each element of the domain when apertures are mechanically generated (A, C, and E) and uniform (B, D, and F), viewing the domain in three different orientations. Fractures are visualized as transparent planes, and opaque-colored arrows show the direction and speed of flow at each element. Flow occurs from the −x boundary to the +x boundary, and flow vectors are scaled by their magnitude to highlight flow in the fractures. Flow occurs in the matrix at much lower magnitudes. The mechanical and uniform aperture cases have different flow magnitudes due to the large apertures which result from the stress conditions, therefore the difference between the spatial arrangement of the flow between the two cases is the main focus of the figure.
Figures 8a, 8c, and 8e show that mechanical apertures cause flow to be intensely channeled on relatively few fractures in the middle of the domain. This effect is particularly visible in Figure 8e, where the arrangement of red arrows indicates a fast flowing path. This is the best connected path in the domain and also includes fractures having the largest apertures. Relatively little flow occurs on most other fractures, particularly the inclined fractures, because they acquire very low apertures. The location of flow channels changes significantly between different growth steps and is consistent with the results of Bisdom et al. (2015), where, for networks with stress dependent apertures, only a small fraction of the network was found to contribute to permeability. Furthermore, these results match observations of deep subsurface networks, where a minority of fractures provide flow pathways (Follin et al., 2014;Tsang & Neretnieks, 1998). When apertures are uniform (Figures 8b, 8d, and 8f), flow occurs at similar magnitudes on the majority fractures in the domain, without any particular preference for one path.
Equivalent permeability tensors under three tensile stress magnitudes (1, 0.5, and 0.25 × 10 7 Pa) are compared to results with uniform apertures in Figure 9. Figure 9a shows a graph comparing P 32 against the total void space, , in the domain, contributed by fractures due to their apertures. With uniform apertures, has a linear relationship between P 32 and . When apertures are mechanically induced, the relationship becomes nonlinear, characterized by an initially lower total void space, and an increasing gradient with higher P 32 values. In later steps, the total area of fracture surfaces added gradually reduces (as observed in Figures 3a  and 3d). However, the total void space against P 32 curve follows a different trend, increasing in gradient rather than decreasing. The lowest applied stress (0.25 × 10 7 Pa) provides the closest match to the uniform aperture case.
The effects of mechanical apertures on permeability for three tensile loads are shown in Figures 9b and 9c, which plot k xx against and P 32 , respectively. Figure 9b shows that linear increases in can lead to order of magnitude increases in k xx . When apertures are uniform, percolation is characterized by large increase in k xx in one step. Under stress, this increase is spread over multiple growth steps. This indicates that the domain has become intersected, but low apertures on the inclined fractures cause parts of the flow paths to have very low permeability, spreading the permeability increase due to percolation over multiple growth steps. This particularly clear for the tensile load of 0.25 × 10 7 Pa, where the percolation jump occurs two steps later than in the uniform aperture case. Figure 9c shows that mechanical apertures significantly modify the relationship between fracture intensity and directional permeability. After a domain percolates, permeability continues to increase in the mechanical aperture cases, because the stress conditions provide an additional way to increase permeability, in addition to further growth and linkage. Similar trends can be observed for k yy and k zz .

Discussion
The geometric and hydraulic properties of the GDFNs are broadly similar between the six data sets explored in this work. This suggests that the formation of fracture networks is primarily the result of the regional stress, fracture toughness, and the mechanical process of growth. The stochastic influences on a GDFN's development are the locations and orientations of its initial flaws. The particular arrangement of initial flaws changes the number of steps it takes for a network to finish growing, indicating that the point at which no growth occurs depends on the space that fractures have available for growth. The properties of the initial fracture distribution (e.g., initial sizes and orientations relative to the applied stresses) also have a significant influence on GDFN equivalent permeability. The fluctuation of the vertical equivalent permeability during network growth implies that, even though percolation in the z direction does not occur, the permeabilities in different directions are influenced by one another.
During equivalent permeability analysis, distinct flow channels do not form if uniform apertures are assumed. The primary factor influencing percolation is the linkage of growing fractures between one another and the boundaries. In contrast, using mechanically induced apertures results in highly channeled flow and gradual transitions at percolation densities in GDFNs. This gradual percolation transition is the result of two factors: the change in the relationship between fracture sizes and apertures, and the absence of aperture on the inclined fractures. In a tensile stress regime, the largest apertures occur at fracture centers, meaning apertures on initially horizontal fractures continually increase with growth step in Figure 9. This 10.1029/2019JB018899 increases the importance of fracture size on flow magnitudes, and means that there is no characteristic fracture permeability toward which the equivalent permeability increases to if the domain becomes connected. The absence of apertures on inclined fractures reduces the hydraulic connectivity between fractures. These factors result in a gradual increase in the fracture intensity-equivalent permeability curves when GDFNs become fully connected between boundaries.
Altering the material properties will significantly influence the final geometry of a GDFN. In particular, properties like the mechanical boundary conditions, the maximum fracture extension per growth step, and the fracture growth rate exponents are known to change fracture interaction patterns (Thomas et al., 2020b). These parameters have been shown to influence the percolation threshold in the present GDFNs, suggesting that geomechanical network hydraulic properties are not as tied to fracture intensity as they are influenced by the growth conditions, which are more difficult to quantify. Domain boundaries and fracture termination upon intersection significantly influence the final domain geometry by limiting the extent of fracture growth. The influence of these factors can be investigated by modeling through-going fractures in detail. Furthermore, fracture interaction patterns will be sensitive to the magnitude and orientation of the regional stress (Pollard & Aydin, 1988), whereas only one stress state is considered here.
Further realism can be incorporated into the process of generating GDFNs within the present numerical approach. In nature, fractures which are oriented perpendicular to the most tensile stress in a given deformational event will be most likely to propagate, with some influence from fracture size and interaction modifying propagation patterns. Under the boundary conditions we apply, a strongly anisotropic network results, similar to a single fracture set within a natural network, which formed due to a specific tectonic stress episode. By applying multiple different boundary conditions in sequence, more fracture sets will arise within one domain, each one associated with different stress conditions. Additionally, incorporating fracture twisting into the growth method will improve understanding of how tip fragmentation and breakdown influences network properties, particularly if growth is caused by shear or mixed mode deformation (Martel & Boger, 1998).
The most direct analogs to the present networks may be deep fracture networks in crystalline rocks, which have fracture intensities of the same order of magnitude as the final steps of the GDFNs (e.g., Cvetkovic et al., 2004;Follin et al., 2014). A future step toward generating realistic representations of the subsurface will be to directly match the properties of a GDFN to a particular subsurface environment. This can be approached in several ways. A GDFN can be matched statistically to a present network, for example, by matching fracture intensity of the whole network or individual fracture sets. In future work, the tectonic stress regimes which formed a subsurface network can be used as sequential boundary conditions to grow fractures. Specific geometric features can be investigated and incorporated into a GDFN, such as the interaction between a previous fracture set and an actively growing set. A GDFN that matches a specific subsurface setting may reveal insights into its mechanical and hydraulic properties.

Conclusions
A comparison of DFNs modeled with and without mechanical considerations reveal some stark differences. Geomechanically grown DFNs contain realistic features which are absent from standard stochastic networks. During growth, fractures orient themselves with respect to the principal stress direction, leading to significant anisotropy in the final network. Fracture sizes are highly heterogeneous, and their shapes are nonplanar due the mechanical process of mutual fracture growth, interaction, and intersection. Despite GDFNs growing from an initially random set of initial flaws, all six GDFNs converge toward similar final fracture intensities, with many geometric similarities between the final domains. This demonstrates that the primary influences on the network's properties are the stress conditions, fracture interaction, and fracture toughness, rather than stochastic processes that initiated them.
Each of the GDFNs has a similar pattern of permeability during their growth, characterized by an initially isotropic equivalent permeability tensor, percolation perpendicular to the stress direction after several growth steps, and increasing permeability toward a similar limit, influenced by the final network density. Compared to three sets of stochastic DFNs (SDFNs) with different flaw insertion schemes, GDFNs are much more anisotropic, more permeable at higher densities, and generally percolate at higher fracture intensities. A group of SDFNs which were set up to mimic the properties of the grown GDFNs more closely match GDFN properties, but lack the important realistic features which emerge during the growth process.