ogst
Home All issues Volume 74 (2019) Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles, 74 (2019) 58 Full HTML
Open Access
Issue
Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles
Volume 74, 2019
Article Number 58
Number of page(s) 16
DOI https://doi.org/10.2516/ogst/2019026
Published online 21 June 2019

© R. Basirat et al., published by IFP Energies nouvelles, 2019

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 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., 2015, 2016; Liu Z. et al., 2014, 2018; Wang 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) (Liu C. et al., 2018) 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.

2 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 (Li et al., 2017).

thumbnail Fig. 1

Failure modes in tessellation media: a) tensile and shear stresses induced by compression force, b) the idealized relation between hydraulic aperture and normal stress for a rock joint, c) tensile rupture mode, and d) shear rupture mode.

Blocks and contacts are the fundamental elements for DEM. In the normal direction, the contact force-displacement relation is assumed to be linear and the failure condition is governed by the limiting strength and normal stiffness (k n):

(1)

where Δσ n is the normal stress increment and Δu n is the displacement increment including opening and closure in the normal direction.

For contacts under tension, if the tensile strength (σ 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 ) (Li et al., 2017):

(2)

(3)

where Δτ is the shear stress increment and is the elastic shear displacement component. However, the shear strength (τ max) is limited by a combination of cohesion (C) and internal friction angle (φ) according to the Mohr-Coulomb criterion.

For analysis of fluid flow through the fractures of a system of impermeable blocks, a fully coupled mechanical-hydraulic 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):

(4)

where kj is a joint permeability factor (whose theoretical value is 1/12μ); μ 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

3.1 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):

    1. The same strength parameter values for both intact rock and the fictitious joints.

    2. Shear modulus to joint shear stiffness JKS ratio (G/JKS) between 0.008 and 0.012.

    3. A joint normal stiffness/JKS ratio (JKN/JKS) between 2 and 3. The most appropriate value in this range may be Young’s modulus/G value (E/G) for this particular rock.

  2. NFs with real properties

    NFs with real strength properties (C, φ, tension, JKN, and JKS) are added to the model after creating FDFN using Voronoi element. The advantages of the used method in this article are:

    1. Simulation of the injection process in a fractured media without defining the roughness parameters of different modes of I and II.

    2. 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.

thumbnail Fig. 2

Generated fracture network before and after processing by common DEM (La Pointe, 1988).

3.2 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.

thumbnail 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.

Table 1

Mechanical properties of synthetic sample and natural fracture filling materials (Fatahi et al., 2017).

3.3 Surveying of numerical model results

In this paper, the impact of the joint is examined from the following points:

  1. The path of HFP.

  2. Stimulated Reservoir Volume (SRV) or affected area by HF treatment.

  3. The shape of the affected area in the process of HFP.

  4. 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 m2 was constructed and a region 160 × 160 m2 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 m3/s and 2 cp, respectively. Table 2 displays considered parameters of intact rock, fictitious joints, and NFs properties.

thumbnail Fig. 4

Numerical model dimensions for investigating the effect of joint sets on HF.

Table 2

Rock mechanical properties of rock mass, FDFN, and NFs.

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:

  1. Anisotropy due to the presence of the NFs or joint set.

  2. Anisotropy in in-situ stresses (maximum to minimum horizontal stress ratio K = σ H/σ h).

For this purpose, some models were simulated with different K of 1, 1.5, and 2 and different approaching angles (θ) 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 θ = 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.

thumbnail Fig. 5

The hydraulic opening of fictitious joints and joint set in times of 400, 800, 1200, 1600, 2000, 2400, and 2500 s in hydrostatic condition and θ = 45°.

3.3.1 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 (σ 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°.

thumbnail Fig. 6

The HFP path in the model of K = 1.5 and different approaching angles at the injection time of approximately 2400 s.

The results are briefly presented in Figure 7 for different stress coefficient K h = (σ H−σ h)/σ h under different approaching angles. Based on 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 θ > 30° and K h > 0.5 (the hatched part of the graph), the HF is propagated in the direction of the σ H and joint sets do not affect the direction of HFP.

thumbnail Fig. 7

Interaction mode between NF and HF under different conditions.

It has to be noted that in some , 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 and Rasouli, 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 θ = 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.

thumbnail Fig. 8

Shear displacements in the model of K = 1.5 and θ = 45°.

3.3.2 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 θ = 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). 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:

