Reactivation of natural fractures network by hydraulic fracturing

. Modelling the process of induced fracture initiation, propagation and interaction with natural fractures is a very challenged task. Significant progress has been made in recent years in the development of complex fracture models to address the needs for more suitable design tools than the conventional planar fracture models. However, some aspects of this complex process are not still fully understood in terms of their impact or importance to the overall fracture geometry, or the complexity of simulating them is beyond the current modelling capabilities or requires computation effort that is not practical for engineering purpose. A technique that has been developed to represent the process of fracture propagation is the Continuous Approximation of Strong Discontinuities, which introduces a special kinematics, capable of representing the process of degradation of the material. One way to implement this approximation is to introduce the effects of a very narrow band of localized deformations within the existing finite elements. In this paper was used a finite element procedure that performs numerical analysis of fluid flow in a deformable porous media in a fully coupled scheme. In this analysis the propagation of the hydraulic fracture occurs specially along the pathway of the natural fractures, due to their lower tensile strength and greater permeability.


Introduction
Hydraulic fracturing is an evolving technology that has been used to increase oil recovery in naturally fractured reservoirs [1], usually all reservoirs are fractured to some degree. Depending on the density, dimension, orientation and cementation of natural fractures and the location where the hydraulic fracturing is executed, pre-existing natural fractures can impact hydraulic fracture propagation and the associated flow capacity [2].
The network pattern of fractures resulting from a hydraulic fracturing operation depends on several factors, among them, the strength, deformability and permeability of the rock mass and the in-situ stress state [3], because natural fractures are formed within geological time scale and the formation itself may have gone through multiple tectonic events, the distribution and orientation of the natural fractures do not necessarily relate to the current insitu stresses [2].
When a hydraulic fracturing operation is being done in a naturally fractured reservoir it is possible to have two scenarios, on one hand, the stimulated natural fractures, if-well connected, can enlarge the reservoir contact area, expand the fluid flow path and enhance productivity. On the other hand, the reactivated natural fracture and weak planes may divert fluid flow from the main hydraulic fracture channel and lead to an unreasonably high proppant concentration and early screenouts.
There are two competing forces that govern the general trend of hydraulic fracture propagation path: one is the propensity to propagate along the direction of the weak planes resulting from the natural fracture sets, and the other is the tendency to propagate along the direction of maximum principal stress and open against the minimum principal stress. Often, the directions of natural fractures sets are not aligned with in-situ principal stresses, so the fracture propagation has to balance these two forces to find the path that has lesser fracturing energy, which is why the hydraulic fracture propagates along a zig-zag path [2].

Fracture Network Modelling
In this paper a Continuous Approximation of Strong Discontinuities is adopted to model the behavior of natural fractures under a hydraulic fracturing process, this approach introduces a special kinematics, capable of representing the process of material degradation, the localization of deformations in narrow bands (discontinuity in the field of deformations, or weak discontinuity) associated to the zone of fracture processing, until the complete degradation which corresponds to discontinuity in the displacement field (strong discontinuity).
This technique allows the representation of the discontinuity effects inside finite element, providing great advantages in computational cost for this type of modeling, since it eliminates the need to use excessively discrete meshes, to remesh several times the geometry of the problem or to use discrete interface elements.

Finite element with embedded discontinuities
A strong discontinuity approach to embed discontinuities into finite elements is considered to represent the behavior of fractures in rock formation. The mathematical development of the continuous approximation of strong discontinuities introduces the effect of a strain band within existing finite elements. This approach, as proposed by [4], was implemented in CODE BRIGHT [5] [6], a finite element procedure that performs numerical analysis of fluid flow in a deformable porous media in a fully coupled scheme [7].
To properly derive embedded discontinuity finite element formulations, fundamental aspects regarding to the kinematics and statics of the discontinuity must be considered. The kinematic enrichment must correctly reflect the position of the interface in the element as well as the relative displacement (opening and sliding) between the two opposite faces of the interface, Furthermore, the traction continuity condition must be correctly imposed to endure a correct relationship between the tractions in the internal interface and the stresses in the surrounding continuum portion.
The implementation of this technique in CODE BRIGHT has been quite convenient and can be done in non-intrusive way in the code. Since the jump in the field of displacements can be approximated as equivalent inelastic deformations, which can be calculated in the constitutive module of the program in finite elements, within the procedure of integration of stresses in the elements. The standard MEF equations (without discontinuity) are used to represent the continuous part. The behavior of the discontinuity interface can be described by a constitutive equation that directly relates tensions and displacements (discrete model) or by a constitutive model of the continuous type, which relates tensions and deformations, according to the continuous approximation of strong discontinuities [3].
Consider the triangular element with three nodes from domain Ω length l , with a band of strain localization S , width ℎ, which divides the element into two parts, isolating the node 1 from nodes 2 and 3, as shown in Fig.   1. Fig. 1 Finite element crossed by a discontinuity band A displacement field jump, ⟦ ⟧, on lead to a rigidbody relative motion, the strain field of the continuum portion and discontinuity can be written as: The corresponding stress field is given by: Where ℎ denotes the bandwidth of strain localization, Σ is the constitutive matrix of the material; M and n are the matrices obtained from the components of vectors m and n respectively and ⟦ ⟧ is the jump in displacement field. Traction continuity is assured as follows:

