# Transient Pressure Behavior of Complex Fracture Networks in Unconventional Reservoirs

Authors: Gou Feifei; Liu Chuanxi; Ren Zongxiao; Qu Zhan; Wang Sukai; Qin Xuejie; Fang Wenchao; Wang Ping; Wang Xinzhu
Journal: Geofluids
 (2021)

---

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


---

## Body

## 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 [5–9] (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 [12–14] (Palagi and Aziz,1994; Skoreyko and Peter,2003; Li et al.,2003). Based on a structured grid system in conventional numerical simulators [15–17], 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 Figure1. 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 Figure1). 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.Figure 1
Physical model and discretized fracture segments of a complex fracture network.From Figure1, we can see that the total number of hydraulic fractures is TH and the number of natural fractures is TN. The distance of the main fractures is ΔLH and of the natural fractures is ΔLN. The half-length of the hydraulic fracture is LH and of the natural fracture is LN. 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 SH segments, and the natural fracture is discretized into SN segments. Thus, the total number of segments in the fracture network is TH×SH+TH×TN×SN. 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 [20–23](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).Figure 2
Schematic of the volume source in closed boundaries.As shown in Figure2, the geometry of this model is described with the following parameters: reservoir dimensions in the x, y, and z directions (xe, ye, and ze, respectively); reservoir permeability along the principal axes (kx, ky, and kz); coordinates of the center point of the source (cx, cy, and cz); and half-lengths of the source in the x, y, and z directions (wx, wy, and wz, 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 AppendixB):
(1)pδDxD,yD,zD,tD=1+∑n=1∞sinnπcxD+wxD−sinnπcxD−wxD2nπwxDcosnπxDe−n2π2kx/kL/xe2tD×1+∑n=1∞sinnπcyD+wyD−sinnπcyD−wyD2nπwyDcosnπyDe−n2π2ky/kL/ye2tD×1+∑n=1∞sinnπczD+wzD−sinnπczD−wzD2nπwzDcosnπzDe−n2π2kz/kL/ze2tD.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 (kx, ky, and kz). In addition, contrary to Gringarten’s source function, Equation (1) can take the dimension of volume source (2wx, 2wy, and 2wz) into account.(2) Establishment of the New DVS Model. The presence of volume source surfaces in the x−y 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.Figure 3
The surface of the volume source is not parallel to the reservoir boundaries.The inclined angle between the volume source and thex-axis is denoted as θx. The endpoint coordinates of the inclined source are (cx1, cy1) and (cx2, cy2). The cx coordinate of the center line for a particular volume source (labeled blue dashed line in Figure 3) is a constant when θx=90o. At other angles, the cx coordinate of the center line is a function of y. 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 cx (shown as Equation (2)).The geometrical relationship betweencx and y is presented as follows.
(2)cx=cx1+y−cy1/tanθx,θx≠π2cx1,θx=π2.Figure4 shows the case in which the volume source forms an angle θy with the y-axis.Figure 4
The surface of volume source is not parallel to the reservoir boundaries.The geometrical relationship betweency and x is similar to Equation (2) and is as follows.
(3)cy=cy1+x−cx1/tanθy,θy≠π2cy1,θy=π2.Therefore, the pressure response for an instantaneous withdrawal from a volume source, which can have arbitrarily angle (inx/y or both x and y directions), is given as follows:
(4)pδDxD,yD,zD,tD=1+∑n=1∞sinnπcxD+wxD−sinnπcxD−wxD2nπwxDcosnπxDe−n2π2kx/kL/xe2tD×1+∑n=1∞sinnπcyD+wyD−sinnπcyD−wyD2nπwyDcosnπyDe−n2π2ky/kL/ye2tD×1+∑n=1∞sinnπczD+wzD−sinnπczD−wzD2nπwzDcosnπzDe−n2π2kz/kL/ze2tD,(5)cxD=cxD1+yD−cyD1/tanθx,θx≠π2cxD1,θx=π2,cyD=cyD1+xD−cxD1/tanθy,θy≠π2cyD1,θy=π2.From Equations (4) and (5) generalized the forms of Valko’s solution, which is specific condition when θ=90o is presented.To obtain the response of the reservoir for a continuous volume source, we numerically integrate the pressure derivative solution over time:(6)pDxD,yD,zD,tD=qD∫0tDpδDxD,yD,zD,τdτ.As for the complex fracture network, shown in Figure1, each segment can be considered a volume source; therefore, the total number of volume sources is
(7)NT=TH×SH+TH×TN×SN.The pressure response at any point in the reservoir (or inside the fracture network) can be calculated by superimposing the source function of theNT segments. Thus, the dimensionless pressure of fracture i can be obtained as follows.
(8)pDi=∑j=1n=NTqDjpDi,j.In Equation (8), qDj represents the source strength of the segment j, and pDi,j represents the dimensionless pressure calculated at the center of segment i if the source is placed in segment j.Applying Equation (8) to all of the fracture segments in the fracture network, NT 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.
(9)kfμ∂2pf∂x2+qfwfh=0.The initial condition is(10)pft=0=pi.The inner boundary condition is(11)pfx=0=pw.The outer boundary condition is(12)∂pf∂xx=LH=0.With the definitions for the dimensionless variables (see AppendixB), the above equations can be adapted as follows:
(13)−∂2pfD∂xD2+LhCfDqfD=0,(14)pfDtD=0=0,(15)pfDxD=0=pwD,(16)∂pfD∂xDxD=LH=0.Solving Equation (13) with the initial condition and boundary conditions, pressure drawdown in the fracture can be obtained as follows:
(17)pfDxD−pwD=LhCfD∫0xD∫0ξqfDududγ−LhCfDqfwDxD.Discretizing Equation (17) gives the dimensionless pressure drawdown in each fracture segment as follows, for hydraulic fracture segments:
(18)pfDixDi−pwD=LhCfDjΔxD8QfDi+∑j=1i−1QfDiΔxD2+xDi−jΔxD−xDi∑j=1SHQfDi.In Equation (22), i=1,2,⋯,TH×SH, ΔxD=LFDH/SH, and xDi=i−1/2ΔxD.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 forQfDi is shown as follows:
(19)QfD=qfD+∑k=1SNqfDkLfDNLfDH,forcrosssegmentsqfD,forindependentsegments.In Equation (19), qfDk represents the flux into the natural fracture segments.For the natural fracture segments, we have the following:(20)pfDixDi−pwD=LhCfDjΔxD8qfDi+∑j=1i−1qfDiΔxD2+xDi−jΔxD−xDi∑j=1SNqfDi.In Equation (20), i=1,2,⋯,TH×TN×SN, ΔxD=LFDN/SN, and xDi=i−1/2ΔxD.Therefore, otherNT equations are obtained. Considering the continuity condition along with the fracture face, the flux and pressure should satisfy the following condition:
(21)pfDi=pDi,qfDi=qDi.The assumption of constant production rate requires(22)∑i=1NTqfDi=1.Currently we have obtained2NT+1 equations from Equations (8), (18), (20), and (22). Similarly, there are 2NT+1 unknowns, including pwD, qDi, and pDi. Using the Gauss-Jordan elimination, the pwD can be obtained by simultaneously solving the system of equations.


## 2.1. Physical Model

For convenience, natural fractures were assumed to develop along the main hydraulic fractures orthogonally, as shown in Figure1. 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 Figure1). 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.Figure 1
Physical model and discretized fracture segments of a complex fracture network.From Figure1, we can see that the total number of hydraulic fractures is TH and the number of natural fractures is TN. The distance of the main fractures is ΔLH and of the natural fractures is ΔLN. The half-length of the hydraulic fracture is LH and of the natural fracture is LN. 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 SH segments, and the natural fracture is discretized into SN segments. Thus, the total number of segments in the fracture network is TH×SH+TH×TN×SN. 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 [20–23](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).Figure 2
Schematic of the volume source in closed boundaries.As shown in Figure2, the geometry of this model is described with the following parameters: reservoir dimensions in the x, y, and z directions (xe, ye, and ze, respectively); reservoir permeability along the principal axes (kx, ky, and kz); coordinates of the center point of the source (cx, cy, and cz); and half-lengths of the source in the x, y, and z directions (wx, wy, and wz, 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 AppendixB):
(1)pδDxD,yD,zD,tD=1+∑n=1∞sinnπcxD+wxD−sinnπcxD−wxD2nπwxDcosnπxDe−n2π2kx/kL/xe2tD×1+∑n=1∞sinnπcyD+wyD−sinnπcyD−wyD2nπwyDcosnπyDe−n2π2ky/kL/ye2tD×1+∑n=1∞sinnπczD+wzD−sinnπczD−wzD2nπwzDcosnπzDe−n2π2kz/kL/ze2tD.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 (kx, ky, and kz). In addition, contrary to Gringarten’s source function, Equation (1) can take the dimension of volume source (2wx, 2wy, and 2wz) into account.(2) Establishment of the New DVS Model. The presence of volume source surfaces in the x−y 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.Figure 3
The surface of the volume source is not parallel to the reservoir boundaries.The inclined angle between the volume source and thex-axis is denoted as θx. The endpoint coordinates of the inclined source are (cx1, cy1) and (cx2, cy2). The cx coordinate of the center line for a particular volume source (labeled blue dashed line in Figure 3) is a constant when θx=90o. At other angles, the cx coordinate of the center line is a function of y. 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 cx (shown as Equation (2)).The geometrical relationship betweencx and y is presented as follows.
(2)cx=cx1+y−cy1/tanθx,θx≠π2cx1,θx=π2.Figure4 shows the case in which the volume source forms an angle θy with the y-axis.Figure 4
The surface of volume source is not parallel to the reservoir boundaries.The geometrical relationship betweency and x is similar to Equation (2) and is as follows.
(3)cy=cy1+x−cx1/tanθy,θy≠π2cy1,θy=π2.Therefore, the pressure response for an instantaneous withdrawal from a volume source, which can have arbitrarily angle (inx/y or both x and y directions), is given as follows:
(4)pδDxD,yD,zD,tD=1+∑n=1∞sinnπcxD+wxD−sinnπcxD−wxD2nπwxDcosnπxDe−n2π2kx/kL/xe2tD×1+∑n=1∞sinnπcyD+wyD−sinnπcyD−wyD2nπwyDcosnπyDe−n2π2ky/kL/ye2tD×1+∑n=1∞sinnπczD+wzD−sinnπczD−wzD2nπwzDcosnπzDe−n2π2kz/kL/ze2tD,(5)cxD=cxD1+yD−cyD1/tanθx,θx≠π2cxD1,θx=π2,cyD=cyD1+xD−cxD1/tanθy,θy≠π2cyD1,θy=π2.From Equations (4) and (5) generalized the forms of Valko’s solution, which is specific condition when θ=90o is presented.To obtain the response of the reservoir for a continuous volume source, we numerically integrate the pressure derivative solution over time:(6)pDxD,yD,zD,tD=qD∫0tDpδDxD,yD,zD,τdτ.As for the complex fracture network, shown in Figure1, each segment can be considered a volume source; therefore, the total number of volume sources is
(7)NT=TH×SH+TH×TN×SN.The pressure response at any point in the reservoir (or inside the fracture network) can be calculated by superimposing the source function of theNT segments. Thus, the dimensionless pressure of fracture i can be obtained as follows.
(8)pDi=∑j=1n=NTqDjpDi,j.In Equation (8), qDj represents the source strength of the segment j, and pDi,j represents the dimensionless pressure calculated at the center of segment i if the source is placed in segment j.Applying Equation (8) to all of the fracture segments in the fracture network, NT 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.
(9)kfμ∂2pf∂x2+qfwfh=0.The initial condition is(10)pft=0=pi.The inner boundary condition is(11)pfx=0=pw.The outer boundary condition is(12)∂pf∂xx=LH=0.With the definitions for the dimensionless variables (see AppendixB), the above equations can be adapted as follows:
(13)−∂2pfD∂xD2+LhCfDqfD=0,(14)pfDtD=0=0,(15)pfDxD=0=pwD,(16)∂pfD∂xDxD=LH=0.Solving Equation (13) with the initial condition and boundary conditions, pressure drawdown in the fracture can be obtained as follows:
(17)pfDxD−pwD=LhCfD∫0xD∫0ξqfDududγ−LhCfDqfwDxD.Discretizing Equation (17) gives the dimensionless pressure drawdown in each fracture segment as follows, for hydraulic fracture segments:
(18)pfDixDi−pwD=LhCfDjΔxD8QfDi+∑j=1i−1QfDiΔxD2+xDi−jΔxD−xDi∑j=1SHQfDi.In Equation (22), i=1,2,⋯,TH×SH, ΔxD=LFDH/SH, and xDi=i−1/2ΔxD.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 forQfDi is shown as follows:
(19)QfD=qfD+∑k=1SNqfDkLfDNLfDH,forcrosssegmentsqfD,forindependentsegments.In Equation (19), qfDk represents the flux into the natural fracture segments.For the natural fracture segments, we have the following:(20)pfDixDi−pwD=LhCfDjΔxD8qfDi+∑j=1i−1qfDiΔxD2+xDi−jΔxD−xDi∑j=1SNqfDi.In Equation (20), i=1,2,⋯,TH×TN×SN, ΔxD=LFDN/SN, and xDi=i−1/2ΔxD.Therefore, otherNT equations are obtained. Considering the continuity condition along with the fracture face, the flux and pressure should satisfy the following condition:
(21)pfDi=pDi,qfDi=qDi.The assumption of constant production rate requires(22)∑i=1NTqfDi=1.Currently we have obtained2NT+1 equations from Equations (8), (18), (20), and (22). Similarly, there are 2NT+1 unknowns, including pwD, qDi, and pDi. Using the Gauss-Jordan elimination, the pwD can be obtained by simultaneously solving the system of equations.


## 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 [20–23](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).Figure 2
Schematic of the volume source in closed boundaries.As shown in Figure2, the geometry of this model is described with the following parameters: reservoir dimensions in the x, y, and z directions (xe, ye, and ze, respectively); reservoir permeability along the principal axes (kx, ky, and kz); coordinates of the center point of the source (cx, cy, and cz); and half-lengths of the source in the x, y, and z directions (wx, wy, and wz, 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 AppendixB):
(1)pδDxD,yD,zD,tD=1+∑n=1∞sinnπcxD+wxD−sinnπcxD−wxD2nπwxDcosnπxDe−n2π2kx/kL/xe2tD×1+∑n=1∞sinnπcyD+wyD−sinnπcyD−wyD2nπwyDcosnπyDe−n2π2ky/kL/ye2tD×1+∑n=1∞sinnπczD+wzD−sinnπczD−wzD2nπwzDcosnπzDe−n2π2kz/kL/ze2tD.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 (kx, ky, and kz). In addition, contrary to Gringarten’s source function, Equation (1) can take the dimension of volume source (2wx, 2wy, and 2wz) into account.(2) Establishment of the New DVS Model. The presence of volume source surfaces in the x−y 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.Figure 3
The surface of the volume source is not parallel to the reservoir boundaries.The inclined angle between the volume source and thex-axis is denoted as θx. The endpoint coordinates of the inclined source are (cx1, cy1) and (cx2, cy2). The cx coordinate of the center line for a particular volume source (labeled blue dashed line in Figure 3) is a constant when θx=90o. At other angles, the cx coordinate of the center line is a function of y. 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 cx (shown as Equation (2)).The geometrical relationship betweencx and y is presented as follows.
(2)cx=cx1+y−cy1/tanθx,θx≠π2cx1,θx=π2.Figure4 shows the case in which the volume source forms an angle θy with the y-axis.Figure 4
The surface of volume source is not parallel to the reservoir boundaries.The geometrical relationship betweency and x is similar to Equation (2) and is as follows.
(3)cy=cy1+x−cx1/tanθy,θy≠π2cy1,θy=π2.Therefore, the pressure response for an instantaneous withdrawal from a volume source, which can have arbitrarily angle (inx/y or both x and y directions), is given as follows:
(4)pδDxD,yD,zD,tD=1+∑n=1∞sinnπcxD+wxD−sinnπcxD−wxD2nπwxDcosnπxDe−n2π2kx/kL/xe2tD×1+∑n=1∞sinnπcyD+wyD−sinnπcyD−wyD2nπwyDcosnπyDe−n2π2ky/kL/ye2tD×1+∑n=1∞sinnπczD+wzD−sinnπczD−wzD2nπwzDcosnπzDe−n2π2kz/kL/ze2tD,(5)cxD=cxD1+yD−cyD1/tanθx,θx≠π2cxD1,θx=π2,cyD=cyD1+xD−cxD1/tanθy,θy≠π2cyD1,θy=π2.From Equations (4) and (5) generalized the forms of Valko’s solution, which is specific condition when θ=90o is presented.To obtain the response of the reservoir for a continuous volume source, we numerically integrate the pressure derivative solution over time:(6)pDxD,yD,zD,tD=qD∫0tDpδDxD,yD,zD,τdτ.As for the complex fracture network, shown in Figure1, each segment can be considered a volume source; therefore, the total number of volume sources is
(7)NT=TH×SH+TH×TN×SN.The pressure response at any point in the reservoir (or inside the fracture network) can be calculated by superimposing the source function of theNT segments. Thus, the dimensionless pressure of fracture i can be obtained as follows.
(8)pDi=∑j=1n=NTqDjpDi,j.In Equation (8), qDj represents the source strength of the segment j, and pDi,j represents the dimensionless pressure calculated at the center of segment i if the source is placed in segment j.Applying Equation (8) to all of the fracture segments in the fracture network, NT 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.
(9)kfμ∂2pf∂x2+qfwfh=0.The initial condition is(10)pft=0=pi.The inner boundary condition is(11)pfx=0=pw.The outer boundary condition is(12)∂pf∂xx=LH=0.With the definitions for the dimensionless variables (see AppendixB), the above equations can be adapted as follows:
(13)−∂2pfD∂xD2+LhCfDqfD=0,(14)pfDtD=0=0,(15)pfDxD=0=pwD,(16)∂pfD∂xDxD=LH=0.Solving Equation (13) with the initial condition and boundary conditions, pressure drawdown in the fracture can be obtained as follows:
(17)pfDxD−pwD=LhCfD∫0xD∫0ξqfDududγ−LhCfDqfwDxD.Discretizing Equation (17) gives the dimensionless pressure drawdown in each fracture segment as follows, for hydraulic fracture segments:
(18)pfDixDi−pwD=LhCfDjΔxD8QfDi+∑j=1i−1QfDiΔxD2+xDi−jΔxD−xDi∑j=1SHQfDi.In Equation (22), i=1,2,⋯,TH×SH, ΔxD=LFDH/SH, and xDi=i−1/2ΔxD.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 forQfDi is shown as follows:
(19)QfD=qfD+∑k=1SNqfDkLfDNLfDH,forcrosssegmentsqfD,forindependentsegments.In Equation (19), qfDk represents the flux into the natural fracture segments.For the natural fracture segments, we have the following:(20)pfDixDi−pwD=LhCfDjΔxD8qfDi+∑j=1i−1qfDiΔxD2+xDi−jΔxD−xDi∑j=1SNqfDi.In Equation (20), i=1,2,⋯,TH×TN×SN, ΔxD=LFDN/SN, and xDi=i−1/2ΔxD.Therefore, otherNT equations are obtained. Considering the continuity condition along with the fracture face, the flux and pressure should satisfy the following condition:
(21)pfDi=pDi,qfDi=qDi.The assumption of constant production rate requires(22)∑i=1NTqfDi=1.Currently we have obtained2NT+1 equations from Equations (8), (18), (20), and (22). Similarly, there are 2NT+1 unknowns, including pwD, qDi, and pDi. Using the Gauss-Jordan elimination, the pwD 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 Figure1). The grid scale in the simulation was 161×161×10, and the volume of the reservoir was 1610×1610×10m3. Details for the parameters used in the calculations are summarized in Table 1.Table 1
Data used for semianalytical and numerical models.
Reservoir length, m1610Permeability of hydraulic fractures, D28Reservoir width, m1610Length of natural fractures, m200Reservoir thickness, m10Width of natural fractures, m0.0028Reservoir compressibility, MPa-10.00022Permeability of natural fractures, D28Reservoir permeability, mD0.1Wellbore radius, m0.15Fluid viscosity, pas0.001Production rate, m3/d0.5Length of hydraulic fractures, m200Reservoir initial pressure, MPa14.6Width of hydraulic fractures, m0.0028Reservoir porosity0.1The results of the semianalytical model and ECLIPSE are compared in Figure5 (where dpD represents the derivative of the dimensionless pressure pD). From the figure, it is observed that there is a good agreement between our model and the ECLIPSE results.Figure 5
Comparison of the semianalytic model and numerical results.Type curves of the transient pressure behavior for the complex fracture network are shown in Figure6. From Figure 6, we can see that six main flow regimes can be recognized as follows:Figure 6
Flow regimes of the complex fracture network.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. Figure6 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.”Figure 7
Comparison of dimensionless derivative pressure for the four reservoirs.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) [26–28].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 ofkfH and kfN are shown in Table 2.Table 2
The combinations ofkfH and kfN.
ReservoirsPermeability of hydraulic fractures, DPermeability of natural fractures, DFormation A281.4Formation B284.4Formation C288.8Formation D2828The effect of natural fracture permeability on transient pressure response of the four formations is shown in Figure7 (the following results are based on data from Table 1).From Figure7, 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 Table1 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.Figure 8
Comparison of dimensionless derivative pressure for four cases.From Figure8, 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 Table1.The results of Figure9 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.Figure 9
Comparison of type curves with different fracture geometry.


### 3.4. Unorthogonal Fracture Network

In this section, the transient pressure behavior of unorthogonal fracture networks is investigated. Figure10 shows a sketch of an unorthogonal fracture network used in the following calculation.Figure 10
Sketch of unorthogonal fracture network.Figure10 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.Figure 11
Comparison between numerical results and semianalysis results.The comparison between numerical results and this new model is shown in Figure11. 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.


## 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 Figure1). The grid scale in the simulation was 161×161×10, and the volume of the reservoir was 1610×1610×10m3. Details for the parameters used in the calculations are summarized in Table 1.Table 1
Data used for semianalytical and numerical models.
Reservoir length, m1610Permeability of hydraulic fractures, D28Reservoir width, m1610Length of natural fractures, m200Reservoir thickness, m10Width of natural fractures, m0.0028Reservoir compressibility, MPa-10.00022Permeability of natural fractures, D28Reservoir permeability, mD0.1Wellbore radius, m0.15Fluid viscosity, pas0.001Production rate, m3/d0.5Length of hydraulic fractures, m200Reservoir initial pressure, MPa14.6Width of hydraulic fractures, m0.0028Reservoir porosity0.1The results of the semianalytical model and ECLIPSE are compared in Figure5 (where dpD represents the derivative of the dimensionless pressure pD). From the figure, it is observed that there is a good agreement between our model and the ECLIPSE results.Figure 5
Comparison of the semianalytic model and numerical results.Type curves of the transient pressure behavior for the complex fracture network are shown in Figure6. From Figure 6, we can see that six main flow regimes can be recognized as follows:Figure 6
Flow regimes of the complex fracture network.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. Figure6 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.”Figure 7
Comparison of dimensionless derivative pressure for the four reservoirs.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) [26–28].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 ofkfH and kfN are shown in Table 2.Table 2
The combinations ofkfH and kfN.
ReservoirsPermeability of hydraulic fractures, DPermeability of natural fractures, DFormation A281.4Formation B284.4Formation C288.8Formation D2828The effect of natural fracture permeability on transient pressure response of the four formations is shown in Figure7 (the following results are based on data from Table 1).From Figure7, 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 Table1 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.Figure 8
Comparison of dimensionless derivative pressure for four cases.From Figure8, we can see that the permeability of the complex fracture network only influences the pressure response in the early stages of the process.


## 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 ofkfH and kfN are shown in Table 2.Table 2
The combinations ofkfH and kfN.
ReservoirsPermeability of hydraulic fractures, DPermeability of natural fractures, DFormation A281.4Formation B284.4Formation C288.8Formation D2828The effect of natural fracture permeability on transient pressure response of the four formations is shown in Figure7 (the following results are based on data from Table 1).From Figure7, 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 Table1 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.Figure 8
Comparison of dimensionless derivative pressure for four cases.From Figure8, 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 Table1.The results of Figure9 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.Figure 9
Comparison of type curves with different fracture geometry.


## 3.4. Unorthogonal Fracture Network

In this section, the transient pressure behavior of unorthogonal fracture networks is investigated. Figure10 shows a sketch of an unorthogonal fracture network used in the following calculation.Figure 10
Sketch of unorthogonal fracture network.Figure10 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.Figure 11
Comparison between numerical results and semianalysis results.The comparison between numerical results and this new model is shown in Figure11. 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 Section2.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 Sections2.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 Figure9).As shown in Figure11, 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.



---

*Source: 6273822-2021-11-03.xml*