Abstract
Understanding the relationship between petroleum recovery and characteristics of hydraulic fracture network is a key component of economic development of tight reservoirs. Owing to the limitations inherent in current reservoir simulators, optimization of fracture network has been simply focused on the parameters of fracture conductivity, fracture number, aperture, and so on. Deeper insight into the effect of decisive parameters, such as fracture density and fracture connectivity on the well production in tight reservoirs, is now required to maximize the petroleum recovery. In this work, a newly developed discrete fracture simulator is applied to comprehensively study the effect of fracture density and fracture connectivity in tight reservoirs. Conceptual models with different fracture densities and different fracture connectivity are firstly designed and simulated to explore how these two parameters affect the reservoir behavior and establish the equations for effect measurement. Then, we simulate models with different well placement strategies and a fixed set of natural fractures to determine the optimal strategy. Finally, simulations are performed on a field-scale reservoir with three long fractured horizontal wells. Results demonstrate that increases in either fracture density or fracture connectivity can significantly improve well production. However, an optimal value exists considering the economic profit. Compared to the fracture density, fracture connectivity plays a more important role in affecting the well production. In a tight reservoir with abundant natural fractures, making the horizontal well parallel to the direction of natural fractures is determined to be the optimal well placement strategy. The heterogeneous distribution of remaining oil in real tight oil reservoirs is mainly caused by the heterogeneous distribution of fracture density and fracture connectivity.
1. Introduction
Nowadays, tight reservoir provides a significant fraction of petroleum production around the world [1–3]. Extraction of fluids from tight reservoirs mainly relies on the large-scale hydraulic fracturing which creates complex fracture network in the reservoir. This fracture network provides significant conduits for fluids to flow at economic rates while bringing great challenges to the reservoir simulation [4, 5]. Since reservoir simulation is an important tool of optimizing stimulation design, completion practice, and development strategy, determination of the simulation-based solution for optimal development of tight reservoir also faces a great challenge. Fracture network optimization is the primary issue of improving recovery of tight reservoir. Understanding the relationship between fracture network parameters and petroleum recovery will greatly benefit the economic development of tight reservoir [6]. Long and Witherspoon were among the first to investigate the significant effect of fracture network parameter on reservoir performance [7]. In their study, through utilizing a numerical code [8] developed to determine permeability of a fractured system, they demonstrated that as degree of fracture interconnection increases, the permeability of the fractured system increases. In current commercial reservoir simulators, a dual-continuum model, known as the sugar-cube model, is widely used to perform optimization for fractured reservoir. In this model, the fractured rock is generally idealized as an equivalent continuum medium through assuming that the matrix is partitioned by an orthogonal fracture network [9–11]. Except for the recognized shortcomings that are inherent in a dual-continuum model, namely, inadequacy of modeling reservoir with a limited number of flow-dominant fractures [4, 12] and having difficulty in correctly evaluating the exchange terms between the matrix and the fractures [13, 14], the dual-continuum model also cannot be used to optimize some specific fracture parameters, such as fracture density and fracture connectivity because of the fixed configuration of fracture network in the models. Hence, in the past decades, optimization of the fracture network was limited to simple parameters such as fracture conductivity, fracture number and fracture aperture. In 2006, Mayerhofer et al. performed a parametric study using a commercial simulator to show how fracture network size and density, fracture conductivity, matrix permeability, and gaps in the network affect well productivity [15]. Warpinski et al. conducted similar studies, and simulation results demonstrated that shale reservoir with ultralow permeability requires an interconnected fracture network of moderate conductivity with a relatively small spacing between fractures to obtain reasonable recovery factor [16]. Through performing simulations with CMG-IMEX, the economic optimal fracture conductivity and well number for the cases with and without geomechanics in Bakken tight reservoir were determined by Yu and Sepehrnoori [17]. In the study of Saputelli et al., simulation results obtained by Nexus revealed that an optimal number of fractures should be determined in tight reservoir since increasing number of fractures will not always improve the short-term economics [18]. Optimization of these fracture network parameters is critical because of the high cost of drilling and fracturing treatment. However, deeper insights into the relationships between more fracture network parameters and the reservoir performances are required for the best understanding of the tight reservoir development.
Recently, a discrete fracture model (DFM) in which the fractures are represented explicitly and individually rather than idealized as continuous media was developed. The DFM makes it possible for reservoir engineers to optimize complex fracture parameters. However, development of DFMs is still at its early stage. Currently, there are two typical categories of DFM: embedded discrete fracture model (EDFM) and unstructured discrete fracture model (UDFM). In the EDFM, a conventional structured mesh is used to represent the matrix while the fractures are discretely embedded in the matrix [19, 20]. And in the UDFM, fractures are represented by discretely distributed entities and the matrix is partitioned in an unstructured manner conforming to the geometry of fracture network [13, 21, 22]. The EDFM is currently not efficient enough when handling large-scale complex fracture network owing to the boundary element approach used to evaluate the effective matrix permeability associated with the small fractures. The evaluation process becomes expensive while the number of fractures increases [23]. The UDFM currently can only simulate incompressible fluid flow problems, which means energy supplement strategy, such as water or gas injection, must be applied in the simulations. However, energy supplement is currently not practical in tight reservoirs since the matrix in these reservoirs is extremely tight and the fractures can bypass oil in the matrix, leading to a very low efficiency of the energy supplement strategy [24]. Oil extraction from the tight reservoirs is currently mainly relied on the natural formation energy and the mechanisms of rock and fluid expansion [25–27]. Hence, rock and fluid compressibility are critical factors that need to be taken account in the simulators developed for tight reservoirs.
Owing to the limitations of current commercial reservoir simulators and DFM, a comprehensive study of the optimizations of fracture density and fracture connectivity in the large-scale complex fractured network based on reservoir simulation has not been reported in the literature so far. In this study, equations characterizing the fracture density and fracture connectivity are firstly established. Then, a newly developed in-house discrete fracture simulator (SC-CFR) which is capable of modeling compressible rock and fluid and handling large-scale complex fracture networks is used to perform reservoir simulations to quantificationally study the effect of fracture density and fracture connectivity on the reservoir behavior [28]. Some novel physical models including conceptual models and large-scale complex fractured models analogous to the real tight reservoirs are designed for the optimizations of fracture density and fracture connectivity. Through analyzing simulation results of these models, we demonstrate the crucial or even decisive effect of fracture density and fracture connectivity on the tight reservoir behavior. Functions for effect measurement are also proposed. Results presented in this work can help reservoir engineers gain an insight into the optimization of fracture network, which will greatly contribute to the improvement of future horizontal well completion and fracturing strategy and the resultant maximization of economic benefit of tight reservoir development.
2. Description of the Simulator
The framework of the SC-CFR is shown in Figure 1. The SC-CFR is mainly composed of five independent but interrelated modules. The edge-constrained Delaunay Triangulation scheme in which fractures are represented individually and explicitly with a lower dimension than the matrix is utilized to grid the physical domain [29, 30]. A transmissibility list that includes connections of matrix-matrix, matrix-fracture and fracture-fracture is evaluated by a two-point flux-approximation (TPFA) method [21]. Spatial discretization of the mathematical model is greatly simplified by three discrete operators that are defined to represent the expressions of divergence, gradient, and average as their original forms and, thus, enable rapid prototyping of the numerical model. The nonlinear discretized equation system is solved by modified Newton’s method in which the Jacobian matrix is computed by the automatic differentiation (AD) technique. Application of the AD can greatly enhance the efficiency and accuracy of simulation models, especially for the models with complex grid systems [31–33].

