Semianalytical Model for Flow Behavior Analysis of Unconventional Reservoirs with Complex Fracture Distribution

The secondary fracture spacing/aperture has a great influence on flow behavior of multifractured horizontal wells. In this paper, the effects of secondary fracture spacing/aperture distribution on porosity/permeability of the stimulated reservoir volume are described. A new fractal composite flow model was proposed, which can comprehensively characterize the secondary fracture spacing/aperture heterogeneity in the stimulated reservoir volume. Results showed that the larger the fractal dimension of the secondary fracture aperture is, the greater the dimensionless flow rate will be. The smaller the fractal dimension of the secondary fracture spacing is, the smaller the dimensionless flow rate will be. The secondary fracture spacing distribution had a greater influence on flow behavior than the secondary fracture aperture distribution. If the heterogeneity of the secondary fracture spacing or aperture is not considered, the bottomhole pressure and production curve will change greatly. At last, the proposed model was applied to match the production data of two actual fields.


Introduction
Unconventional reservoirs, such as tight or shale reservoirs, usually have very low permeability, and it is difficult to obtain economic production [1,2]. The hydraulic fracturing techniques have been widely used to efficiently enhance production of unconventional reservoirs [3]. As of 2012, approximately 85% of shale gas wells in the United States were produced using horizontal wells and multistage fracturing techniques [4]. There are a large number of complex secondary fractures created in the area around hydraulic fractures [5], and the area is called the stimulated reservoir volume (SRV) [6]. In the past two years, domestic scholars have studied the factors affecting the SRV; the results showed that the matrix permeability and the secondary fracture spacing play a key role in the SRV [5]. It can be seen from the 20-year EUR analysis that the contribution of reservoir production in the unstimulated area is negligible [3,7,8].
The distribution of natural fractures in unconventional reservoirs is complex, and reservoir heterogeneity and stress fields make secondary fracture spacing and aperture distribution more complex [9]. Therefore, in order to simulate the secondary fracture distribution, scholars proposed a nontraditional fracture model [10,11], which describes the hydraulic fracture and secondary fracture distribution in detail through the unstructured grid. In contrast, the network model uses two orthogonal planar grids to represent the fracture network, which is a simple and effective fracture modeling method [12,13]. However, the network model idealized the distribution of the secondary fractures and could not represent the actual secondary fracture distribution. The secondary fracture distribution in SRV can be obtained intuitively through microseismic monitoring. The researches have shown that the secondary fractures around hydraulic fractures can be regressed with fractal theory models [14][15][16][17]. A lot of research showed that the fracture spacing distribution accords with the fractal scale distribution, and the fractal dimension was defined to characterize the fracture spacing distribution [18][19][20][21]. The fractal dimension obtained from natural fractures was applied directly to the secondary fracture analysis of fractured horizontal wells, which was widely used to analyze flow behavior influenced by complex fractures in the analytical/semianalytical models [22][23][24][25][26].
In most studies, the secondary fracture aperture is generally assumed to be uniformly distributed, or the secondary fracture aperture is considered to follow a special distribution, such as the lognormal distribution [27]. Kim and Schechter [28] used CT scans to study the secondary fracture, and the results showed that the secondary fracture aperture distribution accords with the lognormal distribution. For the first time, Sheng et al. [21] proposed the fractal dimension of the fracture aperture to characterize the secondary fracture aperture distribution. The results show that the fractal equation is more general then specific distributions, such as exponential, linear, power, and uniform distributions.
In the studies of flow behavior of fracturing horizontal wells, the numerical simulation, analytical, and semianalytical models are widely used. Numerical simulation models are more accurate, but they are computationally consuming and not suitable for actual reservoirs [12]. The analytical/semianalytical method is simple and fast, and the linear flow models, including bilinear flow models [29], trilinear flow models [30], and the five-zone composite linear flow model [8], are commonly used in multifractured horizontal wells. The linear flow model proposed by the above scholars generally describes the SRV region using the dual-media model, while the USRV region is mainly described by the singlemedia model [30]. The trilinear flow model is mainly applicable to the case where the hydraulic fracture spacing is small and the reservoir between the hydraulic fractures is completely stimulated. The five-zone composite linear flow model is suitable for the case where the reservoir between the hydraulic fractures is not completely stimulated. Therefore, the five-zone composite linear flow model is more general. The secondary fracture porosity/permeability and fracture spacing/aperture distribution have a great influence on the pressure response. At present, CDE (uniform fracture spacing and uniform porosity/permeability, [31][32][33][34][35]), MCDE (variable fracture spacing and uniform porosity/permeabil-ity, [36][37][38][39]), FDE (uniform fracture spacing and fractal porosity/permeability, [22][23][24][25][26]40]), and DMFDE (variable fracture spacing and variable porosity/permeability, [41]) are mainly used to simulate flow behavior in the dualmedia model. Although the influence of the secondary fracture spacing/aperture on the porosity/permeability of the fracture system has been studied, respectively, there is currently no analytical/semianalytical model that can simultaneously consider the influence of the secondary fracture spacing and aperture on fluid flow.
In this paper, based on the five-zone composite linear flow model, a fractal composite flow model of multifractured horizontal wells that can comprehensively characterize the secondary fracture spacing/aperture distribution in the SRV region was proposed. The effects of secondary fracture spacing/aperture on fracture porosity/permeability in the SRV region are described by fractal theory. Based on the model, the influence of secondary fracture heterogeneity on the flow behaviors was analyzed, and the model was applied to two actual reservoirs.

