Geomechanical key parameters of the process of hydraulic fracturing propagation in fractured medium

. Hydraulic Fracturing (HF) is a well-stimulation technique that creates fractures in rock formations through the injection of hydraulically pressurized ﬂ uid. Because of the interaction between HF and Natural Fractures (NFs), this process in fractured reservoirs is different from conventional reservoirs. This paper focuses mainly on three effects including anisotropy in the reservoir, strength parameters of discontinuities, and fracture density on HF propagation process using a numerical simulation of Discrete Element Method (DEM). To achieve this aim, a comprehensive study was performed with considering different situations of in situ stress, the presence of a joint set, and different fracture network density in numerical models. The analysis results showed that these factors play a crucial role in HF propagation process. It also was indicated that HF propagation path is not always along the maximum principal stress direction. The results of the numerical models displayed that the affected area under HF treatment is decreased with increasing the strength parameters of natural fracture and decreasing fracture intensity.


Introduction
Nowadays, fracking fluids (also known as Hydraulic Fracturing [HF]) are utilized extensively in fields including low permeability gas formations, weakly consolidated offshore sediments, "soft" coal beds, and naturally fractured reservoirs to stimulate oil and gas wells (Adachi et al., 2007).
In HF process, fluids are injected into a selected section of the wellbore at a pressure that is high enough to induce a fracturing of the formation. Since the fluid flow generates a pressure differential between the wellbore pressure and the original reservoir pressure, it also creates stress in the formation. By artificially increasing the stress, a point will eventually be reached where the stress becomes greater than the maximum stress that can be sustained by the formation. Consequently, the rock physically splits apart (Economides et al., 2007). By this procedure, existing Natural Fractures (NF) can be reopened and new cracks may be created.
HF treatment is performed in both homogeneous (without fractures) and inhomogeneous (with the presence of NFs) reservoirs. The performance of these two types of reservoirs is completely different during the HF process. In the homogeneous reservoirs, HF propagates to maximum stress direction with a uniform geometry. But, it propagates with different geometries when they encounter NFs under different conditions; the initiation and propagation direction of HF do not always occur along the direction of the maximum horizontal principal in-situ stress (Liu Z. et al., 2014). In reservoirs with NFs, the opening fractures control the fluid flow paths so that the production mechanism in naturally fractured reservoirs is significantly different from that in conventional reservoirs. These NFs may close as the reservoir pressure drops and affect the growth and final geometry of the HF applied for enhancing the production level (Lorenz et al., 1996;Teufel and Clark, 1984).
To simplify the research problem, the PKN model (Nordgren, 1972;Perkins and Kern, 1961) and the KGD model (Geertsma and De Klerk, 1969;Khristianovic and Zheltov, 1955) were developed under the assumption that the reservoir is a homogeneous, isotropic, linear elastic continuum and the fracture geometry is planar. Despite the simplifications in the abovementioned models, several researchers have developed analytical (Blanton, 1986;Renshaw and Pollard, 1995;Sarmadivaleh, 2012;Warpinski and Teufel, 1987), experimental, and numerical models in an effort to examine HF performance with NF. The analytical models were mainly established using differential stresses and approaching angles and were based on simplified assumptions. In experimental laboratory studies (Ai et al., 2018;Dehghan et al., 2015Dehghan et al., , 2016Liu Z. et al., 2014Liu Z. et al., , 2018Wang and Li, 2017;Zhou et al., 2010), blocks with artificial or natural fractures have been tested for investigating the interaction between HF and a single NF. These tests revealed that the conditions required for different interaction modes consist of opening, arresting, and crossing modes.
To simulate the HF process in naturally fractured formations while considering the effect of NFs and multiple fracture interactions in two spaces, numerical models including the Finite Element Method (FEM) (Feng et al., 2015), the extended-FEM (Dahi- Taleghani and Olson, 2011;Zhang et al., 2018a), Cohesive Zone Model (CZM) based on the XFEM (Wang et al., 2018a, b), mixed-FE and the Discontinuous Galerkin Methods (DGM) (Hoteit and Firoozabadi, 2006), Displacement Discontinuous Method (DDM) (Abdollahipour et al., 2015;Behnia et al., 2015;Cheng et al., 2017;Zhang and Jeffrey, 2012), Boundary Element Method (BEM) (Vu et al., 2015), BEM by coupling DDM (Jiang and Cheng, 2018), Displacement Discontinuity Analysis (DDA) (Morgan and Aral, 2015), Distinct Element Method (DEM) by using PFC2D software (Han et al., 2018;Zhou et al., 2017), and hybrid Discrete-FEM (DFEM)  have been proposed. With the development of computing power of computers, the interaction between natural fracture and hydraulic fracture in three spaces was studied. Moinfar et al. (2014) developed an Embedded Discrete Fracture Model (EDFM) for an in-house 3D compositional reservoir simulator that borrows the dual-medium concept from conventional dual-continuum models and also incorporates the effect of each fracture explicitly. Filho et al. (2015) also worked in developing a preprocessing code for the EDFM. Dahi Taleghani et al. (2018) presented an integrated methodology that utilizes the cohesive zone model to simulate the propagation of HF and their interactions with pre-existing NFs. Tang et al. (2018) developed a 3D fracture model based on the three-dimensional DDM to simulate the interactions between vertical/slanted HF and frictional discontinuities (e.g., NFs or horizontal/oblique bedding plane segments) with non-orthogonal approach angle.
According to the simulation mechanism, the mentioned numerical approaches can be classified into two categories: continuum-based methods (FEM, XFEM, CZM, and DDM) and discontinues methods (DEM, DFEM, and EDFM). Compared to continuum-based methods, the basic characteristic of DFEM and DEM is that during each calculation step the model contact patterns are continuously updated. Thus, it is more convenient to mimic fracture initiation, propagation, and intersection by the discontinuous method irrespective to the calculation efficiency (Zhang et al., 2018b).
Based on the results of previous studies, HF processes are highly complex and face many unanswered questions due to a large number of influence factors including the anisotropy of fractured medium, rock structure, and in-situ stress state. Also, considering that many petroleum reservoirs are naturally fractured (Aguilera, 1995), investigating the effects of these NFs on the HF process is of great necessity. Hence, this paper focuses on three effects including a) in-situ stress and discontinuities anisotropy, b) strength parameters of discontinuities, and c) fracture density on HF Propagation (HFP) process. To achieve this goal, many models were simulated using DEM to evaluate the effect of the abovementioned parameters. Finally, the obtained results were discussed.

