Three-Dimensional Numerical Simulation of Multiscale Fractures and Multiphase Flow in Heterogeneous Unconventional Reservoirs with Coupled Fractal Characteristics

A simulated reservoir model, based on the permeability fractal model and three-dimensional (3D) Gaussian filter, was established to account for in-layer and interlayer heterogeneity so that the result conforms to the law of geological statistics. Combined with an embedded discrete fracture method (EDFM), a multiscale fracture system was established, forming the numerical simulation method of multiphase flow in horizontal wells in heterogeneous reservoirs with complex fractures. The heterogeneity and saturation of the reservoir mixed five-point pattern of vertical and horizontal wells and the injection and production of horizontal wells were discussed. The results show that it is difficult to characterize complex reservoirs using a homogeneous permeability model. Thus, it is best to use a heterogeneous model that considers permeability differences in tight reservoirs. Formation fluids coexist in multiple phases, and water saturation has a direct effect on the production. Thus, a multiphase flow model is needed and can play a greater role in injection and production technology. The mixed five-point pattern of vertical and horizontal wells can improve productivity to a certain extent, but the dual effects of heterogeneity and fracturing will cause a decline in production by accelerating the communication of injected fluid. The reservoir is heterogeneous between wells, and there are differing effects on adjacent wells. Therefore, near-well natural microfractures are opened because of fracturing in horizontal wells, and the heterogeneity cannot be ignored, especially when multiple wells are simultaneously injected and produced.


Introduction
Unconventional reservoirs are generally heterogeneous and characterized by low porosity and permeability and natural microfractures. Horizontal well fracturing is often used to improve productivity, which makes fractures the main channels for the migration of oil and gas in such reservoirs [1][2][3]. Because of the importance of fractures, it is critical to improve the accuracy and applicability of numerical simulation of porous media with fractures [4][5][6].
In recent years, many scholars have studied and developed numerical simulation of fractured porous media, and the most widely used models are equivalent continuum models and discrete fracture model (DFM) [7]. However, considering that the equivalent continuum model is unable to characterize the differential flows in fractures, the DFM is deficient in the mesh generation of complex fractures [8][9][10]. Lee et al. [11] and Moinfar et al. [12] proposed the embedded discrete fracture method (EDFM). This method can greatly reduce the computational cost and is suitable for complex fracture modeling by embedding fracture meshes into matrix meshes and then constructing the flow exchange between a matrix and fracture system through nonadjacent links. Yan et al. [13] applied the EDFM in numerical simulation with anisotropy of matrix permeability. Shao and Di [14] developed the integrated EDFM. Based on the EDFM, Zhang et al. [15] established a 3D model of multistage hydraulic fracturing of horizontal wells in tight reservoirs, accounting for the influence of gravity and stress sensitivity. Fiallos et al. [16] used the EDFM to study the impact of inter-well interference on productivity. Zhu et al. [17] designed an EDFM based on corner grid and refined the near-well zone grid. Rao et al. [18] developed an EDFM preprocessing algorithm for arbitrary shapes. The EDFM is arguably the most promising method for numerical simulation of porous media with fractures, but it still has problems. Many scholars have reported that the effective simulation of fractures by the EDFM is suitable for complex fractures, while they commonly ignored the fracture scales. Reservoir models and their meshes are usually large in scale (a grid may represent tens of meters). The use of the EDFM can effectively represent fractures at a large scale, but considering natural microfractures in the reservoir will increase the complexity of fracture modeling. In addition, the matrix is regarded as homogeneous in most numerical models, which is quite different from the real system. Hence, it is of great importance to find a method that can reduce calculation difficulty in the numerical simulation of an unconventional reservoir while accounting for the complexity of the reservoir.
In general, there is a lack of practical geological models, so most researchers assign different values of permeability to different layers or replace them with random numbers lacking in physical meaning in the modeling of heterogeneous reservoirs [19]. In practice, reservoir heterogeneity follows the principles of geo-statistics1. Unfortunately, this is often neglected in numerical modeling for the purpose of simplification. Moreover, many researchers have established very detailed models to characterize the relationship between porosity and permeability of porous media [20][21][22][23]. However, most models use a theoretical calculation of permeability, and the corresponding assumptions are too idealistic.
Fractal geometry was put forward by Mandelbrot and Wheeler in the 1970s and developed rapidly after that [24]. The research object have three obvious characteristics: (1) it belongs to nature, (2) its morphology is irregular and unstable, and (3) it has the characteristics of self-similarity [25]. Except for the three features, fractal geometry generally studies structures in nature that cannot be described by traditional geometry. As a fractal structure, the pore structure of reservoir is very suitable for the study of fractal geometry because of its microscopic heterogeneity [26]. Many scholars have found that the pore microstructure of rocks meets the self-similarity of fractal characteristics, and fractal geometry theory has been widely used in the study of the transmission characteristics of porous media (such as permeability) [27]. Yu and Cheng [28] simplified the complex pore size distribution according to the fractal theory and proposed a fractal permeability model of bidispersed porous media, which addressed the defect of equal diameters in the original K-C equation. Chen and Yao [29] proposed a fractal permeability     2 Geofluids model which was theoretically deduced based on the fractal geometry theory with considering the irreducible water saturation. Dong et al. [30] used an improved statistical algorithm with binary image data to estimate the geometric parameters of each pore such as the perimeter and area to determine the fractal parameters of pore microstructure. The traditional relationship model between micropore structure parameters and macropermeability was modified. Therefore, it is a meaningful work to extend the fractal permeability model to numerical simulation.
The main objective of this work was to combine the heterogeneity caused by natural microfractures with a fractal permeability model and a 3D Gaussian filter of the matrix, which results in a model that conforms to statistical laws. Using an EDFM, fractures in large-scale and artificial fractures were established to form a multiscale fracture system, the solution of which is determined based on the framework of the open-source simulator black-oil model in MATLAB Reservoir Simulation Toolbox (MRST) [31]. By comparing with the model established by Zhang, the correctness of the simplified homogeneous version of the proposed model was verified. The effects of permeability and saturation on the productivity of fractured horizontal wells were discussed for homogeneous reservoirs. Then, the effects of areal heterogeneity and interlayer heterogeneity on horizontal well productivity were further discussed. In addition, the injection and production of the mixed five-point pattern of vertical and horizontal wells and adjacent horizontal wells under complex reservoir conditions were simulated to demonstrate the potential of the proposed model. The influence of heterogeneity and multiscale fractures on injection-production was thereby theoretically analyzed, which provides guidance for the reverse interpretation of production data and clarification of the main control factors for production.

