Sensitivity Analysis of the Influence of Fracturing Spacing in the Construction of Complex Fractures Network for Exploration and Production of Shale Gas/Shale Oil

Shales reservoirs have a high degree of anisotropy due to the presence of natural fractures (NFs) and also the orientation of beddings. Thus, hydraulically induced fractures (HFs) interact with natural fractures and generate a network of fractures with complex geometry. The existence of NFs modifies the stress field in the shale and directly influences the geomechanical behaviour of the HFs during the fracturing operation, generating branches in the dominant fracture and contributing to the complex network of fractures. The construction of a network of fractures increases significantly the conductivity of the formation, as it connects previously isolated fractures and pores, thus increasing the productivity index of the wells and providing greater economic viability in the shale gas/oil reservoir designs. This work presents a sensitivity analysis of the influence of fracturing spacing in the construction of the network of complex fractures generated in shales, aiming to understand how this parameter modifies the volume of stimulated reservoir (SRV) and the distribution of propant in the network of fractures, in order to avoid problems in this step of the design and thus, maintain the economic viability of the network. The literature review includes the main published works on this subject and the unconventional fracture models (UFM) used to model the network of complex fractures. Sensitivity analysis will be performed using the MShale software, which uses a stochastic model of the discrete fracture network (DFN) method and numerically solves the equilibrium equations and pore elasticity for shales in terms of effective stresses, in addition to mass conservation equations, linear momentum and energy with viscous dissipation for stokes creeping flow. For the analysis, the other parameters that influence the construction of the network will be kept constant and only the spacing between fracturing will have variation.


INTRODUCTION
The generation of oil occurs in sedimentary basins such as lakes, oceans, rivers and marshes, where sedimentary rocks with large amounts of fine grains, clay and silt are deposited over millions of years and carrying a certain volume of organic matter accumulated in the interior of their pores [1], Such rocks are called source rocks [2]. After the maturation of the organic matter present in the source rock, a natural fracture occurs in the same one, due to the conditions of high pressure and temperature, generating the primary migration of the hydrocarbon for its later imprisonment in a rock with high permo-porous properties, denominated rock reservoir [3]. The next step of the conventional petroleum system generating mechanism is to surround the reservoir by means of traps and a sealant rock, i.e. with low permeability, e.g. evaporite (salt) or shale itself, thereby preventing the hydrocarbon be drained to other layers [3].
The shale gas /shale oil system is different from the conventional reservoirs (CR) mentioned above, since shales are part of a group known in the literature as unconventional reservoirs (NCR), because the primary migration has not yet occurred and also, because the shales have low permo-porous properties. Thus, the drilling operation is done straight on the source rock, so the shales are classified as source-reservoir rocks (SRR) [4]. The structure of the shales is basically formed by clastic or detrital rocks. Its formation is by means of transport by water or air. Such rocks are usually composed of silica, combined with other common minerals such as feldspar and clay minerals, e.g. quartz, limestone and feldspar, which have low granulometry, making it difficult to explore shale gas /shale oil without use stimulation techniques. Another characteristic of the shales is to possess high total organic content, TOC [5].
In the last decades, the advances in directional/horizontal drilling technologies, through rotary steerable systems (RSS), which allowed for better control of the well trajectory, and also in hydraulic fracture stimulation techniques (Figure 1), allowed economical exploration and shale gas/shale oil reserves in the world [6].
The main challenge of petroleum engineering for the exploration and production of shale gas/oil reservoirs is that their geomechanical properties are very different from conventional reservoirs. A big volume of natural gas is still confined to the source rock as the shale and, therefore, requires a stimulation technique so that its extraction is economically feasible [7]. In general, the main characteristics of shales are: • Large variations in mineralogy [8]; • High anisotropy and heterogeneity, due to the laminations of the depositional process, [9]; • Porosity between 2 % and 15 %, [10]; • Presence of high volume of organic matter [11]; • High degree of brittleness [12]; • Presence of natural fractures [13]; • Permeability between nanodarcy and microdarcy [3].
The main objective of this work is to understand how the fracturing spacings influence the support of the network of complex fractures and SRV in the construction of the network of complex fractures, (Figure 2) in order to optimize the flow of shale gas/shale oil to the well [14].