Flow Problem
Flow following Darcy's law into the discontinuity can be decomposed in a normal and tangential component as shown in Fig. 2, the formulation adopted in this paper considers that the fluid flow in a fracture occurs on its direction, therefore only the tangential component is taken into account. Darcy's law can be written, for continuum and fracture, as follows: Since a fracture in a porous medium becomes a preferred path for the flow towards the discontinuity, that can be interpreted as anisotropy induced, in the direction of the fracture, throughout the medium, as a result, an effective permeability tensor for the element in given by: Two different permeability variation laws, one for continuum portion and other for fracture, are adopted as shown in Fig. 3.

Fig. 3 Permeability Variation laws for continuum media and fractures
It is important to highlight that the finite element program CODE BRIGHT solves the mechanical and fluid flow problems in a full coupled scheme, where, each Newton-Raphson interaction, a single system of equations is solved for the entire mesh, where unknowns are pressure (hydraulic problem) and displacements (mechanical problem) [3].

Study Case
Was simulated a case with synthetic natural fractures, submitted to hydraulic fracturing operation, in order to assess the influence of the in-situ stress state on induced fractures orientation and also in the density of fracture network.
The rock mass is under as initial stress state as can be seen in Fig. 4. In this case it was considered that there is no fractures propagation in the rock, so the evolution of the degradation is restricted to natural fractures. This degradation occurs when the effective stress state reaches degradation criterion of constitutive damage model.
The simulated case represents a rock mass of 200x200m, cut by different natural fractures. A constant fluid injection is applied to the lower part of the rock mass, at the end of one existing fractures as is show in Fig.  5.  The boundary domain is impermeable and initial pressure is zero, this initial condition represents the case where in-situ stress state is given in effective stress state. Mechanical and hydraulic properties of this case are in Table 1.

Results
In this section, the results of the numerical simulation for the different stress state scenarios will be presented. Figure 6 shows the pressure distribution in two simulation moments, for the isotropic stress state scenario. For this scenario, the fracture opening occurs initially on two direct fractures where fluid injection is applied, then are opened the connected fractures until reaching the upper section limit. Figures 7 and 8 present the pressure distribution for stress anisotropy scenarios 2 and 3 respectively. For these scenarios, the main fracture where the fluid is injected is initially opened and other fractures connected with it are pressurized during the simulation. In the scenario of less anisotropy, the fractures connected to the main fracture propagate in a distributed way, and those that are in the upper limit of the section are influenced, in the two cases of anisotropy the main fracture that is initially opened is in direction of major stress ( ), when this fracture is completely pressurized occurs a stress state reorientation, which results in the propagation of the other fractures.   Figure 9 shows fractures opening during pressurization, this is due to elements degradation that reached the fault, it is possible to observe that there is a greater opening of the main fracture when the stress anisotropy is greater, contributing from remarkable form in global permeability.
The fracture opening presented in Figure 9 is reflected in a permeability increase, as can be seen in Figure 10 which shows the final distribution of permeability after hydraulic fracturing.
The permeability values presented are referred to macro element permeability (effective permeability), since the technique used incorporates the permeability due to the fracture that enriches element permeability in fracture direction according to the formulation presented in this work.
In the scenario of initial anisotropic stresses, the increase in permeability occurs in fractures connected to the main fracture, thus forming a network of connected fractures.

Conclusions
The simulations presented considered that the hydraulic fracturing operation only reactivated and opened the preexisting natural fractures.
The results presented show a sensitivity of hydraulic fracturing in relation to the initial stress state, the numerical techniques and constitutive models used were able to predict the formation of a fracture network.