Mathematical Model
where a represents different fluids. In this paper, oil phase and water phase are represented by subscripts o and w, respectively. b = 1/B, B are volume factors; ϕ is porosity (%); S stands for fluid saturation (%); μ for fluid viscosity (mPa·s); k for absolute permeability (μm 2 ); k r for relative permeability (μm 2 ); p for pressure (MPa); and q for the source term (m 3 /s). The saturation of each phase satisfies the equation S 0 + S w = 1. The compressibility coefficient of the fluid is where c is the isothermal compressibility coefficient of the fluid (1/MPa). Thus, b α in equation (1) can be expressed as where B αe is the reference volume coefficient and p e is the reference pressure (MPa).

Embedded Discrete Fracture
Model. Due to the influence of rock mechanical properties, tectonic stress, thickness of nonrock mechanical layers, and hydraulic fracturing, the fracture system in tight reservoirs usually has multiscale properties [32]. The influence of fractures of different scales on seepage system is obviously different. According to the length, width, porosity, and permeability of natural fractures, they can be divided into large fractures, medium fractures, small fractures, and microfractures [33]. Due to the difficulty of meshes division of microfractures, the physical property differences caused by them are simplified as matrix heterogeneity. Hydraulic fracturing technology can produce multiscale hydraulic fractures. Both natural and hydraulic fractures are characterized by EDFM. This means that fractures of different sizes are represented by differentiated settings based on physical properties (porosity, permeability, and degree of filling). By using EDFM, the matrix is directly divided into structural meshes; then, the fractures are embedded into the bedrock grid system, and nonneighboring connection (NNC) is formed according to the intersection of fractures and matrix, so as to realize the flow exchange between meshes that are adjacent in the physical model but not adjacent in the computational one. Because this model avoids refinement of  4 Geofluids meshes around fractures, its computing speed is greatly improved, especially when dealing with complex fracture distribution, it has more obvious advantages over the discrete fracture model [4,7,12]. Note that the crack is dimensionally reduced, so k f ∈ Rn − 1. With flow exchange between matrix and fracture system considered, the following equation can be obtained [34]: where subscript m and f , respectively, represent the matrix and fracture system; q mf is the flow exchange between matrix system and fracture system; q f f is the flow exchange among fractures; I and j are the i th and the j th fracture, respectively; T mf is the conductivity coefficient of matrix and fracture; T f f is the conductivity coefficient of fractures; and the general formula for calculating the conductivity coefficient of nonadjacent connection is [15] where K NNC is the permeability of NNC, namely, the effective permeability (μm 2 ); A NNC is the contact area of NNC pair, namely, the flow area (m 2 ); and d NNC is the relevant characteristic distance (m). A more detailed theoretical approach of EDFM can be found in the paper of Moinfar [35]. Figure 1 is a scanning electron microscope image of a block in Changqing tight reservoir of China. Natural microfractures are widely distributed in the reservoirs, and permeability of tight matrix and natural microfractures is significantly different [36,37]. EDFM is suitable for fractures in large scale, while natural microfrac-tures are in relatively small scales and have a complex distribution. In this work, permeability differences caused by natural microfractures are considered while the cost of numerical calculation is minimized to the greatest extent. Therefore, the heterogeneity of natural microfractures is equivalent to that of the matrix to establish a heterogeneous reservoir conforming to the permeability difference between natural microfractures and tight matrix. In order to characterize the porosity and permeability of heterogeneous reservoirs with a certain distribution law, firstly, based on the widely used Gaussian model, the heterogeneous porosity field is established according to the porosity range, average porosity, and standard deviation. Then, in order to characterize the permeability field corresponding to the heterogeneous porosity field, the complex pore structure characteristics are described according to the fractal characteristics. The more detailed the actual 5 Geofluids reservoir data obtained, the more realistic the established heterogeneous field will be. This process is also applicable to other heterogeneous reservoir methods, such as non-Gaussian field.