Discrete element method description
Rock mass is composed of intact materials and discontinuities different sizes that enhance the heterogeneity and anisotropy of rock. The Discrete Element Method (DEM) is a relatively new numerical method in rock mechanics. The essence of this method is to represent the rock as an assemblage of blocks.
In this method, the displacements caused by block motion and rotation are obtained by solving the equations of motion. One significant advantage of this method is that the fracture opening and complete detachment can be explicitly described in the DEM.
For brittle rocks under compressive stress, the fracture initiation and propagation are to a large extent determined by the presence of randomly distributed NFs. The external compressive force is transformed into internal tensile forces between blocks that are perpendicular to the loading direction and shear forces along the boundaries (Fig. 1). Among the two failure modes, the tensile failure mechanism plays a strong role in the fracturing process due to the local in homogeneity of rock .
Blocks and contacts are the fundamental elements for DEM. In the normal direction, the contact forcedisplacement relation is assumed to be linear and the failure condition is governed by the limiting strength and normal stiffness (k n ): where Dr n is the normal stress increment and Du n is the displacement increment including opening and closure in the normal direction. For contacts under tension, if the tensile strength (r t ) is exceeded, the contact is broken and the stress around the crack is redistributed. Similarly, for contact under shearing, the shear stress is proportional to the sliding displacement with respect to the shear stiffness (k s ) : where Ds is the shear stress increment and Áu e s is the elastic shear displacement component. However, the shear strength (s max ) is limited by a combination of cohesion (C) and internal friction angle (u) according to the Mohr-Coulomb criterion.
For analysis of fluid flow through the fractures of a system of impermeable blocks, a fully coupled mechanicalhydraulic analysis is performed, in which fracture conductivity depends on mechanical deformation, and joint fluid pressures affect the mechanical computations. The flow rate is then given by Itasca Consulting Engineers Company (2016): where kj is a joint permeability factor (whose theoretical value is 1/12l); l is the dynamic viscosity of the fluid; a is the contact hydraulic aperture; and l is the length assigned to the contact between the domains.
3 Numerical simulation of HF propagation