The Fractal Composite Flow Model of
Multifractured Horizontal Wells 2.1. Physical Model. It is assumed that the horizontal well is located in the center of the reservoir, and on this basis, it is further assumed that the fluid flow is completely symmetrical about the wellbore. The flow pattern of the single-fracture control region is symmetrically distributed with respect to the hydraulic fracture, so that only a quarter of one fracture control flow area needs to be considered in the research process ( Figure 1). The region I is the hydraulic fracture of a multifractured horizontal well. The hydraulic fracture has a relatively small width, the fracture aspect ratio is large, and the permeability within the hydraulic fracture is relatively large. The region II is the SRV. The natural fractures in this area are further extended, and the secondary fractures are created to form a complex secondary fracture network. The heterogeneity of secondary fracture spacing/aperture distribution in this area has a great impact on fluid flow and needs to be considered. Region III, region IV, and region V are the where k m is the matrix permeability, D. p 5 is the pore pressure in region V, atm. φ m is the matrix porosity, dimensionless. c tm is the matrix compression factor, atm -1 . μ is viscosity, mPa·s. t is time, s.
where p 4 is the pore pressure in region IV, atm. x f is the halflength of the hydraulic fracture, cm.
where p 3 is the pore pressure in region III, atm. For the convenience of calculation, the dimensionless parameters are defined as follows: where p nD is the dimensionless pressure, dimensionless. k f is the fracture permeability of the dual-media model in the SRV region, D. h is the reservoir thickness, cm. p ini is the initial reservoir pressure, atm. q is flow rate, cm 3 /s.
where t D is the dimensionless time, dimensionless. φ 2 is the total porosity of the SRV region, dimensionless. c t2 is the total compression factor of the SRV region, atm -1 .
The fluid flow equation in the unstimulated area is subjected to dimensionless, converted to Laplace space, and combined with the boundary conditions of different regions; the flow equations of different regions can be obtained, which are (i) Region V: (ii) Region IV: (iii) Region III: 2.2.2. Flow Equations in Stimulated Reservoir Volume. In the SRV region, the influence of the secondary fracture spacing/aperture on fracture porosity/permeability is described by the fractal secondary fracture spacing distribution (FFSD) and fractal secondary fracture aperture distribution (FFAD) [21], as follows: where φ f is the fracture porosity of the SRV region, dimensionless. φ w is the fracture porosity at the reference point, dimensionless. x w is the position of the reference point, cm. d f s is the fractal dimension of the secondary fracture spacing, dimensionless.
where k w is the fracture permeability at the reference point, D. d f a is the fractal dimension of the secondary fracture aperture, dimensionless. θ is the fractal index, dimensionless. Another important parameter in the fractal porosity/permeability equation, the fractal index, was also fully discussed in the study of Sheng et al. [21]. Unlike conventional fractal porosity/permeability models, Sheng et al. [21] believe that 3 Geofluids the fractal index is nonuniformly distributed in the actual reservoir. This is closer to the actual situation. It is considered that the fractal index can be calculated by the d f s [21] where s f w is the secondary fracture spacing at the reference point, cm.
In this paper, we combine the fractal fracture porosity/permeability models (Equations (9), (10), (11)) with the flow equations for the SRV region in a conventional five-zone composite linear flow model. Therefore, it can accurately describe the fluid flow in a complex secondary fracture of the SRV region. Considering that there is pseudo-steadystate flow between the matrix and the fracture, the flow equations that take into account the secondary fracture spacing/aperture heterogeneity can be expressed as where α is the shape factor of matrix blocks in the dual-media model, cm -2 . p f is the pore pressure in secondary fractures, atm. p m is the pore pressure in matrix blocks, atm. And where b f is the half-width of the hydraulic fracture, cm. q ″ is the mass transfer from secondary fractures to the matrix, cm 3 /s. ρ is the fluid density, g/cm 3 . Similarly, the fluid flow equation in the SRV is subjected to dimensionless, converted to Laplace space, and combined with the boundary conditions, which are where λ is the crossflow coefficient, dimensionless. ω is the storage ratio of the secondary fracture to matrix blocks, dimensionless.

