Abstract

Unconventional resources have been successfully exploited with technological advancements in horizontal-drilling and multistage hydraulic-fracturing, especially in North America. Due to preexisting natural fractures and the presence of stress isotropy, several complex fracture networks can be generated during fracturing operations in unconventional reservoirs. Using the DVS method, a semianalytical model was created to analyze the transient pressure behavior of a complex fracture network in which hydraulic and natural fractures interconnect with inclined angles. In this model, the complex fracture network can be divided into a proper number of segments. With this approach, we are able to focus on a detailed description of the network properties, such as the complex geometry and varying conductivity of the fracture. The accuracy of the new model was demonstrated by ECLIPSE. Using this method, we defined six flow patterns: linear flow, fracture interference flow, transitional flow, biradial flow, pseudoradial flow, and boundary response flow. A sensitivity analysis was conducted to analyze each of these flow regimes. This work provides a useful tool for reservoir engineers for fracture designing as well as estimating the performance of a complex fracture network.

1. Introduction

Technological advances in horizontal-drilling and multistage hydraulic-fracturing have stimulated a boom in unconventional resource generation throughout the world, especially in North America. Because of the presence of stress isotropy and preexisting natural fractures, complicated fracture networks can be created in unconventional reservoirs when conducting stimulation treatments [1, 2] (Fisher et al., 2002; Maxwell et al., 2002). Knowledge of the fluid flow behavior in these complex fracture networks is essential information to evaluate the performance and effectiveness of stimulation.

Several models have been established by scholars to predict the behavior of fracture networks in unconventional reservoirs in the last decade. These models can be divided into three categories. The first category is the analytical method, based on the well-known dual-porosity model [3, 4] (Warren and Root, 1963; Kazemi, 1969), which is comprised of the fracture network and surrounding matrix. Analytical models [59] (Brown et al.,2011; Ozkan et al.,2011; Xu et al.,2013; Leng et al.,2014; Ting et al.,2015) were developed to investigate the transient pressure behavior of multistage fractured horizontal wells on the basis of the dual-porosity hypothesis. These analytic models have helped engineers gain a comprehensive insight into the performance of fracture networks and have provided practical tools to evaluate stimulation effectiveness. However, the fracture network is very complex, and the dual-porosity medium may not rigorously capture the details of the fracture network characteristics, such as the irregular spatial distribution, complex interconnected scenarios, and conductivity heterogeneity of the fractures.

The second category is semianalytical and was developed using the method of source function [10, 11] (Zhao et al., 2014; Pin et al.,2015). In these models, the complex fracture network can be divided into a proper number of segments. This approach allows one to focus on a detailed description of the network properties and overcomes the shortcomings of the analytical method. Yet, these semianalytical models cannot accurately simulate the fluid flow in a reservoir. These semianalytical models assume that the fluid flow in the reservoir is 2D in an infinite plate reservoir, and the source function cannot consider the effects of fracture geometry. We show an approach that can overcome these shortcomings in the model described below.

The third category is numerical simulation, in which the hydraulic and natural fractures are usually represented by high permeability refined grids [1214] (Palagi and Aziz,1994; Skoreyko and Peter,2003; Li et al.,2003). Based on a structured grid system in conventional numerical simulators [1517], Mayerhofer et al. (2010), Warpinski et al. (2008), and Cipolla (2009) simulated the production of orthogonal networks in shale gas reservoirs. These works qualitatively analyzed how the size and density, fracture conductivity, matrix permeability, and gaps in fracture networks affect the horizontal well productivity. However, the numerical methods are time-consuming and have inherent uncertainties that could cause them to be less accurate.

In view of this, there is still lack of an efficient and rigorous approach to model the flow behavior of complex fracture networks in unconventional reservoirs. The main objective of this paper is to develop a semianalytical model that can evaluate the performance of such networks more rigorously and efficiently compared to the existing methods. This new model is based on the DVS (distributed volume sources) method [18] (Valko et al. 2007), in which the fractures in the network are explicitly represented by discrete segments to concentrate on the details of the network characteristics, such as the complex geometry and varying conductivity. The DVS model can simulate fluid flow in a closed boundary reservoir in 3D. The accuracy of the new model was demonstrated by ECLIPSE. Then, using this method, we defined the flow patterns of the fluid in the reservoir and conducted a sensitivity analysis of the transient pressure behavior.