Proposed model description
In this study, the HFP is simulated by considering the concept of opening fictitious joints. Accordingly, two types of joint systems were used to simulate HF in fractured media: 1. Fictitious Discrete Fracture Network (FDFN) with equivalent properties using Voronoi element A limiting factor in modeling stress-induced fracturing using UDEC is that all potential fracture pathways must be pre-defined. To incorporate this limitation and provide added degrees of freedom for fracture propagation, a Voronoi tessellation scheme can be used to generate randomly sized polygonal blocks (Zangeneh et al., 2015). Mathematically, the Voronoi tessellation is a collection of entities that fill the space with no overlaps and gaps. The intersection of two polyhedrons is a plane called the "tessellation face" (and in this paper is known as fictitious joints) and that of three polyhedrons is a straight line called the "tessellation edge". It is possible to make the fictitious joints behave as an intact rock, in a global sense by selecting the parameter values for the constitutive models as follows (Kulatilake et al., 1993): (i) The same strength parameter values for both intact rock and the fictitious joints.
(ii) Shear modulus to joint shear stiffness JK S ratio (G/JK S ) between 0.008 and 0.012. (iii) A joint normal stiffness/JK S ratio (JK N /JK S ) between 2 and 3. The most appropriate value in this range may be Young's modulus/G value (E/G) for this particular rock.

NFs with real properties
NFs with real strength properties (C, u, tension, JK N , and JK S ) are added to the model after creating FDFN using Voronoi element. The advantages of the used method in this article are: (i) Simulation of the injection process in a fractured media without defining the roughness parameters of different modes of I and II. (ii) Fractures that are not connected to the boundaries of the model are eliminated using the conventional DEM. For instance, Figure 2 displays the generated network in pre-and post-processing by a common DEM. As it is shown, a large fraction length of the NF network is removed after processing. However, using the presented method in this article, less than 5% of the total length of NFs was eliminated.

Verification of the model
In the first stage of simulations, two HF laboratory experiments carried out by Fatahi et al. (2017) were chosen for verifying the proposed method. These experiments were performed on the cubic synthetic mortar samples of 10 cm with white glue as NF filling material in the true triaxial test cell. HF test was carried out by injecting fluid with a viscosity of 97 700 cp at an injection rate of 0.1 cc/min. Table 1 illustrates the mechanical properties of the synthetic sample and NF filling materials. Figure 3 shows the results of two tests in experimental and numerical methods. In the first test ( Fig. 3a), principal stresses were imposed in a way that HF initiated and propagated in a direction of 30°with respect to NF. No interaction was occurred between HF and NF in both models ( Fig. 3a and b). HF also initiated and propagated in the direction of maximum horizontal principal stress. In the second test, the principal horizontal stresses were rotated 90°with respect to test in Figure 3a. It is observed that both wings of HF crossed NFs and HF propagated in the direction of maximum horizontal stress and intersected the NF at about 60°in the experimental and numerical models.

Surveying of numerical model results
In this paper, the impact of the joint is examined from the following points: (d) The required time to propagate the HF to a specific length or area.
The concept of SRV was initially proposed to provide a quantitative assessment of stimulation effectiveness based on the spatial distribution of microseismic events induced by hydraulic injections. However, such assessment provides little insight into critical parameters such as the hydraulic fracture conductivity (Cipolla and Wallace, 2014), which can vary significantly depending on the rock properties, the local stress regime, and the fracture treatment design. In this paper, SRV is described as a joint opening area. Figure 4 presents the numerical model for the simulation of HF. To achieve this aim, a model with a dimension of 300 Â 300 m 2 was constructed and a region 160 Â 160 m 2 in the center of the model was simulated using FDFN. The injection well was considered at the nearest node in the center of the model by using a FISH program in the UDEC software. The injection rate and fluid viscosity were considered 3 Â 10 À3 m 3 /s and 2 cp, respectively. Table 2 displays considered parameters of intact rock, fictitious joints, and NFs properties.
The main instigator of fracture complexity is the interaction between induced fractures and existing fractures during the HF treatment (Wu and Olson, 2015). Therefore, due to the prevalence of existing NFs in the majority of hydrocarbon reservoirs, anisotropy in these reservoirs can be observed in two cases: (i) Anisotropy due to the presence of the NFs or joint set. (ii) Anisotropy in in-situ stresses (maximum to minimum horizontal stress ratio K = r H /r h ).
For this purpose, some models were simulated with different K of 1, 1.5, and 2 and different approaching angles (h) of 0, 15, 30, 45, 60, 75, and 90°. A model without the presence of the joint set was also considered to examine its impact on the process of HF propagation. For example, Figure 5 shows the process of HFP (open fictitious joints) in hydrostatic condition with h = 45°during the injection as well as displacement vectors. According to this figure, the length and effective area increase with increasing injection time (volume of injected fluid). Also, the HF is propagated in the joint set direction. In other words, the direction of the joints determines the path of expansion.

