Development of Fracture Networks Through Hydraulic Fracture Growth in Naturally Fractured Reservoirs

A 2-D numerical study was carried out, using a fully coupled rock deformation and fluid flow hydraulic fracturing model, on fracture network formation by advancing, widening and interconnecting discrete natural fractures in a low-permeability rock, some of which are small enough to be considered as a flaw that acts as a fracture seed. The model also includes fractures connecting into one another to form a single hydraulic fracture. In contrast to previous fracture network models, fracture extension and fluid flow behavior, frictional slip, and fracture interaction are all explicitly addressed in this model. Incompressible Newtonian fluid is injected at a constant total rate into fractures to study viscous fluid effects on the network formation. The algorithm for flow division and coalescence is validated through some examples. Numerical results show that the incremental crack propagation that connects isolated natural fracture sets depends on the current stress state and the fracture arrangement. The newly created connecting fracture segments increase local conductivity since they are oriented along a path that is easier to open when pressurized by fluid and provide a new path for fluid flow. However the hydraulic fracture growth process is retarded by some of the resulting geometric changes such as intersections and offsets, and the growth-induced sliding that can impose a barrier to further fracture growth and fluid flow into parts of the network. Such barriers may eventually result in a fracture branch initiating and growing that results in a relatively shorter and more conductive path through a fracture network zone. We consider a specific fracture arrangement consisting of around 20 conductive pre-existing fractures to study the effective behavior of the hydraulic fracture growth through a natural fracture network. Mechanical responses have been studied for two different fracture and flow scenarios depending on the fluid entry details: one fracture system assumes each of four entry fractures has one quarter of the total injection rate and the other system is defined to maintain the influx rates into each inlet fracture so that the pressure across all four inlet fractures is equal (but not necessarily constant in time). For the latter case, a preferential flow pathway is developed as a result of hydraulic fracture growth and the overall permeability of the fracture system increases rapidly after this hydraulic fracture path develops. The former injection condition results in development of more evenly distributed advancing fractures that provide a more homogeneous flow pattern.


