Next Article in Journal
Modeling and Vulnerability Analysis of Cyber-Physical Power Systems Considering Network Topology and Power Flow Properties
Next Article in Special Issue
Near Wellbore Hydraulic Fracture Propagation from Perforations in Tight Rocks: The Roles of Fracturing Fluid Viscosity and Injection Rate
Previous Article in Journal
Integrated Equivalent Circuit and Thermal Model for Simulation of Temperature-Dependent LiFePO4 Battery in Actual Embedded Application
Previous Article in Special Issue
Numerical Simulation of the Propagation of Hydraulic and Natural Fracture Using Dijkstra’s Algorithm
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Investigation of the Time-Dependent and the Proppant Dominated Stress Shadow Effects in a Transverse Multiple Fracture System and Optimization

1
State Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing University, 400030 Chongqing, China
2
Energy Research Center of Lower Saxony, 38640 Goslar, Germany
*
Author to whom correspondence should be addressed.
Energies 2017, 10(1), 83; https://doi.org/10.3390/en10010083
Submission received: 21 October 2016 / Revised: 8 December 2016 / Accepted: 3 January 2017 / Published: 11 January 2017

Abstract

:
In this paper, a numerical study is conducted to investigate the stress shadow effects (stress reorientation and change) during hydraulic fracturing in a transverse multiple fracture system. A numerical model is used for the numerical study. It is a 3D model and can simulate the fracture operation from injection begin to full closure (fracture contact). Therefore, there is no need to assume the fracture geometry for the investigation of the stress shadow effects (unlike previous studies). In the numerical study, the first and second operations in a fictive transverse multiple fracture system are simulated, meanwhile the stress shadow effects and their influences on the propagation and proppant placement of the second fracture are investigated. According to the results, the following conclusions are discerned: (1) most proppants are located in the lower part of the reservoir, even below the perforation; (2) the stress shadow effects are time-dependent and proppant dominated; (3) the stress shadow effects affect the fracture propagation and the proppant placement of the second fracture, and also the fracture conductivity of the first fracture; (4) the time-dependent stress shadow effects can be divided into four phases, fracture enlargement, closure without proppant contact, closure with proppant contact and full closure; and (5) the superposition effect of the stress shadow in a transverse multiple fracture system exists. According to the conclusions, some optimizations are recommended.

1. Introduction

Hydraulic fracturing is currently a hot topic in petroleum engineering. It is now the key technology to enhance the gas recovery in tight and shale gas reservoirs, especially when combined with horizontal well bores. During fracturing operations, reservoirs are artificially cracked. With the support of injected solid proppant, a huge highly conductive channel is then created; thus, gas can be economically produced. On the other hand, a horizontal well bore makes it possible to set a fracture cluster in a reservoir; then a large stimulated reservoir volume (SRV) is created.
At present, many researchers focus on the optimization of the transverse multiple fracture system. Meyer et al. [1] developed an analytical solution for predicting the behavior of multiple patterned transverse vertical hydraulic fractures intercepting horizontal wells. Yu and Sepehrnoori [2] used a 2D semi-analytical model to investigate the influences of reservoir permeability, porosity, fracture spacing, fracture half length, fracture conductivity, gas desorption and well spacing on the productivity and the related economic analysis. Later, in Yu and Sepehrnoori [3], the influence of geomechanical effects on the fracture conductivity was further considered. In these 2D models, each fracture in the transverse multiple fracture system has an equivalent geometry (length) and conductivity (aperture), as well as an equivalent fracture spacing to its neighboring fracture. The fracture spacing is not limited, which is unrealistic, especially when the fracture spacing is short, because the influences of the neighboring fractures during fracturing operation are significantly great. In Olson [4], it was found that the fracture spacing and operational sequence (simultaneous, sequential or alternative sequential fracturing) have great influences on the fracture pattern. In a simultaneous operation, the fractures in the middle have normally a shorter length, while, in a sequential operation, the influences of the previous fractures may cause reorientation of the subsequent fractures.
As an optimization concept for the transverse multiple fracture system, man-made fractures with limited reorientation and limited propagation resistance from neighboring fractures are generally preferred, hence, the influences of each fracture on the others must be studied in detail. In fact, the main influence is the stress shadow effect caused by stress changes (magnitude and orientation) due to fracturing operations. The stress shadow effect was first reported in Green et al. [5] and Sneddon et al. [6]. In these publications, the stress distribution of a semi-infinite 2D fracture under constant internal pressure was studied. Soliman and Addams [7] used analytical solutions for a semi-infinite 2D and a penny-shaped fracture to investigate the stress shadow effect of a transverse multiple fracture system. It was found that the stress shadow effect increases with an increasing number of sequential fractures and a decreasing fracture spacing. Cheng [8] used a 2D model based on the boundary element method to analyze the stress distributions around multiple transverse fractures (simultaneously operated) and the geometries of these fractures. It was found that the width profile is no longer elliptic, and the width of the middle fractures is strongly reduced. Nagel and Sanchez-Nagel [9] used commercial software based on the finite difference method and the distinct element method to investigate the stress shadow effect in 3D. The influences of the fracture spacing, rock mechanical properties and in-situ stress ratio were studied. Roussel and Sharma [10] used the same method as that in the work of Nagel and Sanchez-Nagel [9], however, the influence of fracture barriers (cap rock and base rock) was further considered. In that paper, it was found that an alternate fracture sequence (1-3-2-4-5, etc.) is better for creating a multiple fracture system with limited reorientation. Rafiee [11] presented an optimization concept named the modified zipper fracture. The concept is based on the reduction of the stress shadow effect calculated by an analytical solution and a 3D numerical model. In Kresse et al. [12] and Wu et al. [13], a new hydraulic fracturing model for simulating complex fracture network propagation in a formation with pre-existing natural fractures was introduced. The stress shadow effect is calculated using a 2D, plane-strain, displacement discontinuity solution with a 3D correction factor.
Although numerous studies of the stress shadow effect in a transverse multiple fracture system have been reported, there are still some problems. In these studies, an assumed final fracture pattern (semi-infinite or penny-shaped) with an assumed geometry (length, height and width) is generally used to investigate the stress shadow effect, implying that the maximum stress shadow effect has the same depth as the middle point of the fracture. Indeed, the real fracture at the closure stage does not have an ideal geometry even in a homogeneous reservoir because the final proppant distribution dominates the final fracture geometry (Zhou et al. [14]). Due to proppant settling, the propped area below the perforation is normally greater than the area above. Therefore, the magnitude and the shape of the stress shadow effect could be much different than that predicted in the previous studies. Secondly, the stress shadow effect varies from the operation stage to the closure stage.
To the author’s knowledge, there is little information available in the literature about the time dependency of the stress shadow effect and the influences of the proppant placement on the stress shadow effect. The understanding of the changing process is greatly important for the spacing, sequence, as well as operation schedule optimization etc. of a transverse multiple fracture system. The objective of this work is to investigate the temporal variation of the stress shadow effect in consideration of proppant transport and placement. A full 3D numerical model is applied, which can describe not only the fracturing operation process from injection to full closure but also its influences on the surrounding rock formations. According to the results, some optimization concepts are recommended.