Heterogeneity of Matrix.
2.3.1. 3D Gaussian Filter. Gaussian filter is a linear smoothing filter, which is widely used in an image noise reducing process [38]. The specific operation of Gaussian filtering is to use a template to scan each element in the matrix (here refers to each mesh) and replace the central mesh value of the template with the weighted average value of the mesh value in the neighborhood determined by the template.
The heterogeneous field established in this work is obtained by convolution of a random normal distribution (template) with Gaussian filter. The 3D Gaussian function is where σ is standard deviation; the default value in this paper is 0.6. According to Equation (7), the heterogeneous porosity field can be obtained.

Fractal Permeability Model.
It is hard to establish a satisfactory model to accurately describe the relationship between the complex pore structure and physical properties of reservoir porous media. The Kozeny-Carman (KC) equation is a commonly used empirical formula to express reservoir permeability, which has been expanded and modified by many scholars [39,40]. In this paper, the permeability fractal model established by Zheng and Li [41], considering the influence of specific surface, is selected to characterize the heterogeneity of reservoir permeability.
where D is the fractal dimension of the porous media pore structure, which reflects the complexity of the pore structure; τ is tortuosity pedantic; λ = r min /r max is the ratio of the minimum and maximum pore radii in a porous medium, representing the difference in equivalent pore radii between tight matrix and natural microfractures; and S is the porous media specific surface (μm 2 /μm 3 ).