Introduction
The understanding of fluid movement through fractured rock masses is essential to improving the success of reservoir stimulations for energy resources such as shale gas and geothermal energy and for stimulation by hydraulic fracturing of any naturally fractured rock mass. As we know, fractures play an important role in the flow of fluid through rock masses by building connected networks that channel flow. These networks develop either through enhanced conductivity of existing fractures or by new fracture growth that connects existing fractures. There are many studies in the literature devoted to characterising the fracture system connectivity in relation to fracture orientation, size and conductivity (see the comprehensive review [1]). Besides these static factors, the time dependent evolution of fluid pressure and stress states can generate different fracture patterns that act to assist or inhibit further fluid flow, in particular forming a preferential fluid pathway in the presence of viscous fluid flow [2,3]. Under some circumstances, the fracture growth driven by pressurized fluids can propogate a hydraulic fracture to connect two isolated fracture clusters. Cross-cutting fractures in the connected region are filled with fluid and pressurized but may open or not. Therefore fluid movement in a fractured rock mass will involve both new fracture growth and permeability enhancement of existing fractures. Clearly, using an equivalent porous continuum model to represent fluid flow and fracture growth would be inaccurate, especially when the flow is dominated by hydraulic fracture processes. The approach here is to study the full coupled process in order to determine the parameters that control fracture development. Simplifications and averaging methods can eventually be realistically employed without degrading the ability to predict fracture growth.
Discrete fracture models, so named because these models treat fractures as discrete entities, are applicable whenever the process involves fractures growth and flow where details such as opening, shearing or growth of the fractures are being studied. If many fractures are considered, these discrete models become computationally demanding. Such a system is very heterogeneous and localized in both fracture growth and fluid flow. Early numerical models treated the hydraulic fractures as single planar fractures and did not consider fracture interaction [4][5][6]. The emergent behaviours associated with fracture propagation under a tensile displacement boundary condition have been described as straight paths using a subcritical failure criterion and a propagation speed exponent. Recently, the effect of curving fractures on fluid flow in both the fractures and the matrix under a tensile displacement boundary condition has been considered [7,8]. However, the pressure distributions used inside the fractures are either uniform [7] or based on the steady-state solution for flow in a porous medium [8]. The effect of fluid viscosity on pressure distribution has therefore not been considered in these fracture network models. Moreover, some simplifications in the fracture geometry changes, such as the details of interactions at intersections and in regions where fractures close, are used in these models and these simplications affect fluid flux distribution. Uncoupling of deformation and fluid flow, as used in some models, limits the application of the results obtained.
Fracture models that do not consider flow in the matrix, are applicable to low-permeability rocks [9]. In general, the fracture aperture is very narrow and would be different in each fracture branch intersecting at the same point, and these differences are expected to be significant in their effect on fluid movement through the intersections. In contrast to the stressdriven uniform fracture nucleation in porous rocks subject to a tensile loading environment, the propagation of a hydraulic fracture through a network of pre-exsiting fractures is dependent on local stress states around fracture tips, which in turn depends on the fluid flow and pressure along the non-planar hydraulic fracture path. Under most circumstances, one dominant hydraulic fracture is generated and the entire injected fluid rate is carried to the outlet through this preferential flow path, while most shorter natural fractures that are intersected by this main fracture remain closed or act as dead end branches [9]. Models that include the coupling of rock deformation and viscous fluid flow provide a means of studying the fracture development and the evolution of distinct preferential flow paths and development of a dominant fracture.
At intersections of two or more fractures, the kinematic deformation transfer between slip and opening of the fractures can induce additional fracture aperture changes. In addition, the viscous effect becomes stronger for narrow channel widths, which are commonly associated with intersections and offsets. The importance of tracking the details of fracture geometry lies in the fact that although the pressure losses may diminish after a long time, initial fracture geometry details may strongly affect the final fracture patterns. Studies that neglect viscous fluid effects, by using uniform pressure or steady-state transport and deformation models, will, in many cases miscalculate the stresses and flow rates, thus producing incorrect fracture and flow patterns as time-dependent pressure responses are not determined accurately.
Mechanical interaction among fractures has not received sufficient attention in the literature involving discrete fracture models especially for cases involving pressurized viscous fluids sufficient to result in hydraulic fracture growth. Any inaccuracy in the calculated fracture pathway may cause incorrect flux redistribution at intersections, as fluid flow behaviour is strongly dependent on local width. Due to intrinsic complexity of the problem, numerical methods appear to be the only approach able to explicitly solve the nonlinear and nonlocal fracture-fluid and fracture-fracture interaction in such fracture network. The Distinct Element Method (DEM) and the Finite Element Method (FEM) have been used for this purpose [10,11]. However, the fracture pathways are confined along the element edges in classical FEM models and the out-of-plane fracture propagation is difficult to accurately simulated because it requires remeshing. We have developed a Boundary Element Method (BEM) based program for treating this coupled problem [12,13]. The validation of the code has been carried out for various simple cases involving both viscous fluid and uniform pressure. Additionally, the program treats rock deformation and fluid flow as a whole in that the field variables are obtained in a single framework, instead of the one-way coupling scheme as used in Reference [11].
The purpose of this paper is to provide some initial results for hydraulic fracture growth through a natural fracture network. Fracture growth is allowed and is based on a local failure criterion [14] and fracture coalescence can take place to form a path through an existing network of natural fractures. The numerical treatment of fracture coalescence has been detailed in [12,13]. The rock mass is assumed to be impermeable and the fluid flow is confined to occur along the pre-existing or newly created fractures. The fracture nucleation sites are embedded as the pre-existing secondary fractures with small sizes to reflect the tensile strength heterogeneities existing along the fracture surface [15]. For this plane-strain model, the strain in the out-ofplane direction is assumed to be zero and the fracture should be visualised as extending uniformly a significant distance in this direction.