The concept of embedded well which is inspired by the embedded discrete fracture model is defined in our simulator to establish the well model (Figure 2 in [19]). In this model, a horizontal well is treated as a fracture, which is embedded in the background grid. The flux between grid and the well is determined by the transport index (TI) [19, 20]. The geometry of the well can be very flexible in the embedded well model since the well lines do not need to conform to the grid edges.

(a)

(b)

(c)

(d)
3. Effect of Fracture Density
As shown in Figure 3, four physical models with different fracture densities of the fracture networks are established. Sizes of the models are . To minimize the influence of fracture connectivity, fracture networks are created in the same pattern; namely, fractures are only distributed either in -direction or -direction, and fractures in each direction are equally spaced. Sizes of the SRV in the four models are set to be the same. We define the fracture density herein as where FD is fracture density of fracture network (m/m2); is sum of lengths of all fractures in the SRV (m); is area of the SRV (m2).

(a)

(b)

(c)

(d)
Simulation parameters for the four models are the same and are listed in Table 1. The matrix permeability and matrix porosity are set to be as low as 0.01 mD and 8%, respectively, in order to capture the characteristics of the in tight reservoirs [34].
We simulate the four models by using SC-CFR. Figure 4 shows the comparison result of well performances of the four models. In tight reservoirs, fractures are the main channels for fluids to flow from the reservoirs to the wells. Therefore, it can be observed in Figure 4 that the well productivity increases significantly with the increasing fracture density. Increasing FD from 0.05 to 0.07 leads to 6.5% growth in 200-day cumulative production. However, Figure 4 also depicts that after the fracture density exceeds 0.09, the increase of production become much slighter. Increasing FD from 0.09 to 0.11 only leads to 1.2% growth in 200-day cumulative production. This is mainly because the reserve in the SRV is limited, and the fluids in the area beyond SRV cannot flow into the SRV owing to the extremely tight matrix (Figure 2). Considering that the fracture network with bigger fracture density requires bigger cost of hydraulic fracturing, the optimization of fracture density is suggested.