Flow Equations in Hydraulic
Fracture. Region I is the linear flow of the hydraulic fracture. Regardless of the effect of the skin and wellbore storage, the flow equation in the hydraulic fracture is where p 1 is the pore pressure in the hydraulic fracture, atm. k 1 is the hydraulic fracture permeability, D. φ 1 is the hydraulic fracture porosity, dimensionless. c t1 is the compression factor of the hydraulic fracture, atm -1 .
Similarly, the fluid flow equation in the hydraulic fracture is subjected to dimensionless, converted to Laplace space, and combined with the boundary conditions, which are where s is the Laplace variable.

Model Solution.
In the semianalytical solution of the fractal composite linear flow model, the region V, the region III, the region IV, the region II, and the region I are solved, and finally, the pressure behavior of the well bottomhole can be obtained: where p wD is the dimensionless bottomhole pressure in the Laplace domain (dimensionless), where According to the relationship between the bottomhole pressure and the flow rate under constant pressure conditions in the Laplace domain, the well production can be calculated as where q D is the dimensionless flow rate in the Laplace domain, dimensionless. Using Stehfest numerical inversion, the production or pressure in the Laplace domain can be inverted to obtain the pressure and flow rate in real space.

Influence of Secondary Fracture Aperture on Pressure and
Production. The effects of the secondary fracture aperture distribution on fluid flow and production are analyzed. Assuming that the fractal dimensions of the secondary fracture aperture are 1.7 and 1.9, respectively, the secondary fracture aperture at different positions can be calculated by the FFAD equation, as shown in Table 1. At the same time, in order to analyze the influence of secondary fracture aperture heterogeneity when the total secondary fracture apertures are the same, we set two cases. The secondary fracture aperture in Case 1 is the average aperture when d f a is 1.7. The secondary fracture aperture in Case 2 is the average aperture when d f a is 1.9. As the four cases show in Table 1, the influence of secondary fracture aperture heterogeneity on pressure and production is analyzed, as shown also in Figure 2.
The fractal dimension of the secondary fracture aperture has little effect on the dimensionless pressure and flow rate in the linear flow period. It can be seen from the figure that the larger the d f a is, the smaller the dimensionless pressure in the middle of production period will be and the greater the dimensionless flow rate will be. The d f a has little effect on the slope of the pressure derivative curve before the pressure response reaches the outer boundary, which means that the equivalent permeability distribution of the fracture medium has little effect on the flow behavior. It can be seen from the figure that the larger the d f a is, the earlier the pressure response reaching the outer boundary will be. It should be noted that if the heterogeneity of the secondary fracture aperture is not considered, the bottomhole pressure and production curve will change greatly. This means that when the mean value is used to characterize the secondary fracture aperture, the difference between the calculated bottomhole pressures or production and the actual data is relatively large.