The model
The basic governing equations and the boundary conditions are provided in our previous work (see References [12,13,16]), dealing with the hydraulic fracture model that fully couples mechanisms of rock deformation and viscous fluid flow. The basis of the model is briefly described here for the sake of completeness.
We allow for the fracture surfaces to be rough and tortuous, which imparts a hydraulic aperture for the closed fracture allowing fluid flow, but causing no stress and deformation changes. Fluid volumetric flux q in a closed fracture segment is described by: and in an opened fracture portion it is where μ ′ = 12μ with μ being fluid dynamic viscocity. w and ϖ are mechanical opening and hydraulic aperture along the fracture surface and are functions of time and location. The former is determined by the stress condition given below, but the latter obeys an evolution equation linearly proportional to fluid pressure change [13]. It is noted that the initial value of hydraulic aperture is denoted as ϖ 0 , which is a reflection of the fracture surface roughness and tortuosity.
Fluid flow in the opened fracture portion is based on the lubrication equation: but fluid flow in the closed fracture segment is based on the pressure diffusion equation where χ 1 is the compressibility of the fracture with units of Pa -1 and it is set as 10 −8 Pa -1 in this paper. When fracture surfaces are separated, the change of the hydraulic aperture ceases, but its contribution to fluid flow is retained as provided in Eq. (3).
The nonlocal elastic equilibrium equations for a system of N fractures are given as: where the coordinates of a location of a material point are denoted as x = (x, y) in a twodimensional Cartesian reference framework, t is time and v is the shear displacement discontinuity, ℓ r is the length of fracture r. σ n is the normal stress and τ s is the shear stress carried by the fracture because of its frictional strength, obeying Coulomb's frictional law characterized by the coefficient of friction λ, which limits the shear stress | τ s | ≤ λσ n , that can act in parts of fractures that are in contact, but vanishes along the separated parts. Along the opened fracture portions, we have σ n = p f .
In addition, σ 1 and τ 1 are the normal and shear stresses, respectively, along the fracture direction at location x caused by the far-field stress, whose normal and shear components are denoted as σ xx ∞ , σ yy ∞ and σ xy ∞ . G ij are the hypersingular Green's functions, which are proportional to the plane strain Young's modulus, The global mass balance requires And the fluid front in the hydraulic fracture will be, in general, not coincident with the fracture tip. The fluid front location is found by using Eqs. (1) and (2) with the following equation The fracture growth is based on using the maximum hoop stress criterion, with the maximum mixed-mode stress intensity factor reaching a critical value [14] 2 3 cos ( cos sin ) 2 2 2 where K I and K II are calculated stress intensity factors, K Ic is tensile mode fracture toughness and Θ is the fracture propagation direction relative to the current fracture orientation. The predicted orientation follows the maximum tensile stress direction, and the near-tip stresses are approximated by the analytical LEFM solutions [14].
The problem must be completed by specifying the imposed boundary conditions at the wellbore, that is, the sum of injection rates of hydraulic fractures connected into the wellbore or to the entry zone, should be equal to the given injection rate Q 0 . At the fracture tip the displacement discontinuities are zero, w(ℓ r , t) = 0 and v(ℓ r , t) = 0. In addition, the entire system is assumed to initially be stationary and unsaturated.
The numerical scheme for the above nonlinear and nonlocal coupled problem has been detailed in our previous paper [16] based on the Displacement Discontinuity Boundary Element Method. The hydraulic fractures and other geological discontinuities like joints and faults, and natural fractures are discretised with constant displacement elements. The model solves the hydraulic fracture problem simultaneously including the effects of viscous fluid flow and coupled rock deformation. The solutions are consistent with existing results. Also, a fracture can intersect another one in its path and a new fracture can be nucleated from a position on the natural fracture. In particular, the new fracture seeds are pre-defined in this paper along some natural fractures. The interested reader is referred to our previous papers for the details on implementation of fracture growth and coalescence [12,16]. One important check here is the satisfication of fluid mass conservation after fluid branch coalescence, due to redirection of newly-created fractures.