2. Establishment of the Theoretical Model

First, we describe the fracture network used to introduce our approach and then provide the mathematic model for the reservoir and fracture flow.

2.1. Physical Model

For convenience, natural fractures were assumed to develop along the main hydraulic fractures orthogonally, as shown in Figure 1. The other assumptions were as follows: (1)The reservoir is anisotropic and homogeneous with closed boundaries (shaped like a box)(2)The horizontal fractured well produces at a constant rate(3)To simplify the flow model of natural fractures, which is in the middle of two main fractures, the bisector of the distance between two hydraulic fractures is a no-flow boundary (labeled in green in Figure 1). The red arrows (in Figure 1) represent the flow directions in natural fractures, which are in the middle of two main fractures. Fluid flows from the natural fractures into the hydraulic fractures, then through the main fractures into wellbore(4)The flow model inner fractures are 1D(5)The impacts of gravity and capillary forces are neglected(6)The orthogonal fracture network is symmetric about the horizontal fractured well.

From Figure 1, we can see that the total number of hydraulic fractures is and the number of natural fractures is . The distance of the main fractures is and of the natural fractures is . The half-length of the hydraulic fracture is and of the natural fracture is . The green dashed lines represent no-flow boundaries in accordance with assumption (3). Part of the fracture network is divided in the dashed box and enlarged for clarity. Each main fracture is separated into segments, and the natural fracture is discretized into segments. Thus, the total number of segments in the fracture network is . Each of the segments can be considered a little fracture, so that the pressure response caused by the complex fracture network can be calculated by the superposition of the pressure effects (a detailed introduction is provided in Section 2.2.1).

2.2. Mathematical Model
2.2.1. Reservoir Flow

Different techniques were developed to solve the single-phase slightly compressible flow problems in porous media in which the fluid is removed or injected from the oil well. One of the most widely used methods is the source function presented by Gringarten and Ramey [19]. From then on, the source function approach was successfully applied to accurately evaluate the performance of a vertical well, horizontal well, horizontal well with hydraulic fractures, and so on [2023](Cinco-Ley et al., 1981; Guppy et al., 1982; Ozkan,1988; Chen et al., 1997). The major disadvantage of this method is the inherent singularity of the solution wherever the point source is placed. Valko et al. (2007) established the DVS method to remove this limitation by assuming a source not in the form of a point but in the form of a rectilinear volume extended inside the surrounding rectilinear porous media. The DVS method has the capability to handle complicated well/fracture configurations. However, the major weakness of the DVS method is the difficulty in calculating the behavior of the complex fracture network because the inner “source box” must be parallelized with the reservoir boundaries. In the following, this shortcoming is eliminated by a new DVS function. The new DVS can calculate the pressure response when the surfaces of the fracture have inclined angles with the reservoir boundary’s directions.

(1) Valko’s (2007) DVS Model. First, the main principle of Valko’s (2007) DVS method will be introduced. A schematic of Valko’s (2007) DVS is shown in Figure 2. The reservoir is homogeneous with closed boundaries (box-shaped). The inner “source box” is assumed to be a smaller rectilinear box with surfaces parallel to the reservoir boundaries (for simplicity, the inner “source box” can be considered a little fracture).

As shown in Figure 2, the geometry of this model is described with the following parameters: reservoir dimensions in the , , and directions (, , and , respectively); reservoir permeability along the principal axes (, , and ); coordinates of the center point of the source (, , and ); and half-lengths of the source in the , , and directions (, , and , respectively). The mathematical model for the volume source in closed boundaries was established by Valko et al. (2007) and is shown in Appendix A.

The pressure response of a rectilinear reservoir with closed boundaries for an instantaneous withdrawal from the volume source was given by Valko (2007) and is as follows (the definitions of the dimensionless variables are shown as Appendix B):