Influence of Secondary Fracture Spacing on Pressure and
Production. The effects of secondary fracture spacing heterogeneity on fluid flow and production are analyzed. Assuming that the d f s are 1.7 and 1.9, respectively, the secondary fracture spacing at different positions can be calculated by the FFSD equation. At the same time, in order to analyze the influence of secondary fracture spacing heterogeneity when the total secondary fracture spacing is the same, we also set two cases. The secondary fracture spacing in Case 3 is the average spacing when d f s is 1.7. The secondary fracture spacing in Case 4 is the average spacing when d f s is 1.9. As the four cases show in Table 2, the influence of secondary fracture spacing heterogeneity on pressure and production is analyzed, as shown also in Figure 3.
It can be seen from the figure that the smaller the d f s is, the larger the dimensionless pressure and the smaller the dimensionless flow rate will be. The slope of the pressure derivative curve increases with the increase of d f s in the early stage of production, close to 1. This is because when the total length of the secondary fracture spacing remains the same, the larger the d f s is, the more uniform the secondary fracture spacing and the equivalent SRV porosity/permeability distribution will be and the larger the secondary fracture spacing near the hydraulic fracture will be. Since less fluid will flow from the secondary fracture to matrix blocks at the early stage, the larger the d f s is, the more pronounced the cross flow will be. It should be noted that when considering the x rD , 5 Geofluids uniform distribution of the secondary fracture spacing (Case 3 and Case 4), the linear flow and cross flow will be more obvious. Similar to the influence of the secondary fracture aperture, if the heterogeneity of the secondary fracture spacing is not considered, the bottomhole pressure and production curve will also change greatly.

Influence of Tortuosity Fractal Index on Pressure and
Production. In the fractal theory, the fractal index is generally used to characterize the influence of fracture tortuosity. Therefore, the distribution of the tortuosity fractal index should also be heterogeneous. The effects of tortuosity fractal index heterogeneity on fluid flow and production are analyzed. Also four cases are given. Case 5 and Case 7 have the same average tortuosity fractal index, and the tortuosity fractal index of Case 5 is calculated when the d f s is 1.7, which considers the heterogeneity of the tortuosity fractal index. Similarly, Case 6 and Case 8 have the same relationship, and their fractal index is calculated when the d f s is 1.9. The distribution of the tortuosity fractal index of the four cases           Table 3. The effects of the tortuosity fractal index on pressure and production are shown in Figure 4. The tortuosity fractal index represents the length of fluid flow distance in the fracture network. The smaller the fractal index is, the shorter the fluid flow path will be. The tortuosity fractal index has a certain influence on the whole flow stage. The smaller the fractal index is, the lower the bottomhole pressure will be and the higher the production rate will be. It can be seen from the figures that the fractal index heterogeneity has a great impact on bottomhole pressure and production, especially in the middle production period.
3.4. The Application of Proposed Model in Actual Field. Here, we apply the proposed model to match the production data of two actual fields. The production data at all times and the average bottomhole flow pressure of the Xinjiang tight oil reservoir are given, and in the other field from Barnett shale, the bottom flow pressure and production data for each time are available. According to the production matching results in Figures 5 and 6, it can be seen that the model proposed in this paper can match well with the actual production data. At the same time, based on the production fitting, we can obtain the distribution of the secondary fracture spacing and aperture, which can help to evaluate the fracturing of the horizontal well. In the tight oil reservoir, the d f s is 1.6 and the d f a is 1.7. In the shale gas reservoir, the d f s is 1.7 and the d f a is 1.9. It can be seen that the multifractured horizontal in the Barnett shale creates more secondary fractures, which is better to enhance production.

Conclusions
A new fractal composite flow model of multifractured horizontal wells that can comprehensively characterize the secondary fracture spacing/aperture heterogeneity was proposed. The influence of secondary fracture heterogeneity on the flow behaviors was analyzed.
(1) The larger the d f a is, the smaller the dimensionless pressure in the middle production period will be and the greater the dimensionless flow rate will be. When the mean value is used to characterize the secondary fracture aperture, the difference between the calculated bottomhole pressure or production and the actual data is relatively large (2) The smaller the d f s is, the greater the dimensionless pressure and the smaller the dimensionless flow rate will be. If the heterogeneity of the secondary fracture spacing is not considered, the bottomhole pressure and production curve will also change greatly (3) The tortuosity fractal index has a certain influence on the whole flow stage; the smaller the fractal index is, the lower the bottomhole pressure will be and the higher the production will be. The fractal index heterogeneity has a great impact on bottomhole pressure and production, especially in the middle production period (4) The model proposed in this paper can match well with the actual production data. The distribution of the secondary fracture spacing and aperture can be obtained from data fitting, and the fracturing can be evaluated

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.  Figure 6: Application of proposed model in a shale gas reservoir from Barnett shale. 8 Geofluids