thumbnail Fig. 9

The opening of fractures in two conditions of θ = 0° and 90° with three different stress ratios of 1, 1.5, and 2.

thumbnail Fig. 10

General conditions of the affected area under the influence of the HF process in a schematic way.

  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 θ < 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 θ > 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 σ H.

  3. K = 2: The path of HFP is in the direction of the σ 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).

3.3.3 The effect of anisotropy on the SRV

To study the effect of anisotropy, two general conditions were considered:

  1. The amount of affected area (SRV) in a constant time;

  2. 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 θ, 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 σ H, the SRV and the maximum opening is approximately equal to each other. This consequence is similar to results of Jaeger and Cook (1979) in relation to the anisotropy effect of bedding on the uniaxial compressive strength and research results of Mousavi Nezhad et al. (2018) and Mighani et al. (2016) in relation to the anisotropy effect of bedding on the tensile strength of the rock samples (Fig. 12). Their results depicted that when the angle between bedding and maximum principal stress is 30°, the sample strength has the lowest value. In additions, when the bedding/joint set direction is 15–30°, the rock mass strength decreases in the condition of from the intact rock or θ > 30°. Therefore, less injection force (less fluid pressure) is required for the opening of the fractures.

thumbnail Fig. 11

The effect of a) joint set dip in SRV and b) maximum opening of NFs at different times in the model of K = 1.5.

thumbnail Fig. 12

The effect of discontinuity dip on uniaxial compressive strength (Jaeger and Cook, 1979) and tensile strength for shale rock (Mighani et al., 2016) and for Lyons sandstone and Pyrophylitic (Mousavi Nezhad et al., 2018).

Figure 13 shows two simultaneous effects of anisotropy on the stress ratio and the angles between joints set and σ 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:

thumbnail Fig. 13

The effect of joints set direction on the required time for injection in different conditions.

  1. 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 θ 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).

  2. At angle θ = 30°, the stress anisotropy effect is eliminated and the time to achieve the specified SRV is identical.

  3. The concave point of the graph is shifted from 60 to 30° by increasing the stress ratio. Referring to the Jaeger and Cook (1979) diagram, the sample has the least resistance at θ = 30° in the uniaxial condition (σ h = 0 → K ≈ ∞). 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°.

  4. 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-θ diagram is vertical asymptotically formed.

  5. 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 m2 is less in the conditions of 60 < θ < 70° and 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 θ range that acts effectively in the HF process is increased. In K = 1.5, this range is within 20 < θ < 35 while in K = 2 this range is 5 < θ < 35 and θ > 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.

3.3.4 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 (φ), 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 φ 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 (φ = 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.

thumbnail Fig. 14

The effect of strength parameters on SRV.

thumbnail Fig. 15

HFP at 1600 s in a model of a) φ = 35° (left) and b) φ = 20° (right).

4 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:

(5)

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 m2 and FDFN and NFs were embedded in a square of 300 × 300 m2 in the center of the model.

thumbnail Fig. 16

Four simulated models with different densities/intensities.

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.

thumbnail Fig. 17

SRV graph with time for four models.

thumbnail Fig. 18

The effect of fracture density/intensity on the HF process at 500 s for different models.

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:

(6)

where E′ is the effective Young’s modulus, E is Young’s modulus of intact rock, v is the Poisson ratio, and N is the number of fractures in the studied area of A. Also, the phrase denotes the fracture intensity per unit area (Li et al., 2013).

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.

5 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.