In addition, the production rate curves in Figure 4 illustrates that the significant influence of fracture density on the well production only occurs in the early period (<60 days) of reservoir development. After that time, the well production is mainly provided by the matrix in the SRV; hence, the fracture density has little effect on the well production and the well produces in a low rate for a long period.
To validate the interpretations above and directly depict the reservoir behaviors of models with different fracture densities, reservoir pressure profiles at a certain simulation time (50 days) for the four models are obtained and shown in Figure 2. As can be seen, in cases of , 0.07, and 0.09, the extent of pressure drawdown in the SRV increases dramatically with the increasing fracture density. However, the pressure in the SRV of is almost already at its lowest value and widest extend, which means increasing fracture density will not make a big difference. Therefore, when FD increases from 0.09 to 0.11, the change in pressure profile is less noticeable. In all the four models, pressures of the areas beyond SRV keep at their initial values.
The relation equation for the effect measurement of fracture density is obtained by curve fitting, as shown in Figure 5. As can be seen, the cumulative production over fracture density shows the law of logarithmic equation: where is cumulative production (m3).

4. Effect of Fracture Connectivity
In this section, we first investigate how the fracture connectivity impacts the reservoir behavior through simulating four conceptual models. Then, we take the well placement optimization as an example to show the important role that fracture connectivity plays in field development. The fracture connectivity is defined as where is number of nodes at fracture intersections and is the number of fracture endpoints.
4.1. Conceptual Models
The four physical models for simulation of fracture connectivity are shown in Figure 6. Each model contains 13 fractures, and fracture No. 1 (F1) can be treated as a dominant fracture which is created by hydraulic fracturing of the well, and other fractures can be treated as nature fractures. From model (a) to model (d), more and more natural fractures intersect with F1. And in model (d), some nature fractures intersect with each other. Simulation parameters herein are the same with parameters listed in Table 1.

(a)

(b)

(c)

(d)
The obtained well productions of the four models are compared in Figure 7, from which we can observe that the bigger the fracture connectivity, the higher the well production (158.9% and 49.2% increase in 200-day cumulative production when increasing FC from 0 to 0.15 and from 0.15 to 0.31, respectively). However, when increasing FC from 0.31 to 0.62, the increase in cumulative production (5.1%) significantly decreases. This is also because of the limited reserves in the SRV.

