A Decomposed Fracture Network Model for Characterizing Well Performance of Fracture Networks on the Basis of an Approximated Flow Equation

The gridless analytical and semianalytical methodologies can provide credible solutions for describing the well performance of the fracture networks in a homogeneous reservoir. Reservoir heterogeneity, however, is common in unconventional reservoirs, and the productivity can vary signi ﬁ cantly along the horizontal wells drilled for producing such reservoirs. It is oversimpli ﬁ ed to treat the entire reservoir matrix as homogeneous if there are regions with extremely nonuniform properties in the reservoir. However, the existing analytical and semianalytical methods can only model simple cases involving matrix heterogeneity, such as composite, layered, or compartmentalized reservoirs. A semianalytical methodology, which can model fracture networks in heterogeneous reservoirs, is still absent; in this study, we propose a decomposed fracture network model to ﬁ ll this gap. We discretize a fractured reservoir into matrix blocks that are bounded by the fractures and/or the reservoir boundary and upscale the local properties to these blocks; therefore, a heterogeneous reservoir can be represented with these blocks that have nonuniform properties. To obtain a general ﬂ ow equation to characterize the transient ﬂ ow in the blocks that may exhibit di ﬀ erent geometries, we approximate the contours of pressure with the contours of the depth of investigation (DOI) in each block. Additionally, the borders of each matrix block represent the fractures in the reservoir; thus, we can characterize the con ﬁ gurations of complex fracture networks by assembling all the borders of the matrix blocks. This proposed model is validated against a commercial software (Eclipse) on a multistage hydraulic fracture model and a fracture network model; both a homogeneous case and a heterogeneous case are examined in each of these two models. For the heterogeneous case, we assign di ﬀ erent permeabilities to the matrix blocks in an attempt to characterize the reservoir heterogeneity. The calculation results demonstrate that our new model can accurately simulate the well performance even when there is a high degree of permeability heterogeneity in the reservoir. Besides, if there are high-permeability regions existing in the fractured reservoir, a BDF may be observed in the early production period, and formation linear ﬂ ow may be indistinguishable in the early production period because of the in ﬂ uence of reservoir heterogeneity.