References

  • Abdollahipour A., Fatehi Marji M., Yarahmadi Bafghi A., Gholamnejad J. (2015) Simulating the propagation of hydraulic fractures from a circular wellbore using the displacement discontinuity method, Int. J. Rock Mech. Min. Sci. 80, 281–291. doi: 10.1016/j.ijrmms.2015.10.004. [CrossRef] [Google Scholar]
  • Adachi J., Siebrits E., Peirce A., Desroches J. (2007) Computer simulation of hydraulic fractures, Int. J. Rock Mech. Min. Sci. 44, 5, 739–757. doi: 10.1016/j.ijrmms.2006.11.006. [CrossRef] [Google Scholar]
  • Aguilera R. (1995) Naturally fractured reservoirs, PennWell Publishing Company, Tulsa. [Google Scholar]
  • Ai C., Li X.-X., Zhang J., Jia D., Tan W.-J. (2018) Experimental investigation of propagation mechanisms and fracture morphology for coalbed methane reservoirs, Pet. Sci. 15, 4, 815–829. doi: 10.1007/s12182-018-0252-z. [CrossRef] [Google Scholar]
  • Behnia M., Goshtasbi K., Fatehi Marji M., Golshani A. (2015) Numerical simulation of interaction between hydraulic and natural fractures in discontinuous media, Acta Geotech. 10, 4, 533–546. doi: 10.1007/s11440-014-0332-1. [CrossRef] [Google Scholar]
  • Blanton T.L. (1986) Propagation of hydraulically and dynamically induced fractures in naturally fractured reservoirs, SPE Unconventional Gas Technology Symposium, Society of Petroleum Engineers, Louisville, Kentucky. doi: 10.2118/15261-MS. [Google Scholar]
  • Cheng W., Jiang G., Jin Y. (2017) Numerical simulation of fracture path and nonlinear closure for simultaneous and sequential fracturing in a horizontal well, Comput. Geotech. 88, 2017, 242–255. doi: 10.1016/j.compgeo.2017.03.019. [CrossRef] [Google Scholar]
  • Cipolla C., Wallace J. (2014) Stimulated reservoir volume: A misapplied concept? SPE Hydraulic Fracturing Technology Conference, Society of Petroleum Engineers, The Woodlands, Texas, USA. doi: 10.2118/168596-MS. [Google Scholar]
  • Dahi-Taleghani A., Olson J.E. (2011) Numerical modeling of multistranded-hydraulic-fracture propagation: Accounting for the interaction between induced and natural fractures, Soc. Pet. Eng. 16, 3, 575–581. doi: 10.2118/124884-PA. [Google Scholar]
  • Dahi Taleghani A., Gonzalez-Chavez M., Yu H., Asala H. (2018) Numerical simulation of hydraulic fracture propagation in naturally fractured formations using the cohesive zone model, J. Pet. Sci. Eng. 165, 42–57. doi: 10.1016/j.petrol.2018.01.063. [CrossRef] [Google Scholar]
  • Dehghan A.N., Goshtasbi K., Ahangari K., Jin Y. (2015) Experimental investigation of hydraulic fracture propagation in fractured blocks, Bull. Eng. Geol. Environ. 74, 3, 887–895. doi: 10.1007/s10064-014-0665-x. [CrossRef] [Google Scholar]
  • Dehghan A.N., Goshtasbi K., Ahangari K., Jin Y. (2016) Mechanism of fracture initiation and propagation using a tri-axial hydraulic fracturing test system in naturally fractured reservoirs, Eur. J. Environ. Civ. Eng. 20, 5, 560–585. doi: 10.1080/19648189.2015.1056384. [CrossRef] [Google Scholar]
  • Dershowitz W.S., Herda H.H. (1992) Interpretation of fracture spacing and intensity. 33rd U.S. Symposium on Rock Mechanics (USRMS), 3–5 June, Santa Fe, New Mexico. [Google Scholar]
  • Economides M.J., Mikhailov D.N., Nikolaevskiy V.N. (2007) On the problem of fluid leakoff during hydraulic fracturing, Transp. Porous Media 67, 3, 487–499. doi: 10.1007/s11242-006-9038-7. [CrossRef] [Google Scholar]
  • Eshelby J.D. (1957) The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proc. Royal Soc. Lond. Ser. A. Math. Phys. Sci. 241, 1226, 376–396. doi: 10.1098/rspa.1957.0133. [MathSciNet] [Google Scholar]
  • Fatahi H., Hossain M., Sarmadivaleh M. (2017) Numerical and experimental investigation of the interaction of natural and propagated hydraulic fracture, J. Nat. Gas Sci. Eng. 37, 2017, 409–424. doi: 10.1016/j.jngse.2016.11.054. [CrossRef] [Google Scholar]
  • Feng Y., Arlanoglu C., Podnos E., Becker E., Gray K.E. (2015) Finite-element studies of hoop-stress enhancement for wellbore strengthening, Soc. Pet. Eng. 30, 1, 1–14. doi: 10.2118/168001-PA. [Google Scholar]
  • Filho J.S., Shakiba M., Moinfar A., Sepehrnoori K. (2015) Implementation of a preprocessor for embedded discrete fracture modeling in an impec compositional reservoir simulator, SPE Reservoir Simulation Symposium, 23–25 February, Houston, Texas, USA. doi: 10.2118/173289-MS. [Google Scholar]
  • Geertsma J., De Klerk F. (1969) A rapid method of predicting width and extent of hydraulically induced fractures, J. Pet. Technol. 21, 12, 1571–1581. [CrossRef] [Google Scholar]
  • Han Z., Zhou J., Zhang L. (2018) Influence of grain size heterogeneity and in-situ stress on the hydraulic fracturing process by PFC2D modeling, Energies 11, 6, 1413. doi: 10.3390/en11061413. [CrossRef] [Google Scholar]
  • Hoteit H., Firoozabadi A. (2006) Compositional modeling of discrete-fractured media without transfer functions by the discontinuous Galerkin and mixed methods, SPE J. 11, 3, 341–352. doi: 10.2118/90277-PA. [CrossRef] [Google Scholar]
  • Itasca Consulting Engineers Company (2016) UDEC 6 Manual. [Google Scholar]
  • Jaeger J.C., Cook N.G.W. (1979) Fundamental of rock mechanics, Chapman and Hall, London. [Google Scholar]
  • Jiang G., Cheng W. (2018) Hydraulic fracture deflection at bedding plane due to the non-orthogonal propagation and the dissimilar material properties, Arab. J. Sci. Eng. 43, 6535. doi: 10.1007/s13369-018-3291-2. [CrossRef] [Google Scholar]
  • Khristianovic S., Zheltov Y. (1955) Formation of vertical fractures by means of highly viscous fluids, in: Proceedings of the 4th World Petroleum Congress, 6–15 June, Rome, Italy, Vol. 2, pp. 579–586. [Google Scholar]
  • Kulatilake P.H.S.W., Wang S., Stephansson O. (1993) Effect of finite size joints on the deformability of jointed rock in three dimensions, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 30, 5, 479–501. doi: 10.1016/0148-9062(93)92216-D. [CrossRef] [Google Scholar]
  • La Pointe P.R. (1988) A method to characterize fracture density and connectivity through fractal geometry, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 25, 6, 421–429. doi: 10.1016/0148-9062(88)90982-5. [CrossRef] [Google Scholar]
  • Li X.F., Li H.B., Zhao J. (2017) 3D polycrystalline discrete element method (3PDEM) for simulation of crack initiation and propagation in granular rock, Comput. Geotech. 90, 96–112. doi: 10.1016/j.compgeo.2017.05.023. [CrossRef] [Google Scholar]
  • Li Y., Wei C., Qin G., Li M., Luo K. (2013) Numerical simulation of hydraulically induced fracture network propagation in shale formation, International Petroleum Technology Conference, 26–28 March, Beijing, China. doi: 10.2523/IPTC-16981-MS. [Google Scholar]
  • Liu C., Jin X., Shi F., Lu D., Liu H., Wu H. (2018) Numerical investigation on the critical factors in successfully creating fracture network in heterogeneous shale reservoirs, J. Nat. Gas Sci. Eng. 59, 427–439. doi: 10.1016/j.jngse.2018.09.019. [CrossRef] [Google Scholar]
  • Liu Z., Chen M., Zhang G. (2014) Analysis of the influence of a natural fracture network on hydraulic fracture propagation in carbonate formations, Rock Mech. Rock Eng. 47, 2, 575–587. doi: 10.1007/s00603-013-0414-7. [CrossRef] [Google Scholar]
  • Liu Z., Wang S., Zhao H., Wang L., Li W., Geng Y., Chen M. (2018) Effect of random natural fractures on hydraulic fracture propagation geometry in fractured carbonate rocks, Rock Mech. Rock Eng. 51, 2, 491–511. doi: 10.1007/s00603-017-1331-y. [CrossRef] [Google Scholar]
  • Lorenz J.C., Warpinski N.R., Teufel L.W. (1996) Natural fracture characteristics and effects, Lead. Edge 15, 8, 909–911. doi: 10.1190/1.1437388. [CrossRef] [Google Scholar]
  • Mighani S., Sondergeld C.H., Rai C.S. (2016) Observations of tensile fracturing of anisotropic rocks, SPE J. 21, 4, 1289–1301. doi: 10.2118/2014-1934272-PA. [CrossRef] [Google Scholar]
  • Min K.B., Jing L., Stephansson O. (2004) Determining the equivalent permeability tensor for fractured rock masses using a stochastic REV approach: Method and application to the field data from Sellafield, UK, Hydrogeology 12, 497–510. doi: 10.1007/s10040-004-0331-7. [CrossRef] [Google Scholar]
  • Moinfar A., Varavei A., Sepehrnoori K., Johns R.T. (2014) Development of an efficient embedded discrete fracture model for 3D compositional reservoir simulation in fractured reservoirs, SPE J. 19, 2, 289–303. doi: 10.2118/154246-PA. [CrossRef] [Google Scholar]
  • Molebatsi T., Galindo Torres S., Li L., Bringemeier D., Wang X. (2009) Effect of fracture permeability on connectivity of fracture networks, International Mine Water Conference, 19th–23rd October, Pretoria, South Africa. [Google Scholar]
  • Morgan W.E., Aral M.M. (2015) Modeling hydraulic fracturing in naturally fractured reservoirs using the discontinuous deformation analysis, 49th US Rock Mechanics/Geomechanics Symposium, 28 June–1 July, San Francisco, USA. [Google Scholar]
  • Mousavi Nezhad M., Fisher Q.J., Gironacci E., Rezania M. (2018) Experimental study and numerical modeling of fracture propagation in shale rocks during Brazilian disk test, Rock Mech. Rock Eng. 51, 6, 1755–1775. doi: 10.1007/s00603-018-1429-x. [CrossRef] [Google Scholar]
  • Namdari S., Baghbanan A., Habibi M.J. (2016) Effects of matrix permeability and fracture density on flow pattern in dual porous rock masses, EUROCK, 29–31 August, Cappadocia, Turkey. [Google Scholar]
  • Nordgren R.P. (1972) Propagation of a vertical hydraulic fracture, Soc. Pet. Eng. J. 12, 4, 306–314. [CrossRef] [Google Scholar]
  • Perkins T.K., Kern L.R. (1961) Widths of hydraulic fractures, J. Pet. Technol. 13, 9, 937–949. [CrossRef] [Google Scholar]
  • Renshaw C.E., Pollard D.D. (1995) An experimentally verified criterion for propagation across unbounded frictional interfaces in brittle, linear elastic materials, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 32, 3, 237–249. doi: 10.1016/0148-9062(94)00037-4. [CrossRef] [Google Scholar]
  • Sarmadivaleh M. (2012) Experimental and numerical study of interaction of a pre-existing natural interface and an induced hydraulic fracture, PhD Thesis, Curtin University, Australia. [Google Scholar]
  • Sarmadivaleh M., Rasouli V. (2014) Modified Renshaw and Pollard criteria for a non-orthogonal cohesive natural interface intersected by an induced fracture, Rock Mech. Rock Eng. 47, 6, 2107–2115. doi: 10.1007/s00603-013-0509-1. [CrossRef] [Google Scholar]
  • Sarmadivaleh M., Rasouli V. (2015) Test design and sample preparation procedure for experimental investigation of hydraulic fracturing interaction modes, Rock Mech. Rock Eng. 48, 1, 93–105. doi: 10.1007/s00603-013-0543-z. [CrossRef] [Google Scholar]
  • Tang J., Wu K., Li Y., Hu X., Liu Q., Ehlig-Economides C. (2018) Numerical investigation of the interactions between hydraulic fracture and bedding planes with non-orthogonal approach angle, Eng. Fract. Mech. 200, 1–16. doi: 10.1016/j.engfracmech.2018.07.010. [CrossRef] [Google Scholar]
  • Teufel L.W., Clark J.A. (1984) Hydraulic fracture propagation in layered rock: Experimental studies of fracture containment, Soc. Pet. Eng. J. 24, 1, 19–32. doi: 10.2118/9878-PA. [CrossRef] [Google Scholar]
  • Turcotte L.D. (1997) Fractals and chaos in geology and geophysics, Vol. 131, Cambridge University Press, Cambridge. [CrossRef] [Google Scholar]
  • Vu M.N., Nguyen S.T., Vu M.H. (2015) Modeling of fluid flow through fractured porous media by a single boundary integral equation, Eng. Anal. Bound. Elem. 59, 166–171. [CrossRef] [Google Scholar]
  • Walsh J.B. (1965) The effect of cracks on the uniaxial elastic compression of rocks, J. Geophys. Res. 70, 2, 399–411. doi: 10.1029/JZ070i002p00399. [CrossRef] [Google Scholar]
  • Wang B., Zhou F., Wang D., Liang T., Yuan L., Hu J. (2018b) Numerical simulation on near-wellbore temporary plugging and diverting during refracturing using XFEM-Based CZM, J. Nat. Gas Sci. Eng. 55, 368–381. doi: 10.1016/j.jngse.2018.05.009. [CrossRef] [Google Scholar]
  • Wang X.L., Shi F., Liu C., Lu D.T., Liu H., Wu H.A. (2018a) Extended finite element simulation of fracture network propagation in formation containing frictional and cemented natural fractures, J. Nat. Gas Sci. Eng. 50, 309–324. doi: 10.1016/j.jngse.2017.12.013. [CrossRef] [Google Scholar]
  • Wang Y., Li C.H. (2017) Investigation of the effect of cemented fractures on fracturing network propagation in model block with discrete orthogonal fractures, Rock Mech. Rock Eng. 50, 1851–1862. doi: 10.1007/s00603-017-1198-y. [CrossRef] [Google Scholar]
  • Warpinski N.R., Teufel L.W. (1987) Influence of geologic discontinuities on hydraulic fracture propagation (includes associated papers 17011 and 17074), J. Pet. Technol. 39, 2, 209–220. doi: 10.2118/13224-PA. [CrossRef] [Google Scholar]
  • Wu K., Olson J.E. (2015) Numerical investigation of complex hydraulic fracture development in naturally fractured reservoirs, SPE Hydraulic Fracturing Technology Conference, 3–5 February, The Woodlands, Texas, USA. doi: 10.2118/173326-MS. [Google Scholar]
  • Xu S., Tang X., Torres-Verdín C., Su Y. (2018) Seismic shear wave anisotropy in cracked rocks and an application to hydraulic fracturing, Geophys. Res. Lett. 45, 11, 5390–5397. doi: 10.1029/2018GL077931. [CrossRef] [Google Scholar]
  • Zangeneh N., Eberhardt E., Bustin R.M. (2015) Investigation of the influence of natural fractures and in situ stress on hydraulic fracture propagation using a distinct-element approach, Can. Geotech. J. 52, 926–946. doi: 10.1139/cgj-2013-0366. [CrossRef] [Google Scholar]
  • Zhang B., Ji B., Liu W. (2018a) The study on mechanics of hydraulic fracture propagation direction in shale and numerical simulation, Geomech. Geophys. Geo-energy Geo-Resour. 4, 2, 119–127. doi: 10.1007/s40948-017-0077-z. [CrossRef] [Google Scholar]
  • Zhang L., Zhou J., Braun A., Han Z. (2018b) Sensitivity analysis on the interaction between hydraulic and natural fractures based on an explicitly coupled hydro-geomechanical model in PFC2D, J. Pet. Sci. Eng. 167, 638–653. doi: 10.1016/j.petrol.2018.04.046. [CrossRef] [Google Scholar]
  • Zhang X., Jeffrey R.G. (2012) Fluid-driven multiple fracture growth from a permeable bedding plane intersected by an ascending hydraulic fracture, J. Geophys. Res. 117, B12402. doi: 10.1029/2012JB009609. [Google Scholar]
  • Zhou J., Jin Y., Chen M. (2010) Experimental investigation of hydraulic fracturing in random naturally fractured blocks, Int. J. Rock Mech. Min. Sci. 47, 7, 1193–1199. doi: 10.1016/j.ijrmms.2010.07.005. [CrossRef] [Google Scholar]
  • Zhou J., Zhang L., Pan Z., Han Z. (2017) Numerical studies of interactions between hydraulic and natural fractures by Smooth Joint Model, J. Nat. Gas Sci. Eng. 46, 592–602. doi: 10.1016/j.jngse.2017.07.030. [CrossRef] [Google Scholar]