Problem statement
The stimulation of a certain formation by means of the hydraulic fracturing technique consists in the pumping of a high pressure fracture fluid and controlled flow in a producing zone that has low permo-porous properties, in order to exceed its limit poroelastic and to generate artificial fractures which aim to increase the productivity index of the well by increasing the conductivity and permeability of the payzone. Usually it is the zone where the oil reservoir is located. In shale exploits, the pay zone is the range of the shale itself being traversed. In conventional reservoirs, the HF propagates in a two-wing shape, however, in shale gas/shale oil reservoirs, due to some geomechanical, hydraulic and design parameters, the activation of natural fractures occurs through the HFs, building a network of complex fractures [15,16].
Thus, the hydraulic fracturing operation has the purpose of: • Connecting isolated pores and existing natural fractures and with this, build a network of complex fractures; • Exceeding the radius of damage caused by the skin effect during the drilling operations, due to the invasion of gravels in the formation pores and also the projectiles in the cannon phase.
According to Gandossi [17], the fracturing fluid consists of: Base Fluid + Additives + Sustaining Agent (AS). The construction of a network of fractures increases SRV and with this also occurs the increase of shale gas production, [14]. After the fracture network is built up, the fracture fluid is replaced by a fluid with greater transport properties of AS for its lodging inside the network and to provide support for it [18].
The main parameters that influence the network behavior of complex fractures in shales are: • Fracturability [9]; • Number of existing natural fractures [19]; • State of stresses [20,21]; • Spacing between fractures [19]; • Shale's degree of brittleness [22]; • HF intersection angle with NFs [20,21]; • Flow speed and pumping pressure of the fracturing fluid [23]; • Dynamic viscosity of the fracturing fluid [22]. Normally, the shales have several natural fractures present in the rocky matrix and a high degree of brittleness. In this way, the hydraulic fracturing design aims to build a network of hydraulic fractures connected to the pre-existing natural fractures, generating multiple channels within the formation, to optimize the flow of shale gas to the well through a large increase [24,25].
In a hydraulic fracturing design in shales, the following parameters should be considered: • Type of fracturing fluid, [8]; • Pressure and pumping flow of the fracturing fluid, because the higher the injection flow, the lower the pressure required in the pumping of the fracture fluid to overcome the resistance limit of the rock and generate the fracture [23]; • Spacing between fractures • Dynamic viscosity of the fracturing fluid, because viscosity variations generate different behavior of the network of complex fractures [26]; • Fracture volume, to estimate the volume of propant to be used; • Breakdown pressure [27]; • Permeability and anisotropy of SRR; • Existence of natural fractures in SRR [28]; • Containment stresses [23]; • Controlling the orientation of the fracture, as the propagation of the same in a disordered way can reach adjacent layers, water tables and neighboring wells [29].
In addition to the parameters mentioned above, for a hydraulic fracturing operation in shales to be successful, it is necessary that the rock has great fracturability property, i.e., the capacity of the shale to propagate the dominant fracture with several secondary fractures connected to it [30], in order to obtain the highest reservoir volume value possible.
High fracturability is associated with friable shale regions, with high modulus of elasticity and low Poisson ratio. In these conditions a wide and stable fracture network is created [31,32].