Model Validation
The open-source MATLAB Reservoir Simulation Toolbox (MRST) was used to solve the problem [31], and the coupling of automatic differential module, black oil model, EDFM fracture module, and fractal model was realized. The governing equations are discretized by the two-point flux approximate finite volume method. The backward Euler scheme is used for time discretization. Then, the NNC operator is embedded into the solution frame of the black oil model to realize flow simulation. Figure 2 shows the main framework and solution process. Firstly, in order to verify the accuracy of the proposed model, it is simplified to a single-phase homogeneous version for comparison with the modeling result in the literature [15]. The model for verification is a rectangular reservoir of 1000 × 500 × 10 m 3 . The horizontal well was subjected to constant pressure production after four-stage fracturing. The main parameters are consistent with those in the literature. Simulation results in Figure 3 show that the results of the model established in this paper agree with those of the conventional homogeneous model, with a comprehensive error of less than 5%. In unconventional reservoirs, horizontal well fracturing is used to perform depleted development with a low overall recovery, during  6 Geofluids which a rapid decline is shown in the early stage; a low production is maintained in the later stage. As a result, some enhanced recovery (EOR) techniques, such as water injection, gas injection huff, and puff, are also adopted for production [42]. However, the works reported in a relevant reference are obviously not applicable for they only take consideration of single-phase homogeneous model. The model presented in this work could be used for simulation of homogeneous and heterogeneous reservoirs, single-phase or multiphase flow and complex EOR technologies, showing a great application potential.

Results and Discussion
In this section, both homogeneous and heterogeneous reservoirs are considered, and there is a deep discussion on physical properties, saturation, in-layer heterogeneity, interlayer heterogeneity, and examples of EOR techniques (the  7 Geofluids injection and production of mixed five-point pattern of vertical and horizontal wells and adjacent horizontal wells, etc.) to illustrate the applicability, robustness, and performance of the model.

Homogeneous Reservoir
4.1.1. Different Permeabilities of Matrix. In order to illustrate the necessity of discussing the matrix heterogeneity, the influence of matrix permeability in homogeneous reservoirs is first discussed. Basic parameters are set as shown in Table 1. The diagram of horizontal well mesh, production curve, and pressure distribution at different times are shown in Figure 4. The effects of three different matrix permeabilities with a small difference on the production of fractured horizontal wells are discussed in homogeneous reservoirs. It can be seen from the results that although the matrix permeability difference is small, there is a significant difference in the production of horizontal wells during the early stage of rapid decline. The main reason is that the flow exchanges between fracture and matrix and the exchanges between fracture and the well bore are the source of energy supply in this stage. The greater the matrix permeability, the greater the production is in the early stage. There is a linear decline in the later stage, in which the main energy supply is from the matrix, and the difference gradually decreases. When the permeability difference is small, the adoption of a single permeability will cause a large error. When the natural microfractures in tight reservoirs are widely distributed, the limitations of the characterization by single permeability will be even more obvious; thus, the heterogeneous model with consideration of permeability difference is much more in line with the actual situation.

Different Saturations of Matrix.
Generally, the fluids in the reservoir coexist in multiple phases, including oil phase, formation water, and dissolved gas (not considered in this paper). During the actual production of oil wells, coproduction of water and oil will occur. Therefore, the irreducible water saturation set in this section is 35%, and the residual oil saturation is 25%. The oil saturation in the reservoir for comparison is 63% (scheme 1) and 60% (scheme 2), the matrix permeability of which is set to be 0:5 × 10 -3 μm 2 , and other parameters are set as shown in Table 1. The production obtained by simulations is shown in Figure 5.
The initial oil production rate of the two schemes is the same. The lower the oil saturation, the faster the initial production decline will be. So that the overall production of low oil saturation is lower. When the water saturation is higher than the irreducible water saturation, water production will occur during oil production. In scheme 1, the water saturation is 37%, which is slightly higher than the irreducible water saturation. In this case, the movable water is only in small amount and the water production is small, but with the increase of water saturation, the water production will significantly increase. The average water cut in scheme 1 is 10% while the average water cut in scheme 2 is 42%. In actual reservoirs, water production is also common [43,44]. The effects of saturation and relative permeability of oil and water on the simulation results cannot be ignored, especially in some EOR methods [45][46][47]. Thus, the model proposed can be applied to the complex situation of coproduction of water and oil during the actual production.