Introduction
Hydraulic fracturing technology has significantly stimulated the production of unconventional resources, such as tight and shale oil/gas, for the last two decades over the world [1][2][3][4][5]. A complex fracture network often ensues after the fracturing treatment because of the stress anisotropy and preexisting natural fractures in the unconventional formations [6,7]. Evaluating the well performance of the fractured reservoirs becomes a fundamental task for the petroleum industries to optimize the production and maximize the profit, but the complex fracture networks accompanied with reservoir heterogeneity pose a severe challenge for the engineers to conduct the reservoir simulations for these reservoirs.
Production of such reservoirs with fracture networks and spatial heterogeneity is often evaluated with numerical simulation methods. The numerical methods discretize the reservoir into grid cells to capture the spatial heterogeneities and the configurations of fracture networks. But the intrinsic deficiencies make these methods difficult to be applied in complex fracture network reservoirs. The use of the structured grid system in conventional simulators is normally constrained to the modeling of planar fractures [8][9][10] and orthogonal fracture networks [11][12][13][14]. Although the unstructured grid system provides a feasible way for simulating both the complex fracture networks and the reservoir heterogeneity [15][16][17], the heavy load of building such a model and its low computation efficiency render this method less attractive.
Analytical and semianalytical approaches, such as the multiple-porosity model, Green function method, and stimulated reservoir volume (SRV) model, have been extensively used to model the production of fracture networks. In the multiple-porosity model, the fractures are upscaled into the multiple-porosity model to take hydraulic fractures and natural fractures into consideration [11,[18][19][20]. But this model cannot reflect the real configurations of the fracture networks. The Green function method simulates the complex fracture network by discretizing the fractures into small panels and then coupling the contribution of these panels [7,21]. This method can treat the fractures explicitly as well as simulate the well performance efficiently. Besides, we can also summarize the effect of fractures as an enhanced permeability within the SRV, which, however, limits our insight into the transient flow between the fractures and the reservoir matrix [22,23]. With their efficiency and stability, these analytical and semianalytical methods enhance our understanding of the transient behavior of fractured wells and help the industry optimize the well operations. However, these approaches bear stringent restrictions because they require the reservoir matrix to be homogeneous in most cases [24,25]. The unconventional reservoirs, however, often exhibit a high degree of heterogeneity, and their productivity can vary significantly along the horizontal wells drilled for producing such reservoirs ( [26]; Miller et al., 2011; Cipolla and Wallace, 2014).
Many authors carried out different attempts to characterize the reservoir heterogeneity. Loucks and Guerrero, [27] obtained an analytical solution to describe the pressure drop in a circular composite reservoir. Based on the reflectiontransmission concept of electromagnetics, Kuchuk and Habashy [28] presented a new general method for solving the pressure diffusion equation in laterally composite reservoirs. Liu and Wang [29] studied the transient behavior of a two-dimensional flow in layered reservoirs. The reservoir heterogeneity was characterized by assigning different permeabilities to different layers. Medeiros et al. [18] divided the reservoir into blocks to model the pressure transient behavior in compartmentalized reservoirs. However, these methods can only model simple cases involving reservoir heterogeneity. They cannot be applied to deal with a heterogeneous reservoir with a fracture network.
Based on the aforementioned arguments, one can conclude that an analytical method which can simultaneously model the fracture network and the reservoir heterogeneity is still lacking. In the real field cases, however, the fracture network, accompanied with reservoir heterogeneity, can be frequently present in unconventional reservoirs. Hence, it is imperative for us to develop analytical/semianalytical methods, which are convenient and stable for computational use, to model the fluid flow in heterogeneous fracture network reservoirs.
From the fracture network model provided in Nagel et al. [30] and microseismic information provided in Warpinski et al. [14], one can see that the fractures and reservoir boundary may divide the reservoir into a set of matrix blocks. Following these real-life scenarios, we propose a new approach to simulate the flow behaviors of fracture networks in heterogeneous reservoirs. Figure 1(a) shows a simple fracture network. We first discretize the fractured reservoir into matrix blocks as shown in Figure 1(b) and model the flow behavior within each block, and then, we couple all of the flow behavior in the matrix blocks to predict the pressure response within the entire reservoir. As the flow behavior in each matrix block is individually characterized in our method, these blocks can have different properties, enabling us to take into account the reservoir heterogeneity. Additionally, the borders of the blocks can reflect the configuration of the fractures, as shown in Figure 1(c); thus, a complex fracture network can be represented by assembling all of the borders of the matrix blocks.

Methodology
When being applied for reservoir simulation purposes, the semianalytical methods have several distinct advantages, such as high computation efficiency and high stability. To obtain a semianalytical solution to predict the well 2 Lithosphere performance in heterogeneous reservoirs, we make the following fundamental assumptions: (i) The reservoir model is assumed to be isotropic and sealed by an upper and lower impermeable layer (ii) The reservoir model is simplified into a twodimensional model, and the temperature variation in the reservoir is neglected (iii) All of the fractures penetrate the formation completely in the vertical direction, and the fluid flows into the wellbore only through fractures (iv) Only single-phase oil flow is considered in this work, and the well is constrained with a constant production rate It is worth noting that, although only oil flow is considered in this work, this proposed method can be readily extended to characterize gas flow or gas condensate flow by use of the concepts of pseudopressure and pseudotime. The definitions of pseudopressure and pseudotime can be found in Singh and Whitson [31].

Fluid Flow in the Matrix.
In real-life scenarios, the matrix blocks that are divided by the fractures and/or reservoir boundary may have different geometries, and it is difficult to obtain a general flow equation to characterize the fluid flow in these blocks. In this work, we propose a new method to overcome this difficulty. The core idea of our method is that we assume that the contours of pressure can be approximated with the contours of the depth of investigation (DOI) within each matrix block. It is noted that, although this assumption cannot capture the physical scenarios, based on this assumption, we derive a new flow equation that can be applied in the matrix block with an arbitrary geometry. This idea is validated with case studies in the validation section, and the proposed flow equation can be expressed as where The derivation details for Equation (1) can be referred to Appendix A.
It is worth mentioning that ωðrÞ is equal to 0, 1/r, and 2/r for linear flow, radial flow, and spherical flow, respectively, such that Equation (1) can be readily transformed into the well-known linear flow equation, radial flow equation, and spherical flow equation. For a hydraulically fractured reservoir, where the reservoir can be discretized and represented with N m matrix blocks, we can have the following equation to characterize the fluid flow in different matrix blocks: For a matrix block with a given geometry, the ω function in Equation (4) can be readily obtained, and we can analytically characterize the fluid flow in this matrix block by solving Equation (4). Subsequently, the fluid flow in the global reservoir matrix can be described by coupling the flows from all the matrix blocks. In this work, we provide the analytical pressure functions for characterizing the fluid flow in some typical geometries, and these analytical functions are tabulated in Tables 1, 2, and 3. Appendix B shows the details of how these pressure functions are obtained. For convenience, the pressure within the n th matrix block is implicitly expressed as p n = f n k n , μ n , α n , ω n , R n , h, q n , r, t ð Þ : ð5Þ

Pressure Equation of
Fractures. In the unconventional reservoirs, where the matrix permeability can be extremely low, the pressure drop along the fractures can be neglected compared with that in the matrix; therefore, we characterize the fracture pressure with an average pressure p f . This approximation is especially reliable when the matrix permeability is in a range of microdarcies or below. Based on our assumption that the fluid flows from the matrix to the fractures and then to the wellbore and there is no direct fluid transportation between the matrix and the wellbore, the material balance equation for the fracture system can be expressed as In the fractures, the convergence flow near the wellbore can be simplified with a radial flow (see Figure 2), and the 3 Lithosphere Table 2: General pressure functions for block geometries that have an inscribed circle.
a /2 4 Lithosphere relationship between the bottom-hole pressure (p w ) and the fracture pressure (p f ) can be expressed as In summary, the semianalytical model for describing the fluid flow in a fracture network reservoir includes the following elements: (i) The analytical solutions of the matrix block flow which can be obtained by solving Equation (4) (ii) The material balance equation of the fracture system given by Equation (6) (iii) The convergence flow equation near the wellbore given by Equation (7) 2.3. Solution Methodology. We discretize the entire production period into N t timesteps. Note that the solution of Equation (4) is based on a constant production rate condition, but the contribution of each matrix block to the total production rate is time-dependent in nature. So, we apply the material balance time to make the constant production rate solution applicable to a transient production rate condition. For the matrix blocks, r = 0 indicates the location of the fracture system; thus, we have the following at r = 0: ð9Þ Table 3: General pressure functions for rectangle blocks.

Lithosphere
Equations (6) and (7) at the j th timestep can be rewritten, respectively, as The governing equation at the j th timestep can be expressed as The unknowns at the j th timestep include the following: The governing equations at the j th timestep include the following: (i) N m matrix flow equations at matrix blocks: Equation (8), which can be explicitly and analytically obtained by solving Equation (4) for each matrix block (ii) N m material balance time equations at matrix blocks: Equation (9) (iii) Material balance equation for the fracture system: Equation (10) (iv) Convergence flow equation near the wellbore: Equation (11) This group of nonlinear equations can be readily solved with the Newton method at each timestep. Furthermore, although the proposed methodology is dedicated to modeling single-phase oil flow, a minor modification on the basis of the pseudopressure concept can be made to enable the proposed methodology to be capable of modeling single-phase gas flow.

Validation of the New Semianalytical Model
We first apply the new semianalytical method to a multistage hydraulic fracture model and a fracture network model, respectively, and then validate the calculated results against Eclipse. Each model considers two scenarios: (1) the reservoir matrix is homogeneous, and (2) the reservoir matrix is heterogeneous. Figure 3 shows the configurations of the multistage hydraulic fracture model, while Figure 4 shows the configurations of the fracture network model. To take heterogeneity into consideration, the shadowed matrix blocks shown in Figures 3(b) and 4(b) are set with a higher permeability than the transparent blocks. As seen in Figure 3, the distance between two neighboring fractures is 200 m, the fracture's half-length is 200 m, and the reservoir thickness is 20 m. As seen in Figure 4, the dimension of the fracture network model is 900 × 700 × 20 m, and the side length of each square matrix block is 100 m. As for the four models mentioned above, the wells keep producing for 2000 days at a constant rate of 15 m 3 /d. All these models have been also built with Eclipse; the grid dimension used is 10 × 10 × 5 m, and the local grid refinement technique is applied to characterize the fractures. Table 4 shows the other reservoir/fracture parameters used in these four models.  Table 1. Adopting the solution methodology introduced in Section 2, we predict the well performance for both the homogeneous case and the heterogeneous case. Figure 5 compares the pressure drops and pressure derivatives calculated from the proposed method and those calculated by Eclipse for the homogeneous multistage hydraulic fracture model. The pressure drop plot calculated by our method is in excellent agreement with that obtained by Eclipse. Similar observation can be also found for the pressure derivative plot. In the pressure derivative plot, the     Figure 6 compares the pressure drops and pressure derivatives calculated from the proposed method and those calculated by Eclipse for the heterogeneous multistage hydraulic fracture model. Different from the pressure derivative plot for the homogeneous model (see Figure 5), the pressure derivative plot for this heterogeneous model exhibits two unit slope periods. The high-permeability blocks exhibit the earlier BDF, corresponding to the first unit slope in the pressure derivative plot. The flow behavior in the highpermeability blocks leads to the appearance of the first unit slope period because the pressure response in the highpermeability blocks tends to be more rapid than that in the low-permeability blocks. The transition period between these two unit slopes represents the superposition of two flow regimes, i.e., the BDF regime in the high-permeability blocks and the linear flow regime in the low-permeability blocks. When the flow regimes in all the blocks enter the BDF period, the second unit slope appears in the pressure derivative plot.

Fracture Network Model.
In this fracture network model, the fractures divide the reservoir into two kinds of matrix blocks, namely, the square matrix blocks sealed by fractures and the outer SRV matrix block sealed by the reservoir boundary and fractures (see Figure 4). Table 2 provides the pressure function that characterizes the fluid flow in the square matrix blocks, while Appendix C shows the pressure function that characterizes the fluid flow in the outer SRV matrix block. Appendix C also gives the details for deriving the pressure function. Figure 7 compares the pressure drops and pressure derivatives calculated from the new semianalytical method and those calculated by Eclipse for the homogeneous fracture network model, while Figure 8 provides the same comparison for the heterogeneous fracture network model. Again, we can achieve an excellent match between the pressure drops and pressure derivatives calculated by our method and those by Eclipse. It can be seen from Figure 7 that, in the homogeneous case, there is no distinguishable half-unit slope or unit slope appearing in the pressure derivative plot due to the compound transient flows of the square blocks and the outer SRV block. In comparison, it can be observed from Figure 8 that a unit slope in the early production period appears in the pressure derivative plot, which indicates that BDF      Figure 9 shows a flowchart detailing the procedures used for applying the new semianalytical approach to heterogeneous fracture network reservoirs. The essential steps are as follows:

Application of the New Semianalytical Model
(i) Discretize the fractured reservoir into matrix blocks and upscale the local properties to the blocks (ii) Predict the contours of DOI within each block and obtain the ω function for each block (iii) Substitute the ω functions into Equation (4) and obtain the pressure function for each block (iv) Divide the entire production duration into discrete timesteps. At each timestep, construct and solve the governing equation (i.e., Equation (12)) that describes the fluid flow in the heterogeneous fracture network reservoir In this section, we apply this new method to a heterogeneous reservoir that is, respectively, equipped with an orthogonal fracture network, a nonorthogonal fracture network, and a compound fracture network. Figure 10 Table 5 lists the fluid/rock properties that are incorporated into the heterogeneous model considered in this section. To illustrate the influence of heterogeneity on the pressure response, we compare the pressure drops and pressure derivatives calculated for the heterogeneous reservoir against those calculated for a counterpart homogeneous reservoir with a uniform permeability of 0.002 mD. Using our newly developed semianalytical model, we perform calculations over a production period of 2000 days. Figures 11-16 show the calculation results obtained for the three aforementioned models.

Fracture Network #1: Orthogonal Fracture Network
Model. This orthogonal fracture network model divides the reservoir into a set of square matrix blocks characterized with different permeabilities. Block geometries (2) and (6) provided in Table 2 and block geometry (3) provided in Table 3 can be used to characterize this fracture network model. By using these block geometries to represent the entire fracture network reservoir, we can obtain the simplified reservoir model equipped with fracture network #1 as shown in Figure 11(a). Blocks with different permeabilities can be observed in this fracture network model. Figures 11(b)-11(d) show the pressure drop distributions on the 10 th , 100 th , and 1000 th production days, respectively, which are calculated by our semianalytical model. From Figure 11, one can observe that the pressure drops much faster in the high-permeability blocks than in the low-permeability blocks. Figure 12 presents the pressure drop and pressure derivative plots for the heterogeneous reservoir and the homogeneous reservoir. On the pressure derivative plots, a half-unit slope can be readily recognized in the homogeneous case, which indicates a linear flow from matrix blocks to the fracture system, whereas in the heterogeneous case, there is no distinguishable flow regime in the early production period; a reasonable explanation is that the coupling of the BDF appearing in the high-permeability blocks and the linear flow appearing in the low-permeability blocks render the flow regime indistinguishable in the early production period. In the late production period, the flow regime in all of the blocks enters BDF, and hence, the pressure derivative plot shows a unit slope in both the homogeneous case and the heterogeneous case.  9 Lithosphere orthogonal, dividing the reservoir into irregular matrix blocks. Fracture network #2 is used to simulate such a case. We can use block geometries (1) and (6) provided in Table 2 and block geometry (3) provided in Table 3 to characterize fracture network #2. Figure 13(a) presents the simplified fracture network model that is built with these simple block geometries. As one can see in Figure 13(a), the fractures divide the reservoir into triangle blocks and unsealed square blocks. Figures 13(b)-13(d) show the pressure drop distributions on the 10 th , 100 th , and 1000 th production days that are calculated with our semianalytical model. A more rapid pressure drop can also be observed in the highpermeability blocks than in the low-permeability blocks. Figure 14 compares the pressure drops and pressure derivatives obtained for the heterogeneous reservoir with those obtained for the homogeneous reservoir; both reservoirs are equipped with fracture network #2. The pressure derivative plots shown in Figure 14 are similar to those shown in Figure 12 which is dedicated to fracture network #1. In the early production period, a linear flow regime can be recognized on the pressure derivative plot in the homogeneous case, while the flow regime is indistinguishable in the heterogeneous case because of the influence of heterogeneity. A unit slope, which indicates BDF, can be observed in both the homogeneous case and the heterogeneous case in the late production period.

Fracture Network #3: Compound Fracture Network
Model. Fracture network #3 simulates such a case that the fracture network divides the reservoir into matrix blocks with different geometries as well as different permeabilities. Figure 15(a) shows the simplified reservoir model after applying the block geometries provided in this work to characterize fracture network #3. As seen from Figure 15(a), the reservoir is divided into triangle, rectangle, and square blocks with different dimensions. We characterize this fracture network with block geometries (1), (2), (5), and (6) listed in Table 2 and block geometries (1) and (3) listed in Table 3. Figures 15(b)-15(d) show the pressure drop distributions  on different production days. Again, the pressure drops much faster in the high-permeability blocks than in the low-permeability blocks. Figure 16 compares the pressure drops and pressure derivatives for both the heterogeneous reservoir and the homogeneous reservoir that are equipped with fracture network #3. The pressure drop plots and pressure derivative plots shown in Figure 16 exhibit similar trends as observed in Figures 12 and 14 that are dedicated to fracture network #1 and fracture network #2, respectively. A linear flow regime can be identified in the homogeneous case, while it cannot be observed in the heterogeneous case.
The BDF leads to a unit slope in the late production period in both cases.

Conclusions
In this work, we propose a novel semianalytical approach to simulate the well performance in heterogeneous fracture network reservoirs. This approach is dedicated to the cases that the hydraulically fractured reservoir can be discretized into matrix blocks which are sealed with the fractures and/or reservoir boundary. By coupling the fluid flow in these matrix blocks, we develop a new semianalytical methodology that can be used to model a fracture network in a heterogeneous reservoir. The key features of this new method can be summarized as follows: (i) By approximating the contours of pressure with the contours of DOI, we derive a novel general flow  11 Lithosphere equation that can be applied to characterize the fluid flow in different block geometries. As such, the well performance in a complex fracture network reser-voir can be simulated by discretizing the reservoir into a series of blocks (ii) Compared with the semianalytical methodologies that discretize the fracture network into fracture panels to capture the fracture network configuration, our method discretizes the fractured reservoir into matrix blocks (iii) As the fluid flow in each matrix block is characterized individually, the properties of different blocks do not have to be uniform, such that we can take reservoir heterogeneity into consideration by upscaling the local properties to the matrix blocks We apply our method to several synthetic reservoirs and validate our method against the commercial simulator Eclipse. The simulated results provide us with the flow regime information for these complex fracture networks. The findings include the following: (i) If there are high-permeability blocks existing in the fractured reservoir, a BDF may be observed in the early production period (ii) Linear flow may be indistinguishable in the early production period because of the influence of reservoir heterogeneity  12 Lithosphere In the unconventional reservoirs, the pressure "front" can reach the furthest tip of the fractures in a few seconds. Compared with the time of a few months or even years it may take for the pressure response to cover the entire matrix block, the pressure response within the fractures can be regarded as instant. Therefore, we divide the propagation of DOI into two periods: first, the DOI expands in the fracture system; second, the DOI travels from the borders of the matrix blocks to the interior of the blocks, which is illustrated in Figure 17. In Figure 17, r is the DOI, t is the time, and R is the minimum DOI that can cover the entire matrix block. Figure 17(a) shows a matrix block in a fracture network. In Figures 17(b)-17(d), the fracture in the two-dimensional plane is considered to be composed of a set of point sources, and the DOIs simultaneously expand outwards from the center of these sources. Figure 17(b) shows that the DOIs just start to expand and the matrix block has not been investigated at all. Figure 17(c) shows some point sources along the fractures; at this time, the matrix blocks have been partially investigated. Figure 17(d) illustrates that the matrix block has been totally investigated when the DOIs reach a certain  It should be noted that, based on the assumption that the contours of pressure can be approximated with the contours of DOI, the parameter r in Equation (A.1) indicates the DOI. The inner boundary (r = 0) represents the block borders (fractures), while the outer boundary (r = R) represents the minimum DOI that can cover the entire matrix block. Figure 18 shows a matrix element (grey region) that is sealed by the contours of DOI. The single-phase mass balance equation for this matrix element can be written as The porosity of the matrix will change during the production as per

B. Derivations of the Analytical Solutions
On the basis of our assumption that the contours of pressure can be approximated with the contours of DOI within each matrix block, we can characterize the fluid flow within each matrix block if we can predict the contours of DOI in that block. In a fractured reservoir, the fractures divide the reservoir into matrix blocks with different geometries. We provide the pressure functions for some typical block geometries, and the derivations are detailed below.
B.1. Block Geometries That Have an Inscribed Circle. In a homogeneous matrix block, the propagation velocity of DOI along the borders is approximated to be identical in this work, enabling us to predict the contours of DOI in some specific block geometries. Figure 19 shows the approximated contours of DOI in a triangle block, a square block, and a quadrilateral block, respectively. In Figure 19, the dark dash lines are the approximated contours of DOI, the solid circles are the inscribed circles for these blocks, and R in are the radii of the inscribed circles. For these block geometries, the investigated region can cover the entire block if the DOI equals R in . According to the principle of similarity, where C 0 is the circumference of the matrix block. 16 Lithosphere Figure 23 shows the approximated contours of DOI in the outer SRV matrix block. We approximate the reservoir boundary with the outmost contour that has the same area as that of the reservoir. The real reservoir area equals the area within the outmost contour, such that where a 1 is the length of the reservoir, b 1 is the width of the reservoir, a 2 is the length of SRV, b 2 is the width of SRV, and R o is the DOI of the outmost contour. R o can be obtained from Equation (C.1): : ðC:2Þ The circumference of the contours in the outer SRV matrix block can be expressed as  Figure 21: The matrix blocks that can be converted into a rectangle block.