A simple test problem
To validate the model on coalescence of fluid flow branches, we present results for one specific case, as shown in Figure 1, where two hydraulic fractures with a spacing 0.8 m driven by the same fluid source (located at the plane with x=-1.0 m on the left of the natural fracture) intersect the natural fracture orthogonally. And the sum of their inflow rates is a constant. A new fracture site is located on the natural fracture that is 2.4 m in length, and it is assumed to be a pre-existing flaw at 0.12 m long. It is anticipated that after intersection, the fluid will invade the natural fracture until it reaches the new fracture in Figure 1(a), which is subject to a less compressive stress. As the new fracture becomes the weakest point to continue crack growth, it will propagate when the fluid pressure reaches a sufficient value, as shown in Figure 1(c). It should be noted that when the middle section of the natural fracture between two hydraulic fractures is filled with fluid, a sudden increase in the pressure occurs as indicated in Figure  1(a) and (b) for two close time steps, in the same way as injection into a closed container. At the time shown in Figure 1(c) the natural fracture appears to be full of fluid, the whole fracture system experiences a similar pressure level, high enough to cause the new fracture to propagate. The material constants used for this problem are listed in Table. 1 if not otherwise specified.
The variations of flow rates into each branch with time are provided in Figure 1 to the fracture nucleation site is contributing more in fluid flux to sustain the new fracture growth. As for the loss of geometric symmetry, it is interesting to note that the outflux from the top hydraulic fracture is also larger than its counterpart since Q 1 > Q 2 clearly shown in Figure 1(d).