4.2.1.
In-Layer Heterogeneity. The in-layer heterogeneity is characterized according to Equations (7) and (8) and the Gaussian filter; the permeability field is established. The main parameters are set as shown in Table 2.
The relation between porosity and permeability under different ratios of the minimum and maximum pore radii in porous media ðλ = r min /r max Þ is shown in Figure 6. r min is taken as the pore radius of tight matrix, and r max is taken as the equivalent pore radius containing natural microfracture characteristics. The stronger the heterogeneity is, the smaller the value of λ, and the higher the permeability. The K max /K min in Table 2 stands for permeability difference, the value of which for the four types in the examples can reach about 10 times or even higher, so it can be reasonably adjusted according to the actual reservoir situation. 3D permeability distribution field generated by the parameters of set 1 and set 3 in Table 2, and the corresponding probability histogram is shown in Figure 7.
The production obtained by simulations corresponding to sets 1, 2, and 3 in Table 2 is shown in Figure 8, with an initial oil saturation of 80%. When λ = 0:08, the permeability field conforms to the situation that the microfractures are sparsely distributed in heterogeneous reservoirs (the overall permeability is less than 1 × 10 −3 μm 2 ), the initial production of which is low, showing a near-linear decline, ending up with a low accumulative production. When λ = 0:04 and 0:03, the distribution range of permeability is wide, suggesting a strong heterogeneity due to wide distributed fractures. The initial production under these circumstances is high, but the decline in production is quite rapid, showing that the influence caused by heterogeneity 8 Geofluids on the production of horizontal well is large in the earlier stage and smaller in the later one.

Interlayer
Heterogeneity. The heterogeneity of monolayer only presents the heterogeneity within one layer, but the heterogeneity in vertical direction cannot be neglected either. Interlayers exist in the vertical direction of actual reservoirs (porosity and permeability of which are significantly lower than the upper and lower layers). Effect on production brought by the following two kinds of interlayers is discussed in this section: (a) The vertical direction of reservoir is divided into 10 layers of meshes, among which the 1 st~3rd and 8 th~1 0 th layers have strong heterogeneity with nature microfractures widely distributed, while the 4 th~7th layers are relatively tight, which altogether is considered the interlayer (b) The vertical direction of reservoir is divided into 10 layers of meshes to represent 5 heterogeneous layers that each contains two layers of meshes. The 1 st , 3 rd , and 5 st layers in the vertical direction have strong heterogeneity with nature microfractures widely distributed, while the 2 nd and 4 th layers are considered the tight layers The permeability distribution of the above two situations and the production obtained accordingly by simulations are shown in Figures 9 and 10.
Although the reservoir simulated in case (a) has only one interlayer, the thickness of which is relatively large and concentrated. An overall near-linear decline shows in the production with a small initial production. The total proportion of the two interlayers accounting for all the layers in case (b) is the same as that in case (a), but since the interlayers are divided into two layers and are alternately arranged with 9 Geofluids strong heterogeneous layers, the initial production obviously increases. The later stage of production of these two shows a little difference.

Waterflooding Simulation Based on Five-Point Well
Pattern. For the tight reservoir with low porosity and low permeability, the recovery percent is low and the production 11 Geofluids drops rapidly by using depleted development of horizontal well. Water injection and gas injection are commonly used to improve the productivity.
In this section, discussions are made on establishing a multiscale fracture system, which contains microfracture, fracture in large scale, and artificial fractures, based on the fractal multilayer heterogeneous reservoir simulated. Four injection wells are set at the four corners (taking water injection as an example, with an injection rate of 20 m 3 /d). And the productivity of horizontal wells with external energy supply is explored by using a five-point well pattern. The meshes of injection and production, permeability distribution, and the water saturation when the initial saturation is 80% after producing for 300 days are shown in Figures 11(a) and 11(b), respectively. It is shown that the differences in the front edge of water flooding are large in the four injection wells due to the heterogeneity. When the water reaches the fracture in large scale or artificial fracture, it will quickly flow into the fractures, which will directly affect the water breakthrough time of the horizontal wells. Figure 11(d) shows the production curves under different saturations. Compared with the depleted developing, the production curve is hump-shaped due to the injection of external fluid to replenish formation energy. When the initial oil saturation is 80%, there will be a water-free period, the length of which depends on whether the water communicates with the production well by the fracture or fault. In this case, fractures are widely distributed; the water cut rises in an S-shaped curve.
When the initial oil saturation is 60%, water production occurs from the well opening, and the water cut reaches to around 43%. Although the production curve is also humpshaped, it is obviously lower and turns to decline after reaching the peak, showing no water-free period in the water cut curve.