LITERATURE REVIEW
The development of numerical methods that simulate the behavior of the network of complex fractures in shale gas/shale oil reservoirs has been a huge challenge among engineers of stimulation and researchers of the petroleum industry in the last decade. As previously mentioned, the traditional models used in conventional reservoirs do not adequately reproduce the behavior of the fracture network in shales [33]. Thus, unconventional models for the adequate modeling of the problem take into account the effects of very low permeability, high heterogeneity and also anisotropy of the formations.
Main published studies about the subject Singh [34] studied a rock intercepted by a single set of uniformly distributed joints footnote and found the equivalent modulus of elasticity for the whole system.
Meyer [35] developed a solution that takes into account the effects of natural convection heat transfer between the fracture fluid and the rock and by means of the Nusselt number, obtained a coefficient to consider the thermal effect during the operation the hydraulic fracturing.
Cho et al. [36] presented the equivalent constitutive equations that describe the mechanical behavior of the rock and the fluid in its interior and with several orientations of the joints. The model assumes that the behavior of the rock is elastic-linear and with isotropic permeability of the intact rock. The theory obtained can be used to take into account the effect of anisotropy due to the orientation of the joints. Huang et al. [37] proposed a model based on the stress-strain relationship that treats the intact rock mass (with very low permeability) and the fractured mass (with high permeability) separately. Thus, the equations of the modulus of elasticity for the rock with 3 joints are obtained as a function of the fracture properties also of the intact rocks.
Fisher et al. [38] carried out hydraulic fracture diagnosis design in Barnett's naturally fractured shales and showed that the fracture half-length was a function of the volume of fluid injected in the time interval during which this half-length stopped propagating, after a significant amount of fluid volume is injected into the formation. The length and width of the SRV were observed using micro-seismic monitoring.
Maxwell et al. [39] observed the micro-seismic events monitored during hydraulic fracturing in a naturally fractured shale also of Barnett and found that, from micro-seismic observations that HF occasionally grew at an oblique angle with the direction of fracture assumed. The results also showed that hydraulic fractures grew at an oblique angle because they crossed the natural network of fractures.
Rutqvist and Stephansson [40] have investigated several equivalent continuous models for the fractured medium. One of the simplest approximations is based on a set of three fractures orthogonal to each other in a cube of the porous medium, in which the equations of continuity and elasticity can be solved in an easier way.

UFM models
Roussel and Sharma [41] presented an unconventional twodimensional numerical model and demonstrated that a transverse fracture to a horizontal well can be diverted due to the redistribution of in situ tensions generated by the existence of natural fractures. Thus, a poroelastic model to simulate the stress interference between fractures in horizontal wells was formulated. In these calculations, an HF is represented by the PKN geometry [42,43].
The model was constructed using FLAC-2D (Fast Lagrangian Analysis of Continua) and is based on the explicit finite difference method in 2-D, where each derivative in the set of equations that model the phenomenon is replaced directly by an expression algebra written in terms of field variables e.g. stress or displacement at discrete points in space.
Meyer and Bazan [35] purposed one of the most commonly used unconventional models for fracture network behavior in shale gas/ oil reservoirs: The MShale, which consists of a discrete fracture network (DFN) simulator in naturally fractured formations, e.g. shales and coal bed methane (CBM). The software uses stochastic modeling based on some functions of probability distributions and is used to verify several existing numerical models since it can be used as a diagnostic tool to compare the numerical results of the DFN model with micro seismic data. Thus, through Mshale, the propagation of the network of complex fractures in several planes in shales can be simulated with the specification of the grid.
Belytschko and Black [44] numerically modeled fracture propagation using the Extended Finite Element Method (XFEM), where spacesolution enrichment of the problem with discontinuous functions was used to represent the surface and asymptotic functions for the representation of the other regions of the fracture. XFEM has the advantage of avoiding the problem of remeshing and can find singularities of the field of tensions using the original system of the mesh. The disadvantage of the method is that additional degrees of freedom modify the original structure of the matrix and the need for longer processing time due to the need to integrate numerically each step of the simulation. XFEM is suitable for problems involving the propagation of fractures in heterogeneous media such as shale [45].

Parameters that influence the construction of the network of complex fractures in shale gas/shale oil reservoirs
Net Pressure: Olson and Dahi-Taleghani [46] showed that the net pressure directly influences the construction of the network of complex fractures. Normally the analysis of the influence of the liquid pressure on the interaction between hydraulic fractures and natural fractures is carried out by means of the net pressure coefficient, expressed by:

Dynamic viscosity of the fracturing fluid
The dynamic viscosity of the fracture fluid in shales has a direct influence on the construction of the fracture network, because as the fluid viscosity increases, the extent of the fracture network is significantly reduced [47].
Experimentally, it has been found that a high viscosity fracturing fluid contributes to the generation of the planar fracture pattern and a low viscosity fluid, such as slickwater fracturing, generates the network pattern of complex fractures [47].

Pumping rate
Gomaa et al. [23] investigated the influence of the fracture fluid pumping rate during the slickwater fracturing operation on shales and realized that increasing the flow gives a reduction in the fracturing pressure of the sample.