All Tables

Table 1

Mechanical properties of synthetic sample and natural fracture filling materials (Fatahi et al., 2017).

In the text
Table 2

Rock mechanical properties of rock mass, FDFN, and NFs.

In the text

All Figures

thumbnail Fig. 1

Failure modes in tessellation media: a) tensile and shear stresses induced by compression force, b) the idealized relation between hydraulic aperture and normal stress for a rock joint, c) tensile rupture mode, and d) shear rupture mode.
In the text
thumbnail Fig. 2

Generated fracture network before and after processing by common DEM (La Pointe, 1988).
In the text
thumbnail 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.
In the text
thumbnail Fig. 4

Numerical model dimensions for investigating the effect of joint sets on HF.
In the text
thumbnail Fig. 5

The hydraulic opening of fictitious joints and joint set in times of 400, 800, 1200, 1600, 2000, 2400, and 2500 s in hydrostatic condition and θ = 45°.
In the text
thumbnail Fig. 6

The HFP path in the model of K = 1.5 and different approaching angles at the injection time of approximately 2400 s.
In the text
thumbnail Fig. 7

Interaction mode between NF and HF under different conditions.
In the text
thumbnail Fig. 8

Shear displacements in the model of K = 1.5 and θ = 45°.
In the text
thumbnail Fig. 9

The opening of fractures in two conditions of θ = 0° and 90° with three different stress ratios of 1, 1.5, and 2.
In the text
thumbnail Fig. 10

General conditions of the affected area under the influence of the HF process in a schematic way.
In the text
thumbnail Fig. 11

The effect of a) joint set dip in SRV and b) maximum opening of NFs at different times in the model of K = 1.5.
In the text
thumbnail Fig. 12

The effect of discontinuity dip on uniaxial compressive strength (Jaeger and Cook, 1979) and tensile strength for shale rock (Mighani et al., 2016) and for Lyons sandstone and Pyrophylitic (Mousavi Nezhad et al., 2018).
In the text
thumbnail Fig. 13

The effect of joints set direction on the required time for injection in different conditions.
In the text
thumbnail Fig. 14

The effect of strength parameters on SRV.
In the text
thumbnail Fig. 15

HFP at 1600 s in a model of a) φ = 35° (left) and b) φ = 20° (right).
In the text
thumbnail Fig. 16

Four simulated models with different densities/intensities.
In the text
thumbnail Fig. 17

SRV graph with time for four models.
In the text
thumbnail Fig. 18

The effect of fracture density/intensity on the HF process at 500 s for different models.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.