Injection and Production of Adjacent Horizontal Wells.
The cases above are examples using constant pressure production, in this section, three adjacent horizontal wells are discussed. The middle one is injection well with a rate of 50 m 3 /d, and the other two are in constant production mode with a production rate of 25 m 3 /d. In this case, the main purpose to compare the effects of heterogeneity on injectionproduction connectivity.
Also, two interlayers are arranged. In Figure 12, the meshes are shown in (a); the permeability is shown in (b), suggesting a red section with relatively high permeability between the injection well and production well 2; the same section has a relatively higher water saturation shown in (c), suggesting a better connectivity between the injection well and production well; the production curves obtained by simulation are shown in (d); it could be directly observed that the water breakthrough time of production well 2 is earlier than that of production well 1. And because of the influence by heterogeneity, the decline of oil production and the rise of water production in curves are fluctuant but not as smooth as those in the previous cases. This illustrates the necessity to consider heterogeneity in the simulation for well patterns, for the communication of the fluid injected will be further accelerated if there exists a fracture zone with connectivity to certain extent.

Highly Heterogeneous Reservoir.
The previous examples show a wealth of cases for unconventional reservoirs, and this section focuses on how to build reservoirs with highly heterogeneous reservoirs with differences in permeability that can reach more than 10 6 , which is a common problem in real reservoirs. It is well known that reservoir information is usually limited. Assuming that information about porosity is known, as in the case of SPE10, and the proposed method was used to establish the permeability field.
As shown in Figure 13(a), there is a significant difference in physical properties between the upper and lower layers of the reservoir, which means a significant difference in permeability. In order to understand the fractal characteristics of the reservoir, a preliminary permeability field is first established for the whole reservoir through Equation (8). It is worth noting that we have segmented the permeability field, which means that the fractal characteristics of permeability field are different for different porosity intervals as shown in Figure 13(b). There is still a large error with the actual. Then, considering the difference of lithology, we also carried out stratified and lithologic permeability field reconstruction. This means that the fractal characteristics of different layers and different lithologies are different. Then, several adjustments are made to finally form the permeability field shown in Figure 13(d) that is close to the actual permeability field. Looking back at this process, the more information we have, the better it is to create realistic heterogeneous reservoirs. However, when the target reservoir information is missing and only has local characteristics, the fractal parameters are obtained by fitting, and the porosity and permeability fields of the unknown region are established by using the method in this paper. Of course, the method in this paper also has limitations. The current method is mainly the unknown Gaussian field, but this process is also applicable to other methods. It is necessary to select a suitable method for the target reservoir. The more information we have about highly heterogeneous reservoirs, the closer we get to the actual geological features.

Conclusion
Based on the model in this paper, the reservoir heterogeneity, saturation, five-point well pattern of mixed vertical and horizontal wells, injection and production of horizontal wells, and highly heterogeneous are discussed. The results show that formation fluids coexist in multiphase and the water saturation has a direct effect on the production; the model proposed is suitable for the complex situation of coproduction of oil and gas during actual production and could exert much more value in the injection and production technology.
A rapid decline occurs during the development and depletion of horizontal wells. The mixed five-point pattern of vertical and horizontal wells can improve the production, but the dual effects of heterogeneity and fracturing will result in cross-communication of injected fluids, leading to a rapid decline in productivity.
The interwell injection-production connectivity influences the numerical simulation of injection-production well patterns. Near-well natural microfractures are opened because of fracturing of horizontal wells, and the heterogeneity cannot be ignored, which directly affects the time of water invasion.
However, the characterization of reservoir, especially the 3D parameters of fractures and the complex physical mechanics of unconventional reservoirs, has always been a challenge. This model is only suitable for some cases, and heterogeneity is not fully considered. Productivity prediction in different reservoir scenarios remains to be explored. Accurate productivity forecasting is an interesting and challenging subject that needs to be worked on.

Data Availability
All data included in this study are available upon request by contact with the corresponding author.