Random fracture geometry
In this section, we present more complex cases where several high-angle joints are defined in a random distribution ahead of hydraulic fractures that grow from a left entry zone (along the plane x=-1 m) toward the right (through the plane x=2 m). The remote stress conditions and some geometric parameters are provided in the caption of Figure 2. The fracture injection occurs into four individual fractures located on the far left, as shown in Figure 2. The fracture segments coloured green denote the initial existing fracture configuration. Each existing fracture consists of a single fracture or of several connected fractures of different sizes and orientations. The hydraulic fractures driven by pressure act to connect these separate fractures to form a conductive path from left to right and the newly created fractures are coloured in red in Figure 2. Two different injection boundary conditions are used in the computations. One is associated with even distribution of the injection rate into the four entry fractures on the left, and the second condition specifies that the pressure at each of these four initial fractures is equal. The latter condition is physically reasonable if we assume a wellbore lies along the yaxis of Figure 2, while the former condition would require isolation and injection into each fracture at the same rate.   Of course, some numerical difficulties exist in the simulation of the process of connecting two intersecting fractures [16]. A mesh sensitivity analysis has been carried out so that the change in fracture orientation does not affect the numerical accuracy as the fracture direction has to be such that the intersecton occurs at the the end of an element on the fracture. A fracture connection event that results in a strong deflection from the calculated direction can influence the subsequent fluid flow. In these cases, the natural fractures must be re-meshed and the job must be run again. In the numerical simulation, the individual element size in the vicinity of intersection point is chosen with care, so that the fracture connection can be completed as a smooth path.
For the sake of convenience in discussion, the above two injection types are respectively named as type I and type II injection. Figure 2 shows the results at the time of breakthrough when the hydraulic fracture emerges on the right side (x=2 m) of the network of natural fractures, for these two injection types. After fracture reorientation to the direction normal to the minimum principal stess, the fracture development through the network zone is complete. As expected, the type I injection condition results in more fracture growth paths and intersections across the network, while the type II condition results in a localized path, with only one hydraulic fracture continuing to grow past the network on the right side. In Figure 2(b) there are two separate unconnected sets of fractures from top to bottom, both of which connect to two of the initial fractures on the left which are connected to (in this case) the equal pressure fluid source. The influxes into the bottom fracture set decrease in time and more of the total volume injected enters the upper fracture set. The upper fracture set is more conductive. We define the plane x>2 m as the exit zone. There is one outlet fracture in this case. However, there are two fluid outlets or extraction sites in Figure 2(a) for type I injection. On the other hand, for the type II injection, the new fracture segments are mainly created along the upper fracture path. Figure 3 shows the outlet flux variations in time for the two cases used in Figure 2. For type I injection, only around 43 percent of injected fluid has passed through the fracture system to the outlet at the end of simulation period and the other 57 percent is contained in fracture branch inflation or growth in the network. It is also found in Figure 3 that in this case the outlet flux (sum of outlet 1 and outlet 2 for case (a)) increases at a very slow rate at the large time.
More fractures under type I injection are connected with each other as a result of crack growth near the fluid source and some fractures continue to grow after the main conductive channels are developed. Fracture segments not opened also store some fluid because of their initial conductivity. The rate that fluid is stored in the fracture network differs between the type I and type II injection conditions. For type II injection, more than 55 percent of injected fluid exits the outlet fracture at the end of the simulation and this outlet flux is increasing at a very large rate as shown in Figure 3 so that the rate of fluid volume stored in the system will become extremely small soon in light of this trend.
In addition, the breakthrough time for these two injection types is provided in this figure. The type II injection has a much earlier breakthrough time so that the fracture growth through the fracture network is more rapid. The rapidly increasing trend in outlet flux is also reflected in the opening profiles as shown in Figure 4(a) with the wider open fracture path corresponding to the path that carries the most fluid. The wider fracture channels are localised along the preferential pathway along the upper fracture set to the right-top fracture outlet under type II injection. It is found that at the given time instant, up to 85 percent of the injected fluid is pumped into this preferential path and the fractures not included in this path are all static at late times in the simulation. At the larger times simulated, most fluid just passes through this highly conductive channel to reach the exit zone. Thus, the localised flow channel provides the lowest resistance pathway for fluid flow. Here, we note the contribution of residual hydraulic aperture on the fluid movement and storage. As stated above, each natural fracture has a pre-existing aperture of 0.01 mm in the computations. Thus, some fluid enters and is stored in this pre-existing aperture.
It is found that fractures with wider opening have had their opening enhanced by the large slip along the longest oblique natural fracture, which is oriented at 45 degrees to the x-axis, as shown in Figure 4(b). The kinematic transfer between slip and opening assists the fracture opening.   Also, it is found that some fracture connections are made even after the preferential fluid pathway is established, since fracture growth can still be generated in the higher pressure region near the entry. After the breakthrough, the injection pressure has to be retained at a higher level, to force fluid through width restrictions that occur at intersections and offsets, as shown in Figure 5. The continuing fracture growth in the network can alter the earlier opening distributions. Comparing Figure 4 and Figure 2(b), it is clear that fracture growth and fracture interconnection can occur inside the network region after the flow breakthrough. Of course, one reason for pressure increasing after breakthrough is attributed to the strong resistance to fluid flow at the last vertical natural fracture on the connected flow path for this specific geometry. This vertical fracture provides the strongest barrier to fluid flow as its opening is highly constrained by the geometry, with the approaching hydraulic fracture making a right angle to the natural fracture. A new fracture would possibly be nucleated at some location along this vertical natural fracture which would reduce this restriction. Although Figure 5 shows the increasing trend of the injection pressure, the pressure will eventually level off or even decrease as the outflux from the network increases to 100 percent of the input rate, as indicated in Figure 3.