The effect of anisotropy due to the presence of joint set on the path of HFP
By examining the results of numerical modeling, the HF pathway can be propagated in one of two directions: 1. In the direction of the maximum principal stress (r H ). 2. In the direction of the existing joint set in the model/reservoir. Figure 6 shows the position of the HFP at the end of the injection time at K = 1.5 for different approaching angles. According to this figure, the direction of joints set determines the HF path at angles of less than 30°while the path of HF is parallel to maximum stress at an angle of more than 30°.
The results are briefly presented in Figure 7 for different stress coefficient K h = (r H Àr h )/r h under different   Figure 7, it is understood that at low angles with different K h ratios and at low K h with different angles, the path of HFP is in the direction of the joints set and it is distorted due to the maximum principal stress. However, at h > 30°and K h > 0.5 (the hatched part of the graph), the HF is propagated in the direction of the r H and joint sets do not affect the direction of HFP. It has to be noted that in some models, it is possible to stop/arrest the HF by NFs. This mode has been discussed via the interaction of single NF with HF in the previous investigations (Renshaw and Pollard 1995;Sarmadivaleh andRasouli, 2014, 2015). This event also occurs in the models with the presence of a joint set. For instance, in the model of K = 1.5 and h = 45°, the arresting mode for HF has been created. Figure 6d shows that the HFP decreased in the same injection time compared to other models; therefore, it stopped. This phenomenon can be explained in the form of shear displacements (Fig. 8). According to Figure 8, on the upper and lower joints of the injection borehole, long shear displacements are observed that indicate a slip on these joints. In other words, a large amount of injection fluid is used to slip the joint instead of crossing and opening, and the extended length of the HF (i.e., SRV) decreases in this situation.

The anisotropic effect in the shape of HFP
Usually, the area under the influence of the fluid injection process in the geological environment is oval (with two diameters of a and b). The region affected in this study is a zone where joints (both FDFN and NFs/joints set) are opened due to the fluid injection process or the pore pressure in the joints are changed due to the initial pressure condition of media.
The numerical modeling performed in this study shows that the stress ratio (K) has a significant effect on the shape of the region. Figure 9 presents the opening of fractures in two conditions of h = 0°and 90°with three different stress ratios of 1, 1.5, and 2. According to this figure, the more the stress ratio is closer to the isotropic (hydrostatic) state, the more circular-from (i.e. a % b) the affected area (SRV). Fig. 3. a) and b) Experimental and numerical results of two samples with white glue as NF filling material for NF at 30°with respect to HF, respectively; c) and d) Experimental and numerical result for NF at 60°with respect to HF, respectively.
When K is far from the hydrostatic state (K > 1), the difference in size of two oval diameters is more, and ultimately it looks like a line. Figure 10 provides a schematic representation of the general states of the affected area under the influence of the HF process. The following results can be generally obtained from this figure: 1. K = 1 (a hydrostatic state): Joints set determines the path of HFP. In this case, the larger diameter is parallel to the direction of the joints set. Also, the difference between the two diameters is small and similar to the circle. 2. K = 1.5: At h < 45°, the direction of HFP is dominated by the direction of the joints set, and the larger diameter of the ellipse is aligned with the joints direction.
However, at h > 45°, the path of HFP is parallel to the maximum principal stress. The difference between the oval diameters decreases and approaches linearly when the joint angle is closer to the r H . 3. K = 2: The path of HFP is in the direction of the r H and the larger diameter of the ellipse is in this direction. The difference in oval diameters is greater than other states (K = 1 and 1.5).