Using Equation (1), we obtained a 3D solution of the instantaneous pressure response in anisotropic reservoirs where the permeability along the three axes is different from each other (, , and ). In addition, contrary to Gringarten’s source function, Equation (1) can take the dimension of volume source (, , and ) into account.

(2) Establishment of the New DVS Model. The presence of volume source surfaces in the plane, which are not parallel to the reservoir boundaries, is a common occurrence in complex fracture networks. The schematic for this case is shown in Figure 3.

The inclined angle between the volume source and the -axis is denoted as . The endpoint coordinates of the inclined source are (, ) and (, ). The coordinate of the center line for a particular volume source (labeled blue dashed line in Figure 3) is a constant when . At other angles, the coordinate of the center line is a function of . As shown in Equation (A.4), cx is a term in the Heaviside unit-step function. The role of the Heaviside unit-step function in Equation (A.1) is to limit the position and geometry of the volume source. Therefore, the solution form for the volume source surfaces that are not parallel to reservoir boundaries is the same as that in Equation (1), with the addition of a formula to calculate (shown as Equation (2)).

The geometrical relationship between and is presented as follows.

Figure 4 shows the case in which the volume source forms an angle with the -axis.

The geometrical relationship between and is similar to Equation (2) and is as follows.

Therefore, the pressure response for an instantaneous withdrawal from a volume source, which can have arbitrarily angle (in or both and directions), is given as follows:

From Equations (4) and (5) generalized the forms of Valko’s solution, which is specific condition when is presented.

To obtain the response of the reservoir for a continuous volume source, we numerically integrate the pressure derivative solution over time:

As for the complex fracture network, shown in Figure 1, each segment can be considered a volume source; therefore, the total number of volume sources is

The pressure response at any point in the reservoir (or inside the fracture network) can be calculated by superimposing the source function of the segments. Thus, the dimensionless pressure of fracture can be obtained as follows.

In Equation (8), represents the source strength of the segment , and represents the dimensionless pressure calculated at the center of segment if the source is placed in segment .

Applying Equation (8) to all of the fracture segments in the fracture network, equations are obtained.

2.2.2. Fracture Flow Model

Following assumptions (3) and (4), the flow model was established for the fracture network. According to (Gringarten, et al. 1974; Cinco-Ley et al., 1988) [24, 25], the diffusivity equation in the fractures can be described as follows.

The initial condition is

The inner boundary condition is

The outer boundary condition is

With the definitions for the dimensionless variables (see Appendix B), the above equations can be adapted as follows:

Solving Equation (13) with the initial condition and boundary conditions, pressure drawdown in the fracture can be obtained as follows:

Discretizing Equation (17) gives the dimensionless pressure drawdown in each fracture segment as follows, for hydraulic fracture segments:

In Equation (22), , , and .

Based on assumption (3), fluid flows from the natural fractures into the hydraulic fractures and then through the main fractures into the wellbore. Therefore, the expression for is shown as follows:

In Equation (19), represents the flux into the natural fracture segments.

For the natural fracture segments, we have the following:

In Equation (20), , , and .

Therefore, other equations are obtained. Considering the continuity condition along with the fracture face, the flux and pressure should satisfy the following condition:

The assumption of constant production rate requires

Currently we have obtained equations from Equations (8), (18), (20), and (22). Similarly, there are unknowns, including , , and . Using the Gauss-Jordan elimination, the can be obtained by simultaneously solving the system of equations.

3. Results and Discussion

3.1. Model Validation

The accuracy of this new model was verified using ECLIPSE (Schlumberger 2010). The orthogonal fracture network for the simulation included four main hydraulic fractures and four natural fractures (similar to Figure 1). The grid scale in the simulation was , and the volume of the reservoir was . Details for the parameters used in the calculations are summarized in Table 1.

The results of the semianalytical model and ECLIPSE are compared in Figure 5 (where represents the derivative of the dimensionless pressure ). From the figure, it is observed that there is a good agreement between our model and the ECLIPSE results.