Fluid-driven fracture nucleation, growth and connection
In the model, we only consider hydraulic fracture growth through a finite set of fractures, some of which are initially very small, but potentially provide a conduit with the help of high fluid pressure. The fracture seeds are pre-assumed in this paper and are represented by these small fractures. Therefore, fracture nucleation in a highly stressed area is not dealt with in this paper. This treatment of fracture nucleation can underestimate the fracture number and the fracture connectivity. For the case shown in Figure 4, the lower fluid flow path cannot move to the right outlet and a higher stress level might create new fractures near the entry zone. It would be interesting to study the impact of crack nucleation based on stress conditions rather than just from these pre-existing fractures.
Without considering fluid loss into the rock matrix, hydraulic fracture growth as the main driving force in connecting fractures can create more new segments at the upstream end of the fracture system than at the downstream. In terms of fracture length, both pre-existing and newly created, we can define fracture density as the fracture number per unit area. Although fracture density is a significant measure for fracture connectivity, the longer fractures including the newly created parts would be more important contributors because they are more compliant and will open wider under the same internal fluid pressure [9]. Predicting the early growth of these hydraulic fractures through a pre-exising network, as modeled in this paper, must account for the effect of viscous fluid or incorrect fracture behaviour is predicted.
In this model, fracture growth occurs when the failure criterion is satisfied at any fracture tip, with the failure condition defined within the framework of linear elastic fracture mechanics. Fracture curving is the natural result of the local stress field around the tip if the growth follows the maximum tensile stress criterion. Normally, fractures will reorient themselves to the maximum compressive stress direction to increase fracture opening, resulting in local conductivity enhancement as indicated in above results. However, the fracture curving can sometimes lead to intersection of two fractures at a small acute angle, which will make it difficult for the subsequent flow to enter some segments. Sometimes, the subsequently developed sliding on one fracture can seal the fracture channels near the junctions. The development of geometric networks, produced by growing hydraulic fractures, are illustrated by the results obtained above. The results imply that not all connected fractures can contribute to overall conductivity of the system which is contrary to conventional percolation model predictions. These geometric factors affecting fracture growth and fluid flow have been mentioned in early studies [7]. Some fracture growth can occur in the wake of the fracture and flow fronts near the higher pressure entry zone. Local reversed flow has also been observed in the results due to the pressure changes.
Actually, in addition to injection conditions, many other factors such as injection rate and in situ stress can affect the crack growth and coalescence. At the elevated pressure and based on the assumed fracture geometries, one can find that, even through a network of natural fractures, the hydraulic fracture average direction tends to align as much as possible with the direction of the maximum stress. This orientation reduces the viscous dissipation and injection pressure. However, an increase in fluid pressure because of increased rate or viscosity can produce higher pressures upstream of a local offset or restriction which can then lead to opening of cross-cutting natural fractures and branching.
In addition, this paper only considers a limited number of specific initial fracture geometries, the results may be different if the starting geometries are changed and more cases are being considered as a way of making our conclusions stronger and more general. Although a method to deal with network development is presented here, there is a need to work on different geometries to extract some useful general responses for rational simplifications of the expected response for hydraulic fracture network growth.

Implications for fracture-controlled flow system
Our numerical results quantify the overall path of discrete hydraulic fractures growing through a network of pre-existing natural fractures, as shown in Figs. 3 and 4. These results give insight to the behaviour of fracture-controlled flow systems, where the fluid flow and the rock deformation and fracturing are strongly coupled. It is clear that the conductivity depends on the stress-dependent fracture aperture through a strong coupling to fluid flow, as opposed to fixed aperture fractures in conventional percolation models. Local areas can exhibit higher effective permeabilities or strong growth barriers. Such enhanced or restricted opening occur at intersections and offsets, and their existence can affect the total system conductivity, producing a higher pressure level as shown in Figure 5. Early time rapid hydraulic fracture propagation and intersection of small natural fractures establishes a path for the fracture through the natural fracture network, and a single fracture connection event can cause a strong change in the hydraulic fracture channel system that develops.
The model has particular application for understanding the hydraulic fracture connection process through a network. By varying parameters, one finds the transition of fracturecontrolled flow pattern from more uniform to more localized and from multi-directional to unidirectional. The hydraulic fractures tend to develop wider and more connective localized fracture channels in establishing a preferred path through the network of natural fractures. In contrast, low rate injection processes that do not involve significant fluid viscous dissipation effects, tend to result in flow occurring along all already connected conductive paths. How to better characterise this difference is still open and to find meaningful parameters in connecting the intricate topological fracture network with diffusion flow patterns requires prediction of propped and unpropped fracture permeability that remains after the hydraulic fracture treatment. The model used here may provide a tool for such parametric studies.