The effect of anisotropy on the SRV
To study the effect of anisotropy, two general conditions were considered: (i) The amount of affected area (SRV) in a constant time. (ii) The required time to propagate the HF in a specific SRV. Figure 11 presents the effect of a joint set dip in SRV, as well as the maximum opening of NFs at times of 1200, 1600, 2000 s in the stress condition of K = 1.5. The SRV has the highest amount at approaching angles of 15-30°and, consequently, the maximum opening at the end of the injection time occurs at this angle. With increasing h, the SRV value decreases sharply at 45°. The reason is that HFP is arrested in the interaction with NF in this situation. Also, when a joint set direction is parallel or perpendicular to the r H , the SRV and the maximum opening is approximately equal to each other. This consequence is similar to results of Jaeger and Cook (1979) Figure 13 shows two simultaneous effects of anisotropy on the stress ratio and the angles between joints set and r H at a specific SRV. This diagram represents the time it takes the SRV of the model/reservoir to reach a specific value in different conditions. The following points can be taken from this chart: (a) In general, with increasing the stress ratio, the time to reach the specific SRV increases. In other words, in the hydrostatic condition, the fluid more easily flows into the fractures system at the same h angle. Therefore, the SRV becomes more at a specific time. This point is derived from the models without joints set (dotted lines in Fig. 13). (b) At angle h = 30°, the stress anisotropy effect is eliminated and the time to achieve the specified SRV is identical.   Jaeger and Cook (1979) diagram, the sample has the least resistance at h = 30°in the uniaxial condition (r h = 0 ? K % 1). Here, at K = 2, the model exhibits the least resistance such that, in this case, the time to reach the specified SRV is minimal. When the stress condition reaches the hydrostatic state, the least strength occurs at about 60°. (d) The arresting mode is created in special conditions because of the interaction of HF with NFs. In this situation, HF cannot easily cross through NFs. Therefore, considerable time (more fluid volume) is used to increase fluid pressure so that SRV increases. In simulated models, this phenomenon happened for conditions of K = 1.5 and 2 under approaching angles of 45-60°. Hence, the time-h diagram is vertical asymptotically formed. (e) Comparing the results of the models without fractures with fractured models reveals that the impact of joints set in the HF process can be either effective or ineffective. It is effective when the SRV increases at a specific time or the injected volume of fluid; i.e., reducing the time to reach a specified length or SRV. Otherwise, it is ineffective.
In a hydrostatic model, the required time to reach the SRV of 2000 m 2 is less in the conditions of 60 < h < 70°a nd without joints set. In other situations, the presence of NFs reduces SRV or increases the time or injected volume of fluid. Thus, it can be claimed that the joints are ineffective in HF process in hydrostatic conditions.
With increasing K, the h range that acts effectively in the HF process is increased. In K = 1.5, this range is within 20 < h < 35 while in K = 2 this range is 5 < h < 35 and h > 78°(about half of the models of K = 2). So, in the combination of high-stress anisotropy (K = 2) and anisotropy due to the presence of the joints set, the interaction of the fluid and the fractures network system in the reservoir has an effective role in the HF process. So, due to the difficulty involved in understanding the interaction of HF with NFs, it is necessary to simulate each reservoir based on factors such as in-situ stress, joints set, and rock properties of the reservoir and fluid properties.

The effect of strength parameters of natural joints set on the HFP process
In this section, the effect of joint strength parameters including cohesion (c), friction angle (u), and tensile strength (T) on the HFP process is investigated. To achieve this goal, the hydrostatic model with approaching angle of 30°was used for the sensitivity analysis. Figure 14 demonstrates the effect of these parameters on the SRV during the injection time of 1600 s. According to this figure, with the increase in each strength parameter, the SRV level decreases, because more force (fluid pressure) is needed to open the NFs. Therefore, the number of opened contact elements is reduced at a specified time. Also, the effect of friction angle decreases when u exceeds 30°. For example, Figure 15 indicates SRV in a specific time of 1600 s for the friction angle of 20 and 35°. In addition to the higher SRV in less strength (u = 20°), it is clear that the fluid due to the need for less energy at the contact surface (NFs) tends to open the NFs; and HF easier propagates along the NF. When the friction angle is increased (i.e., the fractures are welded together), more force is needed for opening the fractures. In this condition, the fracture-opening does not occur in the direction of NFs and, eventually, SRV decreases. Therefore, the fluid injection time should be increased to reach the specified SRV.