Type curves of the transient pressure behavior for the complex fracture network are shown in Figure 6. From Figure 6, we can see that six main flow regimes can be recognized as follows:

Regime I is linear flow. As is commonly known, the typical feature of this flow behavior is that the slope of the dimensionless derivative pressure is equal to 0.5.

Regime II is a relatively rare occurrence in the literature. Figure 6 shows that a “cave” occurs at the end of linear flow. Few published papers have discussed this phenomenon, although published work shows this “cave” phenomenon (Pin et al., 2015). To further examine the “cave” behavior, we conducted calculations to analyze this phenomenon in Section 2.2.1 (Figure 7). The results showed that this “cave” reflects the effects of interference between hydraulic fractures and natural fractures. Therefore, we denote this process as “fracture interference flow.”

Regime III is transitional flow, generally raised at the end of Regime II.

Regime IV is biradial flow, which can be recognized by a one-third slope of the dimensionless derivative pressure. This regime has been observed by several other researchers (Zhao et al.,2013; Luo et al.,2014; Chen et al.,2015) [2628].

Regime V is pseudoradial flow with the derivative pressure stabilized at a value of 0.5.

Regime VI is reservoir boundary response flow. In this stage, transient pressure has spread to the outer closed boundaries. The dimensionless derivative pressure curve tilted up and converged to a straight line with unit slope.

3.2. Effect of Fracture Permeability
3.2.1. Effect of Varying the Natural Fracture’s Permeability

Fracture permeability is a key parameter for fractured wells. Here, the effect of fracture permeability on the behavior of transient pressure and fluid flow regimes is evaluated. We considered the permeability of the hydraulic fractures to be constant and varied the natural fracture’s conductivity. The combinations of and are shown in Table 2.

The effect of natural fracture permeability on transient pressure response of the four formations is shown in Figure 7 (the following results are based on data from Table 1).

From Figure 7, we can see that natural fracture permeability primarily affects the fracture interference flow. Larger permeability values for the natural fracture corresponded to a deeper “cave” in the typical curve. This is a result of the increasing permeability of natural fractures, causing larger fluid flow into natural fractures, leading to a stronger interference between hydraulic fractures and natural fractures. This flow regime is a typical signature of transient pressure behavior in complex fracture network.

3.2.2. The Permeability of the Complex Fracture Network Is Homogeneous

Assuming that the proppant is evenly distributed throughout the network, it is suggested that the permeability of the complex network is homogeneous. The data in Table 1 were used in the following calculations. Four cases that were investigated in which permeability of the fracture network were 8D, 18D, 28D, and 58D. Figure 8 shows the effect of fracture permeability on the transient response behavior of the complex network.

From Figure 8, we can see that the permeability of the complex fracture network only influences the pressure response in the early stages of the process.

3.3. Effect of the Complex Fracture network’s Geometry

The geometry of the fracture network also has an important influence on the pressure response in unconventional reservoirs. For simplicity, it was assumed that all of the fractures have the same length. Five cases were considered in which the fracture lengths were equal to 200 m, 300 m, 500 m, 700 m, and 900 m individually. These values were chosen to investigate the effect of the complex fracture network’s geometry on the pressure behavior. The following results are based on the data in Table 1.

The results of Figure 9 indicate that the geometry of the fracture network primarily affects transitional flow, biradial flow, and pseudoradial flow. It had no significant effect on other flow regimes. With the increase of fracture length, the period of transitional flow is also increased, and the biradial flow and pseudoradial flow gradually faded out.

3.4. Unorthogonal Fracture Network

In this section, the transient pressure behavior of unorthogonal fracture networks is investigated. Figure 10 shows a sketch of an unorthogonal fracture network used in the following calculation.

Figure 10 shows an unorthogonal complex fracture network that is composed of several fracture segments with arbitrary angles. From the dashed box in Figure 10, we can see that the discretization is not perfect in the connection of two fractures. There are gaps in the connections, because the surface of the volume source must be a parallelogram. However, we assumed that the flow inside the fracture network was continuous. The permeability of hydraulic fracture was set to 40D, and that of natural fracture was 10D. The total length of the fracture network was 1000 m, and the other pertinent data are listed in Table 1. The semianalytic model results were verified by the results of ECLIPSE (Schlumberger 2010), shown in Figure 11.

The comparison between numerical results and this new model is shown in Figure 11. The difference of the new model compared with ECLIPSE is relatively large in the early period. The reason for this difference may be a result of assumption (3) and the imperfect connections between two fractures in our model, as mentioned above. Except for early timed, the semianalytical model matches very well with the numerical results.

4. Conclusions

This paper provides a semianalytical model for the transient pressure behavior in unconventional reservoirs that have a complex fracture network. The model is capable of simulating the complex fracture network with varying conductivities and complex geometry. Although most of the results and discussion have been restricted to an orthogonal fracture network, we have demonstrated that the approach can be used for unorthogonal fracture networks in which hydraulic and natural fractures interconnect with arbitrary inclined angles by direct comparison with numerical results. The following conclusions can be drawn: (1)We present a more flexible DVS model based on the work by Valko. Then a semianalytic model was established to describe transient pressure behavior of complex fracture networks in unconventional reservoirs with closed boundaries. The accuracy of the new model was demonstrated by comparison with numerical results. In addition, the model used 3D flow to simulate the reservoir flow (in Section 2.2.1)(2)The process of fluid flow in unconventional reservoirs with complex network can be divided into six flow regimes: linear flow, fracture interference flow, transitional flow, biradial flow, pseudoradial flow, and boundary response flow. Note that the “fracture interference flow” is a new flow regime that requires additional work to more fully describe it. Through the research in Sections 2.2 and 3.3, we determined that the permeability of the complex fracture networks has a significant influence on the fracture interference flow regime(3)The results of a sensitivity analysis show that the permeability of the fractures significantly influences earlier stage fluid flow (linear flow and fracture interference flow). The geometry of the fracture network primarily affects transitional flow, biradial flow, and pseudoradial flow (shown in Figure 9).

As shown in Figure 11, the model deviates slightly from the numerical results in unorthogonal fracture networks. However future work will be focused on the optimization of the model for this case. Even so, the model is a useful tool to investigate the flow behavior of complex fracture networks. With this essential knowledge, we can evaluate well performance and stimulation effectiveness in unconventional reservoirs.

Appendix

A. The diffusivity equation for the volume source model is given by Valko et al. (2007) as follows.

where

The initial condition is

The closed boundary conditions are

in Equation (A.1) is the source function which, for the instantaneous volume source, is written as

In Equation (A.4), and represent the Dirac delta function and the Heaviside unit-step function, respectively. The Dirac delta function makes the source instantaneity, and the Heaviside unit-step function limits the geometry of the source.

B. The dimensionless parameters are defined as follows:

Nomenclature

:Length of reservoir in direction, m
:Width of reservoir in direction, m
:Height of reservoir in direction, m
:Length, m
:Width, m
:Coordinate of midpoint of volume source, m
:Permeability, m2
:Total number of hydraulic fractures, integer
:Total number of natural fractures, integer
:Total segments of one hydraulic fracture, integer
:Total segments of one natural fracture, integer
:Wellbore pressure, Pa
:Pressure, Pa
:Instantaneous pressure, Pa
:Total number of fracture segments, integer
:Viscosity, MPas
:Height of fracture in direction, m
:Flow rate,
:Time, second.
Subscripts
D:Dimensionless
N:Natural fracture
H:Hydraulic fracture
f:Fracture
, :Fracture segments number, integer.

Data Availability

Data is available when required.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work is funded by the State Key Laboratory of Shale Oil and Gas Enrichment Mechanisms and Effective Development, National Natural Science Foundation of China (Grant No. 51804258, 51974255, and 51874241), Natural Science Basic Research Program of Shaanxi Province (Grant No. 2019JQ-807, 2020JM-544, 2018JM-5054), and the Youth Innovation Team of Shaanxi Universities.