Through comparing Figure 4 with Figure 7, we can find that the fracture connectivity plays a more important role in impacting the well production than the fracture density does. Even though the fracture densities in model (a) and model (d) are almost the same, well production dramatically differs from each other because of the different fracture connectivities.
Figure 8 shows the reservoir pressure profiles of the four models obtained at the last step of the simulation. The regions of pressure drawdown in these profiles differ dramatically from each other although sizes of SRV and numbers of fracture in SRV are almost the same. In model (a), none of the natural fractures contributes to improving the well production since they disconnect with F1. We characterize the fracture that has no interconnections with other fractures as dead fracture (DF) in this study. As can be seen from Figure 8, with the decreasing number of DFs, the region of pressure drawdown becomes larger and larger; namely, more and more fluid is produced from the SRV.

(a)

(b)

(c)

(d)
It can be concluded from the interpretations above that the fracture connectivity almost plays a decisive role in the development of fractured tight reservoirs. Hence, in real reservoir development, optimization of fracture connectivity in the SRV should be paid much attention to.
The relation equation for the effect measurement of fracture connectivity is obtained by curve fitting, as shown in Figure 9. As can be seen, the cumulative production over fracture connectivity shows the law of power equation:

4.2. Well Placement Optimization considering Fracture Connectivity
In real reservoirs, natural fractures are preexisting fractures, which means their orientations cannot be adjusted as what we did in Figure 6. Hence, the optimization of fracture connectivity should be focused on optimizing parameters of hydraulic fractures and well placement. In this section, we consider the hydraulic fractures as multistage fractures that are perpendicular to the horizontal well, and numbers and lengths of the hydraulic fractures are considered constants. Under these assumptions, optimization of the well placement becomes our main focus.
To capture the densely distributed natural fractures in tight reservoirs, three physical models with the same set of nature fracture which is composed of 400 fractures (randomly distributed at different positions) are established. Typically, in the real reservoirs, fractures of a same fracture set propagate in the same direction owing to the fixed directions of principal stresses [35, 36]. Therefore, all natural fractures in our models are placed at a 45-degree angle from the -direction. The three physical models are shown in Figure 10. Horizontal wells in the models are placed at different directions. Suppose the angle between the horizontal well and the natural fractures is , then in models (a), (b), and (c), °, 45°, and 0°, respectively. Simulation parameters are also the same with the parameters listed in Table 1.

(a)

(b)

(c)
Well performances of the three models differ significantly from each other, which can be observed in Figure 11. With the decreasing , well production increases dramatically. The reason can be interpreted by Figure 10, in which the natural fracture network is simplified as equally spaced lines. As can be seen, in scenario (a), none of the natural fractures intersect with the hydraulic fracture; hence, fluid can be only produced from a very small region around the hydraulic fracture, resulting in a low production. In scenario (b), each hydraulic fracture intersects with three natural fractures, resulting in a much higher well production. And in scenario (c), each hydraulic fracture intersects with four natural fractures, resulting in the highest production. The FC values obtained by Equation (2) are 0, 0.5, and 0.67 for scenarios (a), (b), and (c), respectively (Figure 12).


(a)

(b)

(c)
The reservoir pressure profiles of the three models obtained at the last step of simulation are shown in Figure 13, which fully validates the interpretations above. With the decreasing , the size of region of pressure drawdown increases significantly, resulting in the increase in well production.

(a)

(b)

(c)
According to the above simulation results and considering that, in Figure 12, no matter what direction the horizontal well is placed, each hydraulic fracture can at most intersect with four natural fractures, we conclude that for densely fractured tight reservoirs, the optimal direction of the horizontal well is being parallel to the direction of natural fractures.
5. Simulation of Large-Scale Densely Fractured Model
In this section, a field-scale model containing three long horizontal wells is established. Each well is stimulated by large-scale hydraulic fracturing that creates complex fracture network around the well. The fracture densities of the fracture networks are set to be different from each other. Through performing simulation on this model, we aim to illustrate how the fracture density and fracture connectivity significantly impact the development of real tight reservoirs.
Figure 14 shows the grid of our model. The size of the model is . The length of each horizontal well is 2600 m. The number of fractures around Well-A, Well-B, and Well-C is 350, 250, and 150, respectively. The fractures are randomly distributed around the wells. The SRV of the three fracture networks is the same, namely, . The fracture density and fracture connectivity around the three wells are calculated by using equations (1) and (3). The density of fractures around Well-A, Well-B, and Well-C is 0.035 m/m2, 0.025 m/m2, and 0.016 m/m2, respectively. The connectivity of fractures around Well-A, Well-B, and Well-C is 0.573, 0.424, and 0.223, respectively. Simulation time herein is days. Other simulation parameters are the same with the parameters listed in Table 1.

The obtained well productions are compared in Figure 15. As can be seen, the trend of production variation with the increasing fracture density is similar with the trend depicted in Figure 4. Increasing the number of fractures from 150 () to 250 () leads to a production increase of 210.2 m3, which is much greater than the production increase resulting from the increase of fracture number from 250 to 350 (). Figure 16 shows the time evolution of reservoir pressure profile. As can be seen, the region of pressure drawdown around each well extends little with time and does not contact with each other, indicating the negligible impact of interwell interference on well production. This phenomenon is attributed to the extremely tight matrix in tight reservoirs.


(a) days

(b) days

(c) days
It can also be observed in Figure 16 that pressure distribution in the SRV of every well is highly heterogeneous. In the traditional sand reservoirs, the main reason will be directed to the reservoir heterogeneity which is caused by the heterogeneous distribution of rock permeability. However, in the densely fractured tight reservoirs, this phenomenon is attributed to the heterogeneous distribution of fracture density and fracture connectivity. Take the region (RA) highlighted by blue dotted line around Well-A (Figure 14) and the region (RB) highlighted by black dotted line around Well-B for example; in RA, fracture density is considerably smaller than that of any other regions around Well-A, resulting in the much smaller pressure drawdown of this region. In RB, although the fracture density is as big as that of other regions around Well-B, the pressure drawdown of this region is still smaller than that of any other regions. This is because most fractures in RB are laid in the same direction (-direction), resulting in an unfavorable fracture connectivity. Since RA and RB are bypassed zones and typically, it is very difficult to extract the remaining fluids from the bypassed zones through conventional technics; it is suggested that the fracturing treatments should be optimized to make the distribution of fracture density and fracture connectivity in the SRV as uniform as possible.
6. Conclusions
A new workflow based on reservoir simulation is designed to quantificationally study the effect of fracture density and fracture connectivity on the well performance in tight reservoir. The simulations are performed using an in-house discrete fracture simulator which is capable of modeling compressible rock and fluid and handling complex fracture networks. The results indicate that fracture density and fracture connectivity are two factors that can dramatically affect reservoir behaviors of densely fractured tight reservoirs. Equations for effect measurement are proposed. Optimizations of these two parameters should be taken into account in the stimulation design. For the optimization of fracture density, larger scale of hydraulic fracturing does not always lead to greater economic benefit. An optimal fracture density needs to be determined when considering the high cost of hydraulic fracturing. The primary issue of optimizing fracture connectivity is to identify the propagation direction of natural fracture since the optimal direction of the horizontal well is being parallel to the direction of natural fractures. In tight reservoirs, we reveal that the heterogeneous distribution of remaining oil is mainly caused by the heterogeneous distributions of fracture density and fracture connectivity. Hence, we suggest that the fracturing treatments should be optimized to make the distribution of fracture density and fracture connectivity in the SRV as uniform as possible.
Data Availability
All data included in this study are available upon request by contact with the corresponding author.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research is supported by the NSFC Basic Research Program on Deep Petroleum Resource Accumulation and Key Engineering Technologies (U19B6003) and Sinopec Research Project (P19031-2).