The effect of fracture density on the process of HFP
The density or intensity of fractures network is one of the most effective geometric parameters that can affect rock mass behavior under external loads, such as injection pressure. This parameter can be defined in one, two, or three dimensions. In two dimensions, density is defined as the number of fractures per unit area (P20) and intensity as the sum of the lengths of discontinuities per unit area (P21) (Dershowitz and Herda, 1992). In this study, the power law (in the form of Eq. (5) is used for fracture distribution based on fractal phenomena: where D is fractal dimension, L min and L max are the smallest and largest length fractures, respectively, F denotes the random probability of a uniform distribution in the range 0 < F < 1; which correlated based on field mapping in Sellafield (Min et al., 2004  Several previous studies have proven the use of fractal law for a variety of geological phenomena, including the distribution of discontinuities, faults, and earthquakes (Turcotte, 1997).
In order to model the fracture distribution, a code was written in the MATLAB program based on the power distribution law and its output was used as input in the UDEC software. Then, it is possible to simulate the HFP using the concept of FDFN and combining it with the NF network generated in the program. Four models were constructed with different densities/ intensities using the described method (Fig. 16). Then, each model was combined with the reference FDFN and, eventually, the fluid injection process was simulated in the reservoir. In these simulations, the dimensions of the models increased to 800 Â 800 m 2 and FDFN and NFs were embedded in a square of 300 Â 300 m 2 in the center of the model. Figures 17 and 18 present the SRV-time graph and the effect of fracture density/intensity on the HF process at   500 s for different models, respectively. As shown in Figures  16 and 17, increasing the fracture density results in an increase in the SRV at a specific time. In other words, to stimulate the reservoir at a specific level or volume, more time is needed in a reservoir with lower density. This can be related to the effective elastic modulus of the fractured rock mass and the hydraulic conductivity of the fractures network. Walsh (1965) correlates the effective elastic modulus for fractured media with the following equation:  From the above relation, it is assumed that by increasing the number of fractures per unit area (i.e., increasing the fracture intensity), the effective Young's modulus of the fractured media decreases. Another study has also shown that with increasing the fracture density, the shear modulus decreases (Eshelby, 1957;Xu et al., 2018). Since the shear modulus has a direct relationship with Young's modulus, by decreasing Young's modulus, the shear modulus also decreases. In additions, an increase in the density of fractures leads to the corresponding increase in the hydraulic conductivity of the reservoir system (Molebatsi et al., 2009;Namdari et al., 2016). Therefore, by simultaneous incorporation of the effects of Young's modulus and hydraulic conductively and with increasing the fracture density, HF can be propagated more easily and SRV will eventually increase with the same injection rate and stress ratio.

Conclusion
Natural Fractures (NFs) have a significant influence on the HFP in fractured reservoirs. Also, understanding the interaction of HF with NFs is a complicated matter that cannot be performed by analyzing models alone and requires numerical modeling for each reservoir.
This paper focuses on the effects of geomechanical key parameters on the HFP process using the DEM numerical simulation. To achieve this goal, the fictitious joints of Voronoi elements were considered in the models with a combination of NFs including joints set and fracture network. The main results of the presented study is related to anisotropy in reservoirs and the effect of fracture density are as follows: The growth direction in the naturally fractured reservoir is a complex phenomenon that is difficult to predict with high certainty. The HFP path is not always in the direction of maximum principal stress. The joints set direction determines the HFP path in the conditions of low in-situ stress ratio (K = 1) and moderate stress ratio (K = 1.5) with approaching angle less than 45°. Generally, the SRV is decreased with increasing the strength parameters of NFs including cohesion, friction angle, and tensile strength. The results of numerical models showed that the SRV is increased with increasing the NF density.