Distribution and orientation of natural fractures
According to previous considerations, during the hydraulic fracturing process in shale gas/shale oil reservoirs, natural fractures are activated in order to allow the gas trapped in the shale pores to flow to the well and optimize the stimulation operation [48]. Field tests have shown that any HF in a naturally fractured medium is influenced by the approach angles of the NFs and generates a network of complex fractures. Guo et al. [49] verified that the amount and pattern of fracture orientation directly influenced the shape of the network of complex fractures.
Beddings's direction: Beddings are planar sedimentary structures generated by the deposition of lithological layers overlapping each other. In rocks such as shale, which have different types of granulometry and mineralogical composition, the visualization of beddings is more evident [50].
Gomaa et al., [23] have verified experimentally that the direction in which Beddings are accommodated directly influences fracturing pressure and pumping time.
Stress shadowing effect: Ren et al. [8] found experimentally, through a triaxial test in a naturally fractured medium that the HF propagates along the NF to low values of horizontal stress difference and crosses the natural fractures to high values of tension difference horizontal, generating modification in the network behavior of complex fractures. Thus, the presence of a natural fracture close to the hydraulic fracture alters the stress field in the domain considered [51]. The presence of natural fractures also alters the geometry of HF and how it propagates. This effect is called stress shadowing and can change the trajectory of the hydraulic fracture. Thus, as mentioned previously, the spacing between fractures also influences the behavior of the HF being generated.
Interaction among HF-NF: Ren et al. [52] realized that different settings of fracturing influenced at the construction of the complex fractures network, and purposed 2 kinds of fracturing sequences: segmented hydraulic fracturing and multi-cluster hydraulic fracturing.

Mineralogical composition, modulus of elasticity and Poisson coefficient
The shale's brittleness is directly correlated with its mineralogy, modulus of elasticity and Poisson's coefficient. As the concentration of feldspar, quartz and calcite containing silicon or calcium increases, clay concentration decreases and increases the brittleness of the rock, i.e, it reduces its capacity to undergo deformation under the action of tensile stress in the elastic regime [53,54].
Zou et al. [55] realized that for the economically viable exploration of shale gas/ oil and the construction of complex fracture networks, the mineral concentration in the rock should be between 46% and 60%. In order to characterize the brittleness of shales, it is necessary to combine the modulus of elasticity and the Poisson coefficient. Where: is called the "Cauchy stresses's tensor".
According to (Sayers, 2005), for an anisotropic formation like the shale, the modified Hooke's law is given by: Making the change of the coordinate system (x, y, z) → (x 1 , x 2 , x 3 ), the Cauchy's stresses tensor becomes: Realizing the matrix product, we have:  a a a a a a   a a a a a a   a a a a a a   a a a a a a   a a a a a a   a a a a a In the matrix form, one has: The constitutive equations of poro-elasticity for anisotropic formations According to Zhu et al. [56], the constitutive equations for an anisotropic poroelastic model, in terms of effective stresses, are: The rotation matrix [R] θ that transforms the stress coordinates system in situ to the well coordinate system is expressed by: The stresses's state in well coordinates is given by: The tensor of elasticity ψ, in coordinates of the bedding plane, is: For the mechanical stability analysis in sloped wells, the tensile tensor above must be modified for well coordinates.
Thus, the modified elasticity tensor ψ w is given by: Lekhnitskii, Amadei, Aadnøy and Ong [57][58][59][60], have proposed exact solutions to the mechanical problem around sloping wells in shales taking into account the effects of anisotropy. The solutions are given in terms of complex variable functions according to the expressions: The functions of complex variables D', E' and F' are expressed by:

Hydraulic formulation for porous media
Mass conservation: The mass conservation equation, or continuity equation, shows that "the sum of the rate of mass change within the control volume and the convective liquid mass flow on the control surface must be zero." ( ) In a porous media with compressible fluid, the tensor form of the continuity equation is: In the differential form, one has: And finally, in index notation: The mass balance in the pumping of the fracture fluid for the formation is: Nordgren [43] considered the leak-off effect in the Perkins and Kern model, arriving at the following integral-partial differential equation for the mass balance in the injection phase: In terms of the aperture width w(x, t) e of the fracture height h f we have the following partial differential equation: ( Entering variable change: With the global leak-off coefficient C t , expressed by: Or: Finally, the area of space is: The momentum's transport's equation The linear momentum conservation equation in the material form, i.e., following the movement of the fluid within the fracture and in terms of stress, is given by: In terms of speed and convective heat transfer: From the Reynolds's lubrication's theory, for the Hagen-Poiseuille flow, in terms of the aperture, we have the following partial differential equation for the flow of the fracturing fluid inside the fracture: The hydraulic diffusivity equation The equation of mathematically modeling the behavior of pressure as a function of distance and time during the single-phase flow of Newtonian fluids within the pores of a rock is termed the hydraulic diffusivity equation obtained by combining the continuity equation, Darcy's law and equations of fluid and rock matrix compressibility and is expressed in the tensor form, as "the divergence of the pressure gradient": consists of two equations: One equation models the liquid phase and another the solid phase (Figures 3a-3c).

Thermal formulation (The energy equation)
During the flow of the fracture fluid inside the rock pores and later, in the network of complex fractures, convection and conduction   • ∇ (52) And, in heat flow terms by: In tensor notation: Or, in index notation:

Numerical simulation
Petrobras full version of MShale was not available to simulate the whole real subjects about their stimulated wells due information security, thus it could not be used. So, for numerical simulation, the following assumptions for the mechanical, hydraulic, thermal and proppant transport problems were considered: • Linear, isotropic and homogeneous elastic porous medium; Compressibility of the negligible porous medium; • Spurt Loss inconsiderate; • Fracturing fluid with Newtonian, single-phase and incompressible behavior; • Proposed transport fluid with non-Newtonian behavior with power law model; • Leak-off effect considered by the Carter model; • Non-isothermal flow; • Incompressible flow, i.e., ∇ • u= 0; • Flow without active body force; • Azimuth symmetry, i.e., u θ = 0; • Darcian flow in the porous medium; • Non-slip condition on the walls of the fracture, i.e., u¯(0) = 0 • Viscous dissipation Φ / ; • Creeping flow extremely laminar flow during the pumping of the fracture fluid in intact rock and deformed rock; • Hagen-Pouseuille flow inside the fracture modeled by the Reynolds lubrication theory.

DISCUSSION
This work showed that for the stimulation of a given payzone of shale gas/shale oil by means of the construction of a complex network of fractures, several geomechanical, hydraulic and geometric parameters are directly correlated (Table 1). In this way it was verified that: • The higher the net pressure coefficient, the higher the SRV; • In fragile shale gas/shale oil reservoirs, with low values of dynamic viscosity of the fracturing fluid, the construction of the network of complex fractures occurs with greater ease; • As the pumping flow of the fracturing fluid increases, the pressure required to fracture the formation decreases, but in the field, high pressure values are applied to ensure that the formation is fractured. With this, there is a need to hire more trucks equipped with high pressure pumps, increasing the CAPEX of the project; • The operational sequence and intervals between fractures influence in the construction of the network and increase of the SRV, since neighboring fractures disturb the field of tension of the formation, generating the stress effect shadowing, causing a curvature in the fractures located in the ends of the spacing between fractures and reducing the size of the fracture located at its midpoint due to the reorientation of the stress field in the vicinity of the generated fractures and reducing the permeability of this stimulated region. This reorientation of the tensile field also provides a change in the propagation direction of the HFs and negatively influences the construction of the network of complex fractures; • The existence of natural fractures and their orientation angles directly influence the construction of the network of complex fractures, since they also alter the field of formation tensions, so that the larger the number of natural fractures and the more disordered its distribution, the greater the extension of the fracture network. For the approximation angle of 90, the HF crosses the NF as if it did not exist. Thus, by means of microsysmica, one can look for zones with great quantity of natural fractures and with orientations favorable to the construction of the complex network to realize the hydraulic fracturing and to optimize the project of stimulation; • In the fracture settling phase of the fracture, the dynamic viscosity of the fracture fluid must be higher than in the fracture network construction phase, since low dynamic viscosity fracturing fluids, e.g. slickwater, have low transport property proponent; • It was verified that the more fragile the shale, the network of complex fractures is constructed more easily, since the formation undergoes breakout in several directions, connecting more pores and natural fractures; • experimentally verified that parallel beddings provide shale breakout at shorter intervals of time and pressure than perpendicular beddings; • By means of the numerical simulation, for the shale gas models this phenomena is the energy equation expressed in the material form and in temperature terms by: