Pressure-Transient Performances of Fractured Horizontal Wells in the Compartmentalized Heterogeneous Unconventional Reservoirs

: In order to investigate pressure performance of multiple fractured horizontal wells (MFHWs) penetrating heterogeneous unconventional reservoir and avoid the high computational cost of numerical simulation, a semi-analytical model for MFHWs combining Green function solution and boundary element method has been obtained, where the reservoir is divided into di ﬀ erent homogeneous substructures and coupled at interface boundaries by plane source function in a closed rectangular parallelepiped. Hydraulic fractures are assumed uniform ﬂux and dual porosity model is used for natural fractures system. Then the model is validated by compared with analytical solution of MFHWs in a homogeneous reservoir and trilinear ﬂow model, which shows that this model can achieve high accuracy even with a small interface discretization number, and it can consider the radial ﬂow around each hydraulic fractures. Finally, the pressure responses with heterogeneous parameters of reservoirs are discussed including heterogeneous permeability, non-uniform block-length and fracture half-length distribution as well as dual porosity parameters like elastic storage ratio and crossﬂow ratio.


Introduction
With the rapid development of unconventional reservoirs, the horizontal well intercepted by hydraulic fractures has become the main technology [1,2], which can significantly increase contact areas between fractures and formation zone. Thus, the transient performances of Multiple Fractured Horizontal Wells (MFHWs) have attracted more and more scholars' attention. There have been extensive studies [3][4][5] on the MFHWs of unconventional reservoir, which generally can be classified into analytical model, numerical model and semi-analytical model.
To consider the effect of formation heterogeneity, one of the most common analytical models is the trilinear flow model with three different flow areas of outer region, inner region and hydraulic fracture, which is proposed by Brown [6]. To describe the hydraulic fractures surrounding by the stimulated area, Stalgorova [7,8] subdivided the trilinear flow regions into five flow regions model. Guo [9] extended the flow area into seven parts including various fracture distribution along the horizontal wellbore. Later, on account of the fractal nature of fractured media, the new formulations are implemented in the trilinear flow model using fractal theory [3,10,11]. The fundamental petro-physical characteristic of unconventional reservoirs is considered, but all multi-linear models are based on

Modeling of MFHWs
Firstly, we will present basic the physical model and assumption for hydraulic fractures and stimulated reservoir volume (SRV). Then, the Green's function formulation of the pressure-transient solution for a locally homogeneous reservoir substructure will be given with inner and outer boundaries.

Physical Model and Assumptions
A Multiple fractured horizontal well passing through a heterogeneous reservoir with different properties is shown in Figure 1, where each block is the subsections of the reservoir characterized by uniform average properties. Except for the interfaces between the blocks, other boundaries of the blocks are assumed to be impermeable. The basic assumptions of the model are as follows: (1) The reservoir is heterogeneous along the horizontal wellbore with n f blocks of different length x ei (i = 1, 2, 3, . . . , n f ), and there is one hydraulic fracture distributed in each block. The thickness of the stimulated reservoir volume is uniform h and width is y e . The initial pressure is assumed p i . (2) The artificial hydraulic fractures are vertical to the horizontal wellbore, with fracture half-length L fi , fracture height h fi , (i = 1, 2, 3, . . . , n f ), and uniform flux distribution assumed in each hydraulic fracture. (3) Each block of unconventional reservoir was established by Warren and Root's (1963) dual porosity model, where natural fractures and matrix follow the Darcy's flow model. For each blocks, different natural fracture permeability k fi , crossflow coefficient λ i , and elastic storage ratio ω i (i = 1, 2, 3, . . . , n f ) can be defined. (4) The fluid assumed slightly compressible with constant viscosity and compressibility. (5) Gravity and capillary effect are negligible.
Energies 2020, 13,5204 3 of 16 of the stimulated reservoir volume is uniform h and width is ye. The initial pressure is assumed pi.

Pressure Solution for a Reservoir Substructure.
Based on the Green's function, the pressure transient solution of dimensionless space and time in porous media can be written as follows [18,31,32], where the domain boundary is divided into the inner and outer boundaries, represented by BwD, and BeD, respectively: where ( , ) is the pressure drop, ( , ′ , − ) is the Green's function, ( ′ , ) is the pressure. The corresponds to well surface and e is the outer surface of the dimensionless domain, the points and ′ represent the observation and source locations in dimensionless space, respectively, and tD is dimensionless time, and defined as: where, is the fluid viscosity, k is the permeability, ℓ is the characteristic reference length, is the total compressibility, is the porosity and t is time. Then, the homogeneous Neumann boundary condition is used to eliminate partial derivative of Green function. The point flux term is obtained by Darcy's law to represent the normal pressure derivatives on the boundary surfaces. Therefore, we can rewrite Equation (1) as follows:

Pressure Solution for a Reservoir Substructure
Based on the Green's function, the pressure transient solution of dimensionless space and time in porous media can be written as follows [18,31,32], where the domain boundary is divided into the inner and outer boundaries, represented by B wD , and B eD , respectively: where ∆p(M D , t D ) is the pressure drop, G M D , M D , t D − τ is the Green's function, p M D , τ is the pressure. The B wD corresponds to well surface and B eD is the outer surface of the dimensionless domain, the points M D and M D represent the observation and source locations in dimensionless space, respectively, and t D is dimensionless time, and defined as: where, µ is the fluid viscosity, k is the permeability, is the characteristic reference length, C t is the total compressibility, φ is the porosity and t is time. Then, the homogeneous Neumann boundary condition is used to eliminate partial derivative of Green function. The point flux term is obtained by Darcy's law to represent the normal pressure derivatives on the boundary surfaces. Therefore, we can rewrite Equation (1) as follows: where, q w M D , τ , q e M D , τ are the point flux of well surface and outer surface, respectively. We divide the inner and outer boundary surface into n and m segments of uniform flux, and the source function are defined as B wD = n j=1 B wDj j = 1, 2, 3, · · · , n (4) Therefore, the Equation (3) can be rewritten as follows: Ozkan and Raghavan [21] applied the Laplace transformation into Equation (7), the variable flow rate condition and naturally fractured media can be easily implemented by Green's function. Convert the convolution integrals into algebraic expressions. Then, Equation (7) can be rewritten as The function with bar signs in Equation (8) indicates their function in Laplace transforms with respect to the time t D , and s is the Laplace-transform variable.

Source Functions for the Inner and Outer Boundaries
To obtain the transient pressure solution of fractured horizontal well in a rectangular parallelepiped, the source functions of inner and outer boundaries are evaluated by the Green's function. For the partial and full penetration of the hydraulic fractures in the reservoir, 3D plane source and 2D line source should be considered.
Corresponding to the point source in the rectangular parallelepiped, the Green's function in dual porous media can be expressed as [21] cos kπ y D y eD cos kπ y D y eD chε k,n x D1 +chε k,n x D2 ε k,n sh(ε k,n x eD ) In Equation (9), we have defined  perpendicular to the x-axis, so a planar source in the y-z plane can be calculated by integration of Equation (9) (shown in Figure 2) √ + ( / ) 2 + ( /ℎ ) 2 and = ( ), ( ) = (1− )+ + (1− ) for naturally fractured reservoir. The inner and outer boundaries consist of the planar surfaces in the rectangular parallelepiped, which represented hydraulic fractures and boundary planes. The source function for a plane is perpendicular to the x-axis, so a planar source in the y-z plane can be calculated by integration of Equation (9) (shown in Figure 2) When the reservoir is fully penetrated by hydraulic fractures, the boundary surface is at full penetration, namely ℎ = ℎ ; thus, Equation (10) can simplify to (the second item of Equation (10) is equal to 0) When the reservoir is fully penetrated by hydraulic fractures, the boundary surface is at full penetration, namely h wD = h D ; thus, Equation (10) can simplify to (the second item of Equation (10) is equal to 0)

Solution of Coupling Multiple Blocks
To obtain the pressure drop solution of MFHW in multiple reservoir blocks, we will demonstrate the coupling process of n f rectangular reservoir blocks in series in the x direction as shown in Figure 1. Then, the pressure solution process of the linear system of equations is presented.
The fractured horizontal well extends through a heterogeneous reservoir with rectangular parallelepiped blocks. All block boundary surfaces are assumed to be impermeable, except for the interface between different blocks. We will assume that the hydraulic fractures are uniform flux distribution. For the infinite conductivity condition, the pressure response can be achieved by subdividing each fracture or assumed by 0.732 point pressure of fracture, and finite conductivity fractures can be obtained by incorporating pressure drop with hydraulic fractures [24]. Here, each hydraulic fracture as one segment represents inner boundary of each block, while boundary Energies 2020, 13, 5204 6 of 15 surfaces are discretized into m e segments. From Equation (8), the pressure drop at the center of the ith fracture can be given by where where k = 1, 2, 3, . . . ., m e , M Dei,1k , M Dei,2k indicate the center point of kth segment of left and right surface in ith block, so q ei,1k q ei,2k is the according flux of point M Dei,1k and M Dei,2k , respectively. Except the first and last block of SRV, each block has 1 + 2m e points with unknown parameters of pressure and flux, therefore there are (1 + 2m e ) × 2 × n f − 4m e unknown. Based on Equations (12)-(14), we yield n f + 2 n f − 1 × m e linear equations. In order to get the same number of equations as the unknowns, the pressure and flux continuity condition can be applied.
(1) Continuity conditions of pressure and flux at the interface The pressure drop of the same interface should be equal for each segment k of different blocks, and the flux of interfaces should be continuous, so we have (2) Horizontal wellbore conditions Assuming that the horizontal wellbore is infinite conductivity, the pressure drop of each fracture mid-point M Dwi should be equal to the pressure drop of horizontal wellbore ∆p(M Dw , s), so we can have n f − 1 additional equations Energies 2020, 13, 5204 Moreover, the sum of the fluxes from all hydraulic fractures should be equal to the total production rate q total .
The linear system defined by Equations (12)- (17) has (1 + 2m e ) × 2 × n f − 4m e equations for (1 + 2m e ) × 2 × n f − 4m e unknowns, which can be written in the matrix-vector form AX = B. Here, take a multiple fractured horizontal well with n f fractures and boundary segment m e = 1; for example, the components of coefficient matrix A is where to simplify writing, s). The solution vector, X, has the following components: X = q w1 , q e1,2 , q w2 , q e2,2 , q w3 , q e3,2 , · · · · · · , q w(n f −1) , q e(n f −1),2 , q wn f , ∆p wD T 2n f ×1 (19) and the components of the right-hand side vector, B, are B = [0, 0, 0, 0, · · · · · · , 0, q total ] T (2n f )×1 (20) The pressure and flux solution of this linear system can be solved. Then, the solution in the Laplace transform domain can be inverted into real-time domain by the Stehfest Algorithm [33]. The dimensionless pressure and dimensionless fluxes of ith fracture are defined as:

Comparison and Verification
In order to verify the accuracy of our model, a relatively simple model of MFHW is discussed, which located in the center of a closed, homogenous-and single-porosity media, rectangular reservoir with three blocks. The transient pressure performance of MFHW in homogeneous reservoirs is calculated by [24] using the source function and the principle of superposition method. The properties of the hydraulic fractures and reservoir are given in Table 1, and the results with different discretized number m e (Figure 3) are shown in Table 2 and plotted in Figure 4. Table 2 presents the comparisons of our models with Zerzar's model about pressure and pressure derivatives. Four cases of different interface discretization between the blocks are shown in Figure 3. Table 2 demonstrates that there are almost no differences among the four cases. All of them can match the Zerzar's model very well. The maximum error (m e = 1) between the Zerzar's model and our model is 0.0664% of dimensionless pressure drop, 0.2043% of dimensionless pressure derivative, which is too small to be distinguished on the log-log plot of Figure 4, which shows high correctness of our model. To simplify the calculation process, we assumed m e = 1 for the following discussion. reservoir with three blocks. The transient pressure performance of MFHW in homogeneous reservoirs is calculated by [24] using the source function and the principle of superposition method. The properties of the hydraulic fractures and reservoir are given in Table 1, and the results with different discretized number me (Figure 3) are shown in Table 2 and plotted in Figure 4.   Table 2 presents the comparisons of our models with Zerzar's model about pressure and pressure derivatives. Four cases of different interface discretization between the blocks are shown in Figure 3. Table 2 demonstrates that there are almost no differences among the four cases. All of them can match the Zerzar's model very well. The maximum error (me = 1) between the Zerzar's model and our model is 0.0664% of dimensionless pressure drop, 0.2043% of dimensionless pressure derivative, which is too small to be distinguished on the log-log plot of Figure 4, which shows high correctness of our model. To simplify the calculation process, we assumed me = 1 for the following discussion.   p D dp D /dlnt D p D dp D /dlnt D p D dp D /dlnt D p D dp D /dlnt D p D dp D /dlnt D In Figure 4, the flow stages of MFHW can be identified including first linear flow, first radial flow and late pseudo-steady boundary flow, which are different from traditional MFHWs flow regimes. As hydraulic fractures are close to the impermeable boundaries in our model, the second linear and radial flows do not appear. In addition, we compare our model with trilinear flow model [6], which is widely used for unconventional reservoir because of its concise formula and convenient calculation. In Figure 4, it shows that the pressure and pressure derivative of our model can be matched perfectly with trilinear model at first linear flow and late pseudo-steady boundary flow period. However, radical flow around each hydraulic fracture is missing in trilinear flow model, which is not in line with real situation. Therefore, our model is more accurate than trilinear flow model.

Results and Discussion
The semi-analytical pressure solution for MFHWs within some rectangular reservoir blocks of different parameters is obtained. This section presents pressure behavior and flux distribution of each hydraulic fracture with some heterogeneous parameters. We assume that there are three compartments of the reservoir, the horizontal well penetrates these three compartments and there In Figure 4, the flow stages of MFHW can be identified including first linear flow, first radial flow and late pseudo-steady boundary flow, which are different from traditional MFHWs flow regimes. As hydraulic fractures are close to the impermeable boundaries in our model, the second linear and radial flows do not appear. In addition, we compare our model with trilinear flow model [6], which is widely used for unconventional reservoir because of its concise formula and convenient calculation.
In Figure 4, it shows that the pressure and pressure derivative of our model can be matched perfectly with trilinear model at first linear flow and late pseudo-steady boundary flow period. However, radical flow around each hydraulic fracture is missing in trilinear flow model, which is not in line with real situation. Therefore, our model is more accurate than trilinear flow model.

Results and Discussion
The semi-analytical pressure solution for MFHWs within some rectangular reservoir blocks of different parameters is obtained. This section presents pressure behavior and flux distribution of each hydraulic fracture with some heterogeneous parameters. We assume that there are three compartments of the reservoir, the horizontal well penetrates these three compartments and there are three hydraulic fractures in the middle of each block. Besides the basic properties given in Table 1, the parameters in dual porosity media includes elastic storage ratio ω f = 0.1 and crossflow coefficient λ f = 2. The heterogeneous parameter ratio values of three blocks are shown in Table 3.

Heterogeneous Permeability
We consider five cases of heterogeneous permeability distribution, and the natural fracture permeability ratio of these three blocks is shown in Table 3. In cases 2 and 4, the permeability increases in second and third compartment, while in cases 3 and 5, the permeability decreases in second and third compartment. Figure 5 shows the pressure and pressure derivative responses as well as the flux distribution of three different hydraulic fractures.
Energies 2020, 13, 5204 11 of 16 second and third compartment. Figure 5 shows the pressure and pressure derivative responses as well as the flux distribution of three different hydraulic fractures. In Figure 5a, for the pressure-transient responses, we can see that the heterogeneous permeability distribution of different drainage blocks affects both arriving time and the pressure drop of the whole MFHW flow regimes including first linear flow, crossflow regime, first radial flow and late pseudo-steady boundary flow. The higher the permeability of cases 2 and 4, the earlier different regimes arrive, and the less the pressure drop compared with reference case 1 (homogeneous In Figure 5a, for the pressure-transient responses, we can see that the heterogeneous permeability distribution of different drainage blocks affects both arriving time and the pressure drop of the whole MFHW flow regimes including first linear flow, crossflow regime, first radial flow and late pseudo-steady boundary flow. The higher the permeability of cases 2 and 4, the earlier different regimes arrive, and the less the pressure drop compared with reference case 1 (homogeneous permeability), which is beneficial to production. On the contrary, with the lower permeability of cases 3 and 5, the higher-pressure transient curve means more pressure drop demands for constant flow rate condition of horizontal wellbore. Figure 5b shows the effect of heterogeneous permeability on the flux of each fracture. As shown in Figure 5b, there are great differences between different fracture fluxes with unequal permeability distribution. In homogeneous permeability case 1, the flux of each fracture is equal to q total /3. For the other heterogeneous permeability cases, the flux of fractures in higher permeability blocks is larger during the early linear flow period. Then, the difference will decrease because of crossflow between each block. Finally, the differences will become stable in the last pseudo-steady boundary period. Simultaneously, for lower-permeability case 5, the crossflow period appears earlier compared to the higher-permeability case 4, which is the same as the pressure response in Figure 5a.

Non-Uniform Length of Blocks
To consider the non-uniform length distribution effect, we use a heterogeneous permeability model with k f1 :k f2 :k f3 = 1:2:1. Three different cases with unequal length distributions of each compartment are considered in Table 3 under the constant whole drainage length condition. Then, the pressure response and flux distribution are shown in Figure 6. From Figure 6a, we can see that if the whole drainage length is constant, the dimensionless pressure drop and pressure derivative curves have no change when the drainage length of each fracture is different. However, in Figure 6b, the flux distribution of each fracture is different after considering the block length heterogeneity. Non-uniform block length only affects the final flux proportion of each fracture. It has no effect on the arriving time to the crossflow and last pseudosteady flow period. In other words, there is no obvious effect on the wellbore pressure drop when the partial section area is damaged or stimulated in the same drainage area, which only affect the flux proportion of each fracture.

Different Elastic Storage Ratio
To consider the different elastic storage ratio effect of each block in drainage area, three cases are considered in Table 3, and corresponding pressure and flux curves are shown in Figure 7. From Figure 6a, we can see that if the whole drainage length is constant, the dimensionless pressure drop and pressure derivative curves have no change when the drainage length of each fracture is different. However, in Figure 6b, the flux distribution of each fracture is different after considering the block length heterogeneity. Non-uniform block length only affects the final flux proportion of each fracture. It has no effect on the arriving time to the crossflow and last pseudo-steady flow period. In other words, there is no obvious effect on the wellbore pressure drop when the partial section area is damaged or stimulated in the same drainage area, which only affect the flux proportion of each fracture.

Different Elastic Storage Ratio
To consider the different elastic storage ratio effect of each block in drainage area, three cases are considered in Table 3, and corresponding pressure and flux curves are shown in Figure 7.
fracture is different. However, in Figure 6b, the flux distribution of each fracture is different after considering the block length heterogeneity. Non-uniform block length only affects the final flux proportion of each fracture. It has no effect on the arriving time to the crossflow and last pseudosteady flow period. In other words, there is no obvious effect on the wellbore pressure drop when the partial section area is damaged or stimulated in the same drainage area, which only affect the flux proportion of each fracture.

Different Elastic Storage Ratio
To consider the different elastic storage ratio effect of each block in drainage area, three cases are considered in Table 3, and corresponding pressure and flux curves are shown in Figure 7.  Figure 7a illustrates that the elastic storage ratio between different drainage blocks mainly affects the first linear flow, first radial flow and crossflow regimes. When the elastic storage ratio increases from case 1 to case 2 and case 3, which means the storage capacity of natural fracture system is greater and more fluid is stored in the natural fractures, the dimensionless pressure derivative curves show smaller groove because of less fluid transfer between natural fractures and matrix. From Figure 7b, it shows the same performance as pressure response. The flux distribution is only affected in the early stages, while the flux of each fracture converges to the same qtotal/3 at the later pseudo-steady flow  Figure 7a illustrates that the elastic storage ratio between different drainage blocks mainly affects the first linear flow, first radial flow and crossflow regimes. When the elastic storage ratio increases from case 1 to case 2 and case 3, which means the storage capacity of natural fracture system is greater and more fluid is stored in the natural fractures, the dimensionless pressure derivative curves show smaller groove because of less fluid transfer between natural fractures and matrix. From Figure 7b, it shows the same performance as pressure response. The flux distribution is only affected in the early stages, while the flux of each fracture converges to the same q total /3 at the later pseudo-steady flow period. Similarly, the less natural fracture storage capacity ω, the less corresponding flux of hydraulic fracture.

Different Crossflow Ratio
Another important parameter of dual porosity for naturally fracture system is crossflow ratio. Three cases are investigated in Table 3, and the pressure responses and flux distribution are shown in Figure 8. period. Similarly, the less natural fracture storage capacity , the less corresponding flux of hydraulic fracture.

Different Crossflow Ratio
Another important parameter of dual porosity for naturally fracture system is crossflow ratio. Three cases are investigated in Table 3, and the pressure responses and flux distribution are shown in Figure 8. The crossflow ratio represents the fluid flowing ability from the matrix to the natural fracture system. As shown in Figure 8, the different crossflow ratio only affects crossflow stage in the pressure responses and flux distribution. The larger value of the block is, the earlier the crossflow appears between different cases and the higher the flux is during the crossflow period.  The crossflow ratio represents the fluid flowing ability from the matrix to the natural fracture system. As shown in Figure 8, the different crossflow ratio only affects crossflow stage in the pressure responses and flux distribution. The larger λ value of the block is, the earlier the crossflow appears between different cases and the higher the flux is during the crossflow period.

Non-Uniform Half-Length Distribution of Fractures
With the constant total fracture length, we consider three cases of non-uniform fracture half-length distribution. The ratio of three hydraulic fracture half-lengths is shown in Table 3. Figure 9 gives the pressure response and flux distribution. The crossflow ratio represents the fluid flowing ability from the matrix to the natural fracture system. As shown in Figure 8, the different crossflow ratio only affects crossflow stage in the pressure responses and flux distribution. The larger value of the block is, the earlier the crossflow appears between different cases and the higher the flux is during the crossflow period.

Non-Uniform Half-Length Distribution of Fractures
With the constant total fracture length, we consider three cases of non-uniform fracture halflength distribution. The ratio of three hydraulic fracture half-lengths is shown in Table 3. Figure 9 gives the pressure response and flux distribution. From Figure 9a, we can see that if the total fracture half-length is constant, pressure response curves have nothing to do with different half-length distribution. However, In Figure 9b, the flux From Figure 9a, we can see that if the total fracture half-length is constant, pressure response curves have nothing to do with different half-length distribution. However, In Figure 9b, the flux distribution of each fracture is changed by half-length heterogeneity. The flux difference between each fracture becomes smaller as the flow time increases, and becomes stable gradually after crossflow stages.

Conclusions
This work investigates the pressure transient solutions of MFHWs in heterogeneous systems where the sections of the naturally fractured reservoir display different characteristics. A semi-analytical method is presented which combines source function and boundary element method. Major conclusions are summarized below: (1) This model only discretizes hydraulic fractures and interface boundaries. With a small number of segments interface discretization, it can match verified model very well on the log-log plots. As a result, this model can greatly reduce fine gridding required by numerical approach. (2) Compared to the trilinear flow model, this model can consider the radial flow around the hydraulic fracture, and the interference effect between hydraulic fractures with different formation properties. (3) As the permeability increases along the horizontal well direction, the pressure drop decreases and flow regimes become earlier. The fractures flux fluctuates with flow regimes, where the difference between high-and low-permeability blocks decreases firstly, then increases gradually, and tends to be stable finally. When the entire drainage length keeps constant, the pressure responses have nothing to do with non-uniform heterogeneous block length while flux distribution has been changed obviously. (4) The elastic storage ratio mainly affects the early flow regimes including first linear flow, first radial flow and crossflow periods. The bigger the elastic storage ratio, the smaller the groove and