2. Theoretical Background

Generally, hydraulic fracturing involves the following main physical processes: mechanical deformation induced by fluid pressure changes in fractures and pores; fluid flow within the fracture and the rock formation, including their interactions; fracture propagation; and proppant transport and settling inside a fracture. In the mathematical model, the rock deformation is described by linear elasticity theory, including the equation of motion (Equation (1)), the continuum equation (Equation (2)) and the constitutive equations (Equations (3) and (4), derived from Darcy law and mass conservation, is used to describe the porous media flow process. The mathematical model of the fracture flow (Equation (5)) is similar to the one for porous media flow. The flow equations for the two processes have the same formulation but different characteristic parameters. In the fracture flow, the conductivity is derived from the simplified Navier–Stokes equation for two parallel plates. For the proppant transport, the mass conservation equation (Equation (6)) is used, in which the settling is implicitly expressed in terms of the proppant velocity. The numerical formulation is described in Appendix A.
σ i j , j + ρ ( b i d v i d t ) = 0 ,
Δ ε i j = 1 2 ( Δ u i , j + Δ u j , i ) ,
Δ σ = D Δ ε ,
Q l e a k + α e t + 1 M b P p t = [ K m μ ( P p + ρ f g Z ) ] ,
w ( Q i n j e . Q l e a k ) + w t = [ w 3 12 μ ( P f + ρ f g Z ) ] ,
( C w ) t + ( C v p w ) + w C i n j e . Q i n j e . = 0 ,
where σ = σ′ + αIPp and σ is the total stress (Pa); ρ is the rock density (kg/m3); bi is the volumetric acceleration (m/s2); vi is rock the mass velocity (m/s); Δε is the strain increment (unitless); u is the displacement (m); Δσ′ is effective the stress increment (Pa); α is the Biot coefficient (unitless); I is the unit matrix; D is the physical matrix, i, j ∈ (x,y,z); Qleak is the leak off rate (1/s); Qinje. is the fluid injection rate (1/s); Km is the permeability (m2); e is the volumetric strain (unitless); t is the time (s); Mb is the Biot modulus (Pa); μ is the viscosity (Pa*s); Pp is the pore pressure (Pa); Pf is the fracture pressure (Pa); ρf is the fluid density (kg/m3); g is the gravity (m/s2); C is the proppant concentration (%); w is the fracture width (m); vp is the proppant velocity (m/s); and Cinje. is the proppant injection concentration (unitless).
Most hydraulic fracturing models and commercial software used today, for instance FracPro based on a Pseudo 3D model (Cleary [15] and Cleary et al. [16]), MFrac based on a modified Pseudo 3D model (Meyer [17]) and StimPlan3D based on a Planar 3D model (Peirce et al. [18] and Siebrits et al. [19]) do not consider all of the physical processes mentioned above, e.g., fluid leak off, pore fluid flow, and stress variation in the surrounding rocks. Furthermore, they cannot describe the contact conditions (fracture–fracture contact and fracture–proppant contact) under full fracture closure.
Zhou [14,20] developed a full 3D model based on the governing equations. The simulations of hydraulic fracturing from fracture propagation to full closure (including fracture contact) and proppant transport were verified through comparison of a semi-analytical solution and the results of a proppant transport experiment. Therefore, it can be used to investigate the time-dependent and the proppant dominated stress shadow effect in a transverse multiple fracture system. In the model, these hydro-mechanical coupled equations are numerically solved in a sequential schema (Figure 1), in which the coupling terms have been modified to maintain the conservation of the fluid mass in the fracture. In a computational loop of this method, Equations (4) and (5) are solved independently using the finite volume method at first. In the numerical formulation (Appendix A), the variation of the fracture width in Equation (5) is substituted by a developed constitutive equation (Equation (7)), which describe the aperture variation induced by pressure and stress variation normally acting on the fracture plane Zhou et al. [14]. Then, the fluid pressure and velocity in the fracture and pores as well as the fracture width are estimated. A contact condition is introduced to determine the fracture contact (Equation (7)). If the fracture width is smaller than the given residual width or the proppant concentration is greater than the given maximum value, then the fracture surfaces make contact with either themselves or the proppant, and a contact stress is calculated through the over-reduced fracture width.
Δ w = P f ( t + 1 ) + σ c o n ( t ) + σ n ( t ) α 1 l c ,
where lc is characteristic length of a fracture element (m); α1 = K + 4G/3 is the elastic parameter (MPa); K is bulk modulus (MPa); G is the shear modulus (MPa); σn is the normal stress (Pa); σcon is the contact stress { σ c o n ( t + 1 ) = 0 if   C C max   and   w w r e s σ c o n ( t + 1 ) = σ c o n ( t ) + α 1 Δ ε o if   C > C max   or   w < w r e s (Pa); Cmax is the maximum proppant concentration (unitless); wres is the residual width (m); and Δεo is the over reduced strain (unitless).
With the calculated fluid velocity, as well as the proppant settling and relative velocity models, the proppant velocity can be computed according Equation (8) (Zhou et al. [14]). Then, the distribution of the proppant concentration can be estimated through solving Equation (6). The final step is the mechanical calculation. After substituting the fracture width and the pore pressure into the mechanical model, the new stress, strain and displacement fields will be computed (see Appendix A).
v p = k w c v l + v p s t c c o r r w c o r r cos θ = k w c v l + v p s cos θ ,
where vl is the liquid velocity (m/s), vp is the proppant velocity (m/s), kwc is the relative velocity coefficient (unitless), vpst is the corrected Stoke’s settling velocity for a single particle (m/s), ccorr is the influence factor of proppant concentration (unitless), and wcorr is the influence factor of well retardation (unitless).

3. Numerical Investigation

In this paper, we conduct a numerical study of a fictive transverse multiple fracture system. The main objective is to investigate the influence of time and proppant distribution on the stress shadow effect. Therefore, other influence factors such as rock mechanical properties and in-situ stress ratio are not concerned.

3.1. Modelgeneration

A half reservoir model is generated and illustrated in Figure 2. It has a geometry of 200 m (x) × 550 m (y) × 200 m (z) and consists of three horizontal layers: cap rock, pay zone and base rock. The mechanical and hydraulic properties of each layer are listed in Table 1. Figure 3 shows the primary stress and pore pressure. The direction of the minimum horizontal stress is the same as that of the y-axis. A horizontal well is drilled through the middle part of the model and along the direction of the minimum horizontal stress, indicating a transverse multiple fracture system.

3.2. Numerical Results of the First Fracturing Operation

The first injection position (0, 250, −3100) is set 300 m away from the left side of the reservoir model (Figure 2). An 80 min operation with an injection rate of 6 m3/min is designed. At t = 30 min, proppant is added with a stepwise increasing concentration from 0 to 500 kg/m3 (Figure 4). The properties of the fluid and proppant are listed in Table 2.
Figure 5 shows the temporal development of the fracture and injection volumes. The difference between them is the leak off volume. During the injection phase, the fracture volume increases linearly, then decreases from 80 min to ca. 710 min with a decreasing rate, and finally remains constant, indicating that the final proppant placement and full fracture closure are achieved after that time.
Figure 6 shows the fracture geometry in terms of the width contour, the distribution of the proppant mass per area and the proppant volume fraction. At t = 80 min, the fracture geometry is almost regular and limited between the top and the bottom barriers, which have much higher closure stress. The maximum width is located at the perforation. Unlike the fracture geometry, the proppant has an irregular distribution due to the settling effect. At t = 1000 min (920 min after shut in), the fracture width and the proppant mass per area have the same contour distribution, implying that the magnitude of the fracture width and the area of the propped fracture after full closure are dominated by the proppant placement. The maximum width is reduced from 2.53 cm to 0.95 cm and its position is 34 m deeper than the perforation. The (half) effective fracture area is also reduced from 14,998 m2 to 4196 m2.
Figure 7 shows the spatial contour of the minimum horizontal stress with 2.5° and 5° reorientations at different time. The contour has two symmetric wings and its shape changes with time. In general, the shape of the contour enlarges during the injection and shrinks after shut in. At t = 1000 min, the contour with 5° reorientation disappears. For a transverse multiple fracture system, the maximum influence point of the stress reorientation with a certain grade is defined as the point on the contour furthest to the fracture plane (Figure 7). In general, its projection point is inside the created fracture (except t = 1 min).
The position of the maximum influence point with a 2.5° stress reorientation varies with time. Figure 8 shows its temporal development. The x-coordinate increases quickly during the injection and slowly during the first 420 min after shut in (till 500 min). The tendency is similar to that of the fracture propagation (fast propagation during injection and slow post-propagation after shut in). Then, it decreases until full closure. The decrease is due to the full closure of the front region, where no proppant is placed. The y-coordinate also increases during the injection, decreases after shut in (from ca. 300 m to 262.5 m) and remains unchanged after full closure. The tendency is similar to that of the fracture volume (Figure 5), indicating that the y-coordinate of the maximum influence point is strongly dependent on the created fracture volume. The z-coordinate remains unchanged before t = 300 min because the fracture has an enlargement and closure without influence of proppant until that time. Then, it decreases rapidly followed by a slight increase. The final z-coordinate is 34 m below the perforation because of the irregular fracture closure due to the proppant redistribution and placement. As a result of the settling effect, the proppant is located in the lower part of the fracture. It prevents the created fracture from fully closing. Conversely, the upper part without proppant can be fully closed. When the proppant in the lower part contacts with the fracture, closure with proppant contact occurs.
Figure 9 shows the change of the minimum horizontal stress at t = 80 min and 1000 min. In general, the stresses increase in the region perpendicular to the opened fracture and, meanwhile, decrease outside the fracture front. The contour shape of the stress increase (i.e., 2 MPa in Figure 9) is similar to a cone whose bottom is perpendicular to the fracture plane. The cone-like contour shrinks during the closure process, implying that the influence decreases.

3.3. Numerical Results of the Second Fracturing Operation after Full Closure of the First Fracture

In a transverse multiple fracture system, the previous fracture may cause new fractures to deviate from the trajectory of the previous fracture. A minimum fracture spacing to ensue transverse multiple fracture growth is desirable for efficient gas production. To avoid strong direction deviations, a minimum fracture spacing must be characterized, i.e., the distance from the maximum influence point with a 2.5° stress reorientation to the previous fracture represents the minimum fracture spacing with a fracture reorientation smaller than 2.5°.
Two fracture spacings (25 m and 50 m) are chosen to investigate the influences of the stress shadow effect on the second operation. The maximum stress reorientation on the plane with 25 m and 50 m spacing are 1.5° and 0.8°, respectively. Due to the small stress reorientation and for simplification, the propagation direction of the second fracture is assumed to be parallel to that of the first operation. The stress shadow effect is then the change of the minimum horizontal stress shown in Figure 10.
Figure 11 shows the comparison of the fracture propagation before shut in and the propped fracture after final closure of the first and the second operations. Although the created fracture areas are similar, the propagation of the second fracture before shut in becomes irregular due to the stress shadow effect. The lower part has a smaller propagation rate and a smaller fracture width than the upper part because the minimum horizontal stress increased in the lower part and decreased in the upper part (Figure 10). This is more obvious with the 25 m fracture spacing because the stress change is greater. The finial fracture patterns are also similar, but they have different areas and width distributions (Figure 11). In general, the propped area increases (Table 3), the width decreases and its distribution becomes uniform, and the final fracture pattern shifts upward. Due to the stress shadow effect, the lower part of the second fracture has a higher closure stress and a smaller fracture width, which leads to the reduction of the proppant settling, the acceleration of the fracture closure in the lower part and finally upward-shift of the proppant location. The upward-shift effect decreases with increasing fracture spacing (25 m spacing with a 10 m upward-shift and 50 m spacing with a 4 m upward-shift).
Figure 12 shows the comparison of the spatial contours with a 2.5° stress reorientation of the first and the second operation after full closure. In general, the stress reorientation effect between the two fracture planes (inside the fracture system) is strongly reduced (25 m) or even disappears (50 m); conversely, its range of the influence on both sides of the fracture system increases. The results can be explained by the superposition principle: between the two fractures, the stress reorientation effect from the two fractures has opposite directions, whereas it has the same direction outside the fractures. The distance from the maximum influence point to the fracture plane increases from 12.5 m to 24.5 m and 20.1 m for the 25 m and 50 m spacings, respectively.
Figure 13 shows the change of the minimum horizontal stress after fully closure of the second fracture. The contour changes from asymmetric to symmetric with increasing fracture spacing. The minimum horizontal stress increases due to the superposition effect. The maximum change is on the first fracture plane. It increases from 5.82 MPa to 7.55 MPa and 6.44 MPa for 25 m and 50 m spacings, respectively, indicating that the created conductivity of the first fracture can also be influenced by the second fracture (conductivity reduction caused by increasing contact stress).

3.4. Numerical Results of the Second Fracturing Operation during Closure Process of the First Fracture

According to the results of the first operation, the stress shadow effect is time dependent. If the second fracturing operation starts during the closure process of the first fracture, then the stress shadow effect dynamically influences the second fracture. In this study, the second operation with a 25 m and 50 m spacing is set to start 250 min and 160 min after closure of the first fracture, respectively. After these times, the stress reorientation on the plane of the second fracture is smaller than 2.5° (Figure 8), which ensures a transverse multiple fracture system.
The changes of the minimum horizontal stress on the fracture plane for the 25 m and 50 m spacing are shown in Figure 14. The stress change during the closure process of the first fracture is different from the one after full closure. In general, the stress in the whole pay-zone increases. The maximum stress increase for both is ca. 15 MPa. As mentioned in Section 3.2, the first fracture begins to irregularly close after 300 min (220 min after shut in). Therefore, the stress change at t = 240 min (160 min after shut in) seems to be regular, and at t = 330 min (250 min after shut in), it seems to be irregular.
Figure 15 and Figure 16 show comparisons of the fracture propagation before shut in and the propped fracture after final closure of the first and the second operations with 25 m and 50 m spacings, respectively. The time dependent influences of the stress shadow effect are clearly obtained. In general, the fracture propagation is hindered due to the increasing minimum horizontal stress in the whole pay zone. The finial fracture area for 25 m and 50 m spacings are 13,622 m2 and 14,598 m2 (Table 4), respectively, which are smaller than those of the second operation after full closure of the first fracture. The hindering effect decreases with increasing spacing. The width contour also changes. Due to the hindering effect, a high fluid pressure is needed for fracture propagation, which leads to an increase of the fracture width.
The propped areas of the cases differ. With the 25 m spacing, the propped area did not change significantly. Although the length is shorter due to the hindering effect, the height is greater. The growth in the height occurs mainly in the lower part of the proppant placement. This is due to the higher settling velocity caused by the increasing fracture width. With the 50 m spacing, the propped area is less than that of the second operation after full closure of the first fracture, but the fracture pattern and area are similar to those of the first fracture, indicating that the stress shadow effect does not have a significant influence.

4. Discussion

Through the numerical investigation of the hydraulic fracturing in a fictive reservoir, the following results are obtained:
  • Both the final fracture width and the final fracture conductivity are dominated by the final proppant placement. Although the created fracture appears to be regular, the propped fracture is irregular. In general, the upper part of the created fracture is difficult to prop, and the proppant is concentrated in the lower part due to the proppant transport and settling effects.
  • The stress shadow effect (stress reorientation and change) varies with time. Its range enlarges during the injection and shrinks during the closure process. The range depends on the fracture width distribution. Therefore, it is also dominated by the proppant placement after final closure. The maximum influence point moves firstly along with the fracture front and away from the fracture plane during the injection. Then, it moves backward to the fracture plane and the propped fracture tip during the closure process; meanwhile, it moves down to a location even below the perforation due to the placement of the proppant at the bottom of the fracture.
  • The stress shadow effect affects the fracture propagation and the proppant placement of the second fracture. The influences depend on the time and the fracture spacing. In general, the stress shadow effect decreases with increasing spacing. The time dependent influence can be divided into four time periods: fracture enlargement, closure without proppant contact, closure with proppant contact and full closure. These periods can be determined by the position trajectory of the maximum influence point with a certain reorientation grade.
  • After full closure, the second fracture propagates more in the upper part of the reservoir. The final fracture and propped area increase, whereas the width decreases and its distribution becomes uniform. The fracture pattern shifts upward. During closure without proppant contact, the fracture propagation is strongly hindered, which causes both an increase of the fracture width and the settling effect. Therefore, there is little upward-shift effect, and the propped fracture length as well as the propped area may be reduced. The influences of the stress shadow effect during closure with proppant contact are between those without proppant contact and after full closure.
  • Due to the superposition principle, the stress reorientation effect between the two fractures is reduced or even disappears, whereas it increases on both sides of the fracture system. Conversely, the change of the minimum horizontal stress also increases, especially between the two factures. The second fracture also has influences on the first fracture in the form of a conductivity reduction caused by increased contact stress on the proppant in the first fracture.
Based on these results, the following optimizations are recommended:
  • Because of the settling effect, most proppant is located in the lower part of the reservoir. Therefore, it is better to drill the horizontal bore in the lower part of the reservoir, whereby a good connectivity between the perforation and the propped fracture will be created. The precondition is that the lower barrier must have enough resistance to prevent the fracture from going through, especially for a thin pay zone.
  • The spacing optimization of the transverse multiple fracture is mainly based on the stress shadow effect (stress reorientation and change). In general, a minimum spacing must be determined, that will not result in any strong reorientation (i.e., 5° in Roussel et al. [10]) or great increase of the contact stress on the previous fractures.
  • In a sequential transverse multiple fracture system, the next operation is preferred to begin after full closure of the previous fracture because the area of the hydraulic fracture and the propped fracture are greater, and the fracture width becomes uniform.
  • Regarding proppant transport in boreholes, a horizontal borehole with a small inclination is better than one with no inclination because the gravitational force is positive. In a sequential transverse multiple fracture system, the upward-shift effect plays a positive role for an inclined borehole (Figure 17). The trajectory of the horizontal bore can be optimized according to the calculated fracture pattern.
  • If the reservoir is thin and the fracture spacing is short; then the upward-shift and the superposition effect will reduce the resistance of the upper barrier; thus, it may lose its function as a height control. To avoid this situation, a suitable spacing must be designed because the upward-shift effect decreases with increasing spacing.
  • In a simultaneous multiple fracture system, the propagation of each fracture will be hindered because the stress shadow effect during the injection is greatly stronger than after full closure. To minimize the stress shadow effect, the designed spacing is generally greater than that of a sequential system. This indicates that a simultaneous multiple fracture system can save operation time, but the fracture number must be reduced.
  • An alternative fracture sequence (1-3-2, etc.) is recommended in Roussel et al. [10]. With this sequence, the stress reorientation between the first und the third fracture is minimized as shown in Figure 12. Therefore, there is no reorientation problem for the second fracture, and, thus, the fractures can be placed much closer to each other. However, the stress change must also be considered. There are two problems with such an alternative sequence. One is the increasing stress on the second fracture plane. It hinders the propagation of the second fracture in the horizontal direction and may cause a loss of height control when the reservoir is thin and the stress difference between the pay zone and the barrier is small. On the other hand, it is difficult for a horizontal well to cross the middle part of each fracture due to the double upward-shift effect (Figure 18).
The stress shadow effects contain stress reorientation and variation. In this paper, the temporal development of the stress shadow effects during the operation and after shut in was investigated in detail. As the developed simulator is limited only for fractures with the orientation perpendicular to the minimum horizontal stress, it is not possible to simulate the fracture reorientation due to stress changes. In fact, a 3D model for simulation of arbitrary propagation during hydraulic fracturing is still a challenge. The arbitrary reorientation was already realized in 2D models, e.g., the model developed by Kresse et al. [21] based on the DDM theory. In Kresse et al. [21], a simultaneous multiple fracture operation was modeled. On both sides, a small reorientation was observed. The reorientation is probably due to the small stress anisotropy (1.4 MPa), which is not normal, especially for a deep reservoir. Zhou et al. [22] developed a 2D model using the XFEM theory. In his study, the fracture reorientation does not seem to be significant. On the other hand, the effect of stress variation is great both in Kresse and Zhou’s simulation. In this study, it was also found that the fracture propagation and proppant placement are strongly influenced by the stress variation depending on the fracture geometry of the previous fracture. In future work, efforts should be paid to develop a more realistic 3D model considering fracture reorientation, thus the 3D reorientation effect can be studied.

5. Conclusions

In this paper, a numerical study was conducted to investigate the time dependency of the stress shadow effect and the influences of the proppant placement on the stress shadow effect. According to the study, the following conclusions are discerned: (1) The magnitude and range of the stress shadow effect varies with variation of the fracture width. The variation of the fracture width has four phases, enlargement, closure without proppant contact, closure with proppant contact and full closure, revealing that the stress shadow effect is time-dependent; (2) The closure with proppant contact and full closure phases is influenced by proppant transport and placement, implying that the stress shadow effect is proppant dominated in that two phases; (3) The stress shadow effects affect the fracture propagation and the proppant placement of the next fracture, and also the fracture conductivity of the previous fracture; (4) The superposition effect of the stress shadow in a transverse multiple fracture system exists. According to the conclusions, some optimizations are recommended (Section 4).

Acknowledgments

The work presented in this paper is funded by the National Natural Science Foundation of China (Grant no. NSFC51504043) and the National Basic Research Program of China (No. 2014CB239206).

Author Contributions

Yang Gou proposed the idea; Lei Zhou contributed the calculation program and performed the simulation; Junchao Chen analyzed the data; and Wentao Feng wrote the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Numerical Formulation of the Developed Model to Simulate Hydraulic Fracturing in Tight Gas Reservoirs

In the developed numerical model, hexahedron elements and quadrilateral elements are used to discretize the rock matrix and the rock fracture, respectively. The mechanical solution of the commercial software FLAC3D (Cundall [23]) was adopted in this paper. The motion for the elastic continuum (Equation (1)) can be transformed into a discrete form of Newton’s law at each grid points (for more details, see Cundall [23]):
m l ( d v i d t ) l = T a i l 3 + m l b i + f i l = F i l   for   the   l - th   grid   point ,
where v is the velocity (m/s); l is the index of grid point; i indicates the coordinate direction; Td is the inertial force derived from stresses (N); f is the external force (N); F is the total force acting on a grid point (N); and m is the mass belongs to a grid point (kg).
In order to get a quasi-static solution, the mechanical damping should be considered. Therefore, a damping force is extra added in Equation (A1), and then:
m l ( d v i d t ) l = F i l + F d i l   for   the   l - th   grid   point ,
where F d i l = α | F i l | s i g n ( v i l ) is the damping force (N); α is the damping coefficient and smaller than one (unitless); and s i g n ( y ) = { + 1 if   y > 0 1 if   y < 0 0 if   y = 0 .
Equation (A2) can be solved using an explicit finite difference method (FDM):
v i l ( t + Δ t ) = v i l ( t ) + Δ t m l ( F i l ( t ) + F d i l ( t ) )   for   the   l - th   grid   point ,
With the calculated node velocity at t + Δt, the variation of the displacement in Δt can be estimated, which will be further used to calculate the strain and the stress variation in this step according to Equations (2) and (3), respectively.
To solve the fluid flow in the fracture, the finite volume method (FVM) was adopted. In this method, the pressure is averaged at the center of one control volume. The first derivative of pressure is approximated by a center difference schema and the time derivative by a forward difference schema. Through the finite difference approximation, Equation (5) in combination of Equation (7) can be rewritten in form of a linear equation system, which is solved using iteration methods:
L l ( t ) P f l ( t + Δ t ) + L r ( t ) P f r ( t + Δ t ) + L b ( t ) P f b ( t + Δ t ) + L t ( t ) P f t ( t + Δ t ) ( L l ( t ) + L r ( t ) + L b ( t ) + L t ( t ) + l c α 1 Δ t ) P f o ( t + Δ t ) = Q s ( σ n ( t ) + σ c o n ( t ) ) l c α 1 Δ t
where L i ( t ) = w i ( t ) 2 12 μ a s i is the transmissivity between two neighbor elements; l, r, b, and t are the indexes for the four neighbor elements; si is the distance between the center of two neighbor elements (m); and Qs is the source (m/s).
Then, the pressure and the fluid velocity field are obtained. The driven forces for proppant movement results from fluid velocity and gravity. In Zhou et al. [14], a mathematical model to describe the proppant velocity was introduced (Equation (8)), which considered the influence of the proppant concentration and the ratio of the proppant diameter to the fracture width. For solving the proppant transport (Equation (6)), a hybrid finite volume and upwind schema was applied. Equation (6) is then numerically linearized into Equation (A5).
w r ( t ) v p r ( t + 1 ) C r ( t + 1 ) w l ( t ) v p l ( t + 1 ) C l ( t + 1 ) Δ s l r + w t ( t ) v p t ( t + 1 ) C t ( t + 1 ) w b ( t ) v p b ( t + 1 ) C b ( t + 1 ) Δ s t b + C o ( t + 1 ) w o ( t + 1 ) C o ( t ) w o ( t ) Δ t + C i n j q s = 0
where { v p r C r ( t + 1 ) = C o ( t + 1 ) [ v p r , 0 ] C r ( t + 1 ) [ v p r , 0 ] v p l C l ( t + 1 ) = C l ( t + 1 ) [ v p l , 0 ] C o ( t + 1 ) [ v p l , 0 ] v p t C t ( t + 1 ) = C o ( t + 1 ) [ v p t , 0 ] C t ( t + 1 ) [ v p t , 0 ] v p b C b ( t + 1 ) = C b ( t + 1 ) [ v p b , 0 ] C o ( t + 1 ) [ v p b , 0 ] , [ x , 0 ] = { x x > 0 0 x 0 .
The transport equation of the porous flow (Equation (4)) is also linearized through the FVM in an implicit (Equation (A6)), thus the pore pressure can be obtained.
l = 1 l = n L l P p ( t + 1 ) l ( l = 1 l = n L l + 1 M b Δ t ) P p ( t + 1 ) o = Q s p ( t + 1 ) P p ( t ) o M b Δ t + α e ( t + 1 ) e ( t ) Δ t

References

  1. Meyer, B.R.; Bazan, L.W.; Jacot, R.H.; Lattibeaudiere, M.G. Optimization of multipile transverse hydraulic fractures in horizontal wellbores. In Proceedings of the SPE Unconventional Gas Conference, Pittsburgh, PA, USA, 23–25 February 2010.
  2. Yu, W.; Sepehrnoori, K. Optimization of multiple hydraulically fractured horizontal wells in unconventional gas reservoirs. J. Pet. Eng. 2013, 2013, 151898. [Google Scholar]
  3. Yu, W.; Sepehrnoori, K. Simulation of gas desorption and geomechanic effects for unconventional gas reservoirs. Fuel 2014, 116, 455–464. [Google Scholar] [CrossRef]
  4. Olson, J.E. Multi-fracture propagation modeling: Applications to hydraulic fracturing in shales and tight gas sands. In Proceedings of the 42nd US Rock mechanics Symposium (USRMS), San Francisco, CA, USA, 29 June–2 July 2008.
  5. Green, A.E.; Sneddon, I.N. The distribution of stress in the neighborhood of a flat elliptical crack in an elastic solid. Math. Proc. Camb. 1950, 46, 159–163. [Google Scholar] [CrossRef]
  6. Sneddon, I.N.; Elliot, H.A. The Opening of a Griffith Crack under internal Pressure. Quart. Appl. Math. 1946, 3, 262–267. [Google Scholar] [CrossRef]
  7. Soliman, M.Y.; East, L.; Adams, D. Geomechanics Aspects of multiple fracturing of horizontal and vertical wells. In Proceedings of the SPE International Thermal Operations and Heavy Oil Symposium and Western Regional Meeting, Bakersfield, CA, USA, 16–18 March 2004.
  8. Cheng, Y. Boundary Element analysis of the stress distribution around multiple fractures: Implications for the spacing of perforation clusters of hydraulically fractured horizontal wells. In Proceedings of the SPE Eastern Regional Meeting, Charleston, WV, USA, 23–25 September 2009.
  9. Nagel, N.B.; Sanchez-Nagel, M. Stress shadowing and microseismic events: A Numerical Evaluation. In Proceedings of the SPE Annual Technical Conference and Exhibition, Denver, CO, USA, 30 October–2 November 2011.
  10. Roussel, N.P.; Sharma, M.M. Optimizing fracture spacing and sequencing in horizontal well fracturing. Soc. Petr. Eng. Prod. Oper. 2011, 26, 173–184. [Google Scholar] [CrossRef]
  11. Rafiee, M.; Soliman, M.Y.; Pirayesh, E. Hydraulic fracturing Design and Optimizaiton: A Modification of Zipper Frac. In Proceedings of the SPE Annual Technical Conference and Exhibition, San Antonio, TX, USA, 8–10 October 2012.
  12. Kresse, O.; Cohen, C.; Weng, X.; Wu, R.; Gu, H. Numerical modeling of hydraulic fracturing in naturally fractured Formations. In Proceedings of the 45th U.S. Rock Mechanics/Geomechanics Symposium, San Francisco, CA, USA, 26–29 June 2011.
  13. Wu, R.; Kresse, O.; Weng, X.; Cohen, C.; Gu, H. Modeling of interaction of hydraulic fractures in complex fracture networks. In Proceedings of the SPE Hydraulic Fracturing Technology Conference, The Woodlands, TX, USA, 6–8 February 2012.
  14. Zhou, L.; Hou, M.Z.; Gou, Y.; Li, M.T. Numerical investigation of a low-efficient hydraulic fracturing operation in a tight gas reservoir in the North German Basin. J. Pet. Eng. 2014, 120, 119–129. [Google Scholar] [CrossRef]
  15. Cleary, M.P. Comprehensive Design formulae for hydraulic fracturing. In Proceedings of the SPE Annual Technical Conference and Exhibition, Dallas, TX, USA, 21–24 September 1980.
  16. Cleary, M.P.; Barr, D.T.; Willis, R.M. Enhancement of real-time hydraulic fracturing models with full 3D simulation. In Proceedings of the SPE Gas Technology Symposium, Dallas, TX, USA, 13–15 June 1988.
  17. Meyer, B.R. Three-Dimensional Hydraulic Fracturing Simulation on Personal Computers: Theory and Comparison Studies. In Proceedings of the SPE Eastern Regional Meeting, Morgantown, WV, USA, 24–27 October 1989.
  18. Peirce, A.P.; Siebrits, E. Uniform asymptotic approximations for accurate modeling of cracks in layered elastic media. Int. J. Fract. 2001, 110, 205–239. [Google Scholar] [CrossRef]
  19. Siebrits, E.; Peirce, A.P. An efficient multi-layer planar 3D fracture growth algorithm using a fixed mesh approach. Int. J. Numer. Methods Eng. 2002, 53, 691–717. [Google Scholar] [CrossRef]
  20. Zhou, L.; Hou, M.Z. A new numerical 3D-Model for simulation of hydraulic fracturing in consideration of hydro-mechanical coupling effects. Int. J. Rock Mech. Min. Sci. 2013, 60, 370–380. [Google Scholar] [CrossRef]
  21. Kresse, O.; Wend, X.; Gu, H.; We, U. Numerical modeling of hydraulic fractures interaction in complex naturally fracture formations. Rock Mech. Rock. Eng. 2013, 46, 555–568. [Google Scholar] [CrossRef]
  22. Zhou, L.; Gou, Y.; Hou, Z.; Were, P. Numerical modeling and investigation of fracture propagation with arbitrary orientation through fluid injection in tight gas reservoirs with combined XFEM and FVM. Environ. Earth Sci. 2015, 73, 5801–5813. [Google Scholar] [CrossRef]
  23. Cundall, P.A. FLAC 3D Manual: A Computer Program for Fast Lagrangian Analysis of Continua, version 4.0; ITASCA Consulting Group, Inc.: Minneapolis, MN, USA, 2008.
Figure 1. Flow chart of the computing schema in the numerical model developed by Zhou et al. [20].
Figure 1. Flow chart of the computing schema in the numerical model developed by Zhou et al. [20].
Energies 10 00083 g001
Figure 2. Illustration of the fictive reservoir model.
Figure 2. Illustration of the fictive reservoir model.
Energies 10 00083 g002
Figure 3. Primary stress and pore pressure.
Figure 3. Primary stress and pore pressure.
Energies 10 00083 g003
Figure 4. Injection fluid rate and proppant mass concentration during the operation.
Figure 4. Injection fluid rate and proppant mass concentration during the operation.
Energies 10 00083 g004
Figure 5. Temporal development of the fracture volume.
Figure 5. Temporal development of the fracture volume.
Energies 10 00083 g005
Figure 6. Fracture geometry, distribution of the proppant mass per area and the proppant volume fraction at t = 80 and t = 1000 min (on the plane y = 250 m).
Figure 6. Fracture geometry, distribution of the proppant mass per area and the proppant volume fraction at t = 80 and t = 1000 min (on the plane y = 250 m).
Energies 10 00083 g006
Figure 7. Spatial contour of the minimum horizontal stress with 2.5° and 5° reorientations at different time.
Figure 7. Spatial contour of the minimum horizontal stress with 2.5° and 5° reorientations at different time.
Energies 10 00083 g007
Figure 8. Temporal development of the position coordinates of the maximum influence point with a 2.5° reorientation.
Figure 8. Temporal development of the position coordinates of the maximum influence point with a 2.5° reorientation.
Energies 10 00083 g008
Figure 9. Change of the minimum horizontal stress at t = 80 min and t = 1000 min.
Figure 9. Change of the minimum horizontal stress at t = 80 min and t = 1000 min.
Energies 10 00083 g009
Figure 10. Change of the minimum horizontal stress on the plane for 25 and 50 m fracture spacings.
Figure 10. Change of the minimum horizontal stress on the plane for 25 and 50 m fracture spacings.
Energies 10 00083 g010
Figure 11. Comparison of the fracture propagation before shut in and the propped fracture after final closure of the first and the second operations.
Figure 11. Comparison of the fracture propagation before shut in and the propped fracture after final closure of the first and the second operations.
Energies 10 00083 g011
Figure 12. Comparison of the spatial contour with a 2.5° stress reorientation after full closure of the second fracture.
Figure 12. Comparison of the spatial contour with a 2.5° stress reorientation after full closure of the second fracture.
Energies 10 00083 g012
Figure 13. Change of the minimum horizontal stress after fully closure of the second fracture.
Figure 13. Change of the minimum horizontal stress after fully closure of the second fracture.
Energies 10 00083 g013
Figure 14. Changes of the minimum horizontal stress on the plane with 25 m and 50 m spacings after full closure at 250 min and 160 min after shut in of the first operation.
Figure 14. Changes of the minimum horizontal stress on the plane with 25 m and 50 m spacings after full closure at 250 min and 160 min after shut in of the first operation.
Energies 10 00083 g014
Figure 15. Comparison of the fracture propagation before shut in and the propped fracture after final closure of the first and the second operations with a 25 m spacing.
Figure 15. Comparison of the fracture propagation before shut in and the propped fracture after final closure of the first and the second operations with a 25 m spacing.
Energies 10 00083 g015
Figure 16. Comparison of the fracture propagation before shut in and the propped fracture after final closure of the first and the second operations with a 50 m spacing.
Figure 16. Comparison of the fracture propagation before shut in and the propped fracture after final closure of the first and the second operations with a 50 m spacing.
Energies 10 00083 g016
Figure 17. Demonstration of the borehole placement in a sequential transverse multiple fracture system.
Figure 17. Demonstration of the borehole placement in a sequential transverse multiple fracture system.
Energies 10 00083 g017
Figure 18. Schematic demonstration of a transverse multiple fracture with alternative sequence.
Figure 18. Schematic demonstration of a transverse multiple fracture with alternative sequence.
Energies 10 00083 g018
Table 1. Mechanical and hydraulic properties of the rock formations in the calculation model.
Table 1. Mechanical and hydraulic properties of the rock formations in the calculation model.
Rock FormationRock Density (ρ) (kg/m3)Young’s Modulus (E) (GPa)Poisson Ratio (υ) (unitless)Tensile Strength (σt) (MPa)Porosity (φ) (unitless)Permeability (Km) (m2)
Cap Rock2600250.31.00.0251 × 10−17
Pay Zone2600300.251.00.11 × 10−15
Base Rock2600250.31.00.0251 × 10−17
Table 2. Fluid and proppant parameters.
Table 2. Fluid and proppant parameters.
ParametersSymbolUnitValues
Proppant densityρpkg/m33500
Fluid densityρfkg/m31040
Fluid viscosityμlPa∙s0.1
Proppant diameterdpmm0.667
Power law coefficientn-1
Parameter for viscosity correctiona-1.8
Maximum proppant concentrationCmax-0.65
Table 3. Comparison of the fracture and the propped area of the first and the second operations.
Table 3. Comparison of the fracture and the propped area of the first and the second operations.
First OperationSecond Operation: Spacing 25 mSecond Operation: Spacing 50 m
Fracture area (m2)14,99815,00615,016
Propped area (m2)419643494367
Table 4. Comparison of the fracture and the propped area after full closure and during the closure process.
Table 4. Comparison of the fracture and the propped area after full closure and during the closure process.
First OperationSecond Operation: Spacing 25 m (after Fully Closure)Second Operation: Spacing 25 m (250 min after Shut In)Second Operation: Spacing 50 m (after Fully Closure)Second Operation: Spacing 50 m (160 min after Shut In)
Fracture area (m2)14,99815,00613,62215,01614,598
Propped area (m2)41964349436443674173

Share and Cite

MDPI and ACS Style

Zhou, L.; Chen, J.; Gou, Y.; Feng, W. Numerical Investigation of the Time-Dependent and the Proppant Dominated Stress Shadow Effects in a Transverse Multiple Fracture System and Optimization. Energies 2017, 10, 83. https://doi.org/10.3390/en10010083

AMA Style

Zhou L, Chen J, Gou Y, Feng W. Numerical Investigation of the Time-Dependent and the Proppant Dominated Stress Shadow Effects in a Transverse Multiple Fracture System and Optimization. Energies. 2017; 10(1):83. https://doi.org/10.3390/en10010083

Chicago/Turabian Style

Zhou, Lei, Junchao Chen, Yang Gou, and Wentao Feng. 2017. "Numerical Investigation of the Time-Dependent and the Proppant Dominated Stress Shadow Effects in a Transverse Multiple Fracture System and Optimization" Energies 10, no. 1: 83. https://doi.org/10.3390/en10010083

APA Style

Zhou, L., Chen, J., Gou, Y., & Feng, W. (2017). Numerical Investigation of the Time-Dependent and the Proppant Dominated Stress Shadow Effects in a Transverse Multiple Fracture System and Optimization. Energies, 10(1), 83. https://doi.org/10.3390/en10010083

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop