Integrating Embedded Discrete Fracture and Dual-Porosity, Dual-Permeability Methods to Simulate Fluid Flow in Shale Oil Reservoirs
The oil recovery factor from shale oil remains low, about 5 to 7% of the oil in place. How to increase oil recovery from shale oil reservoirs is attracting more and more attention. CO2 huff-and-puff was historically considered one of the best approaches to improve the oil rate.
Weirong Li1, Zhenzhen Dong2 and Gang Lei1
1Beijing Innovation Center for Engineering Science and Advanced Technology (BIC-ESAT), Peking University, Beijing 100193, China. 2Petroleum Department, Xi’an Shiyou University, Xi’an 710065, China.
Received: 8 August 2017 / Accepted: 18 September 2017 / Published: 23 September 2017
Most previous simulation studies have been based on dual porosity, but simulation results from dual-porosity models have not been as accurate as discrete fracture models in composition modeling.
This study proposes a new model that integrates the embedded discrete fracture model and the dual-porosity, dual-permeability model (DPDP). The newly developed method could explicitly describe large-scale fractures as flow conduits by embedded discrete fracture modeling and could model the flow in small- and medium-length fractures by DPDP modeling. In this paper, we first introduce four different non-neighboring connections and the way to calculate the transmissibility among different media in the new model.
Then, the paper compares the performance of the new method, discrete fracture modeling, DPDP modeling, and embedded discrete fracture modeling for production from oil reservoirs. Following, the paper carries out a series of simulations to analyze the effects of hydraulic fracture stages, hydraulic fracture permeabilities, and natural fracture permeabilities on the CO2 huff-and-puff process based on the new method. In addition, the injection cycle and soaking time are investigated to optimize CO2 huff-and-puff performance.
This study is the first to integrate embedded discrete fracture modeling with DPDP modeling to simulate the CO2 huff-and-puff process in a shale oil reservoir with natural fractures. This paper also provides detailed discussions and comparisons on the integrated strategy, embedded discrete fracture modeling, discrete fracture modeling, and dual-porosity, dual-permeability modeling in the context of fracture simulation with a compositional model. Most importantly, this study answers the question regarding how fractures affect CO2 huff-and-puff and how to optimize the CO2 huff-and-puff process in a reservoir with natural fractures.
Complex fracture networks in shale oil reservoirs make the accurate modeling of oil production in these reservoirs very challenging. Three classes of models, equivalent continuum models (ECMs), dual-continuum models, and discrete fracture models (DFMs) are applied to model fractured systems.
Based on the equivalent continuum theory, Snow  first proposed the ECM which describes flow in fractured media mathematically. Since the introduction of this model, several scholars have conducted research on this method, such as Oda  and Tian .
Dual-continuum modeling is the conventional method for simulating fractured reservoirs and is widely used in the industry. In the 1960s, Barenblatt et al.  first established dual-porosity model to simulate matrix-fracture flow systems. Afterward, Warren and Root  developed a more complete dual-porosity model. Kazemi et al.  and Saidi  developed dual-porosity simulators, which extended the Warren and Root approach to multiphase flow.
Recently, Pruess and Narasimhan  subdivided matrix rock gridding and proposed a multiple interacting continua (MINC) model in 1985. Furthermore, Moinfar et al.  proved that the dual-continuum model could not provide accurate solutions in the high-localized anisotropy and presence of large-scale fractures in 2012.
The DFM method, in which an element and a control volume explicitly represents each fracture by, was established to model fluid flow in individual fractures and provide more realistic representations of fractured reservoirs than dual-continuum models. Most DFMs honor the geometry and the location of fracture networks with relying on unstructured grids. Kim and Deo  adopted a finite element method to perform discretization for DFM and to combine matrix and fractures based on the superposition principle.
Later, Karimi-Fard and Firoozabadi  applied DFM to solve the two-phase flow problem in the fractured media. Likewise, Monteagudo and Firoozabadi  and Matthäi and Belayneh  used control-volume finite element methods for multiphase flow in fractured media and developed numerical simulators. Karimi-Fard et al.  extended DFMs compatible with multiphase reservoir simulators based on an unstructured control-volume finite difference formulation.
In addition, Li and Lee  and Moinfar et al.  developed embedded discrete fracture models (EDFMs), which use a structured grid to represent the matrix and introduce additional fracture control volumes by computing the intersection of fractures with the matrix grid. Recently, more attention is paid to the EDFM because of its flexibility and accuracy. Hadi Hajibeygi et al.  proposed the iterative multiscale finite volume (i-MSFV) approach to improve the efficiency and the accuracy for fractured porous media, and the results are very promising.
Ţene et al.  proposed a projection-based EDFM to deal with cases where the fracture permeability lies below that of the matrix. The fracture-crossflow-equilibrium method was proposed in 2017 by Ali Zidane and Abbas Firoozabadi  to study compositional two-phase flow in the fractured media; simulation results showed that central processing unit (CPU) time with fractured media and without fractured media differs by less than a factor of two. Sander Pluimers  introduced hierarchical fracture modeling (HFM) in 2015, upscaling small- and medium-scale fractures into an effective matrix permeability and studying upscaling criteria for the HFM.
In order to balance accuracy, computational efficiency, and field practice, we propose a new method that integrates EDFM and dual-porosity, dual-permeability (DPDP) concepts to model the production process in shale oil reservoirs. The developed method could explicitly describe large-scale permeable fractures as flow conduits and simulate natural fracture networks that connect the global flow in stimulated areas of shale oil reservoirs. It may take years to reach the pseudo-steady state in the matrix systems, so the traditional dual-porosity approach can result in large inaccuracies.
With the newly developed simulator, we perform comprehensive simulation studies to determine the key factors of fractures and reservoir which affect the ultimate oil recovery in shale oil reservoirs. Different engineering factors and injection strategies are also analyzed and compared to provide guidance for production optimization during the CO2 huff-and-puff process.
Naturally fractured reservoirs typically have a wide range of fracture length-scales, ranging from micrometers up to several hundred meters. Figure 1 shows the correlation between fracture density and length. The small- to medium-scale fractures usually have higher densities than large-scale fractures in the reservoir. These fractured reservoirs are geologically too complex to be fully represented by a DFM because it individually represents each fracture.
Figure 1. Fracture density and fracture length correlation.
Therefore, a new model is developed combining the EDFM with the DPDP model to obtain a simplified reservoir model. As shown in Figure 2, large-scale fractures are assumed to be main fluid conduits and to have a big impact on flow pattern. For this reason, large-scale fractures are kept explicitly in the EDFM to model them with great accuracy. Medium-scale fracture lengths are about one to five times of the grid size and are modeled by the dual-porosity model. Small-scale fractures are expected to have a small impact on flow field. These fractures are approximated by effective “damaged matrix rock” properties, which are obtained by upscaling methods. This way, the small and medium fractures are deducted from the EDFM, resulting in a big reduction in the simulation effort.
Figure 2. The process to build dual porosity dual permeability (DPDP) model and embedded discrete fracture +DPDP model (EDFM+DPDP).
The key aspect of the new model is the calculation of connection transmissibility for flux interaction among different continua. Non-neighboring connections (NNCs) are subsequently determined as adding cells to represent fracture segments. Figure 3 illustrates the connection list of continua in the computational domain for a simple scenario. We adopted and extended the definition of NNCs, Four types of NNCs, which were adopted and extended from Moinfar , are defined as follows:
NNC type I: connection between a fracture segment and the fracture cell in the dual-porosity model;
NNC type II: connection between fracture segments in an individual fracture;
NNC type III: connection between intersecting fracture segments;
NNC type IV: connection between a fracture segment and the matrix cell in the dual-porosity model.
Figure 3. The example of connection list for EDFM+DP method.
NNC Type I: Connection between a Fracture Segment and the Fracture Cell in the Dual-Porosity Model
The NNC transmissibility factor between the fracture segment and fracture cell in the dual-porosity model depends on the fracture geometry and fracture permeability. As a fracture segment fully penetrates a matrix cell, the fracture-fracture cell transmissibility factor is:
where Affe is the area of the fracture segment on one side and kffe is the permeability tensor, which is calculated as:
kf is the embedded fracture permeability tensor, kfe is the fracture cell permeability tensor, and df-fe is the average normal distance from fracture to embedded fracture, which is calculated as
where V is the volume of the fracture cell, dV is the volume element of matrix, and xn is the distance from the volume element to the fracture plane.
NNC Type II: Connection between Fracture Segments in an Individual Fracture
The transmissibility factor between a pair of neighboring segments 1 and 2 is calculated using a two-point flux approximation scheme as:
where kfei and kfej are the fracture permeability, Ac is the area of the common face for these two segments, and dseg1 and dseg2 are the distances from the centroids of segments 1 and 2 to the common face, respectively.
NNC Type III: Connection between Intersecting Fracture Segments
It is very challenging to model fracture intersection accurately and efficiently for DFM due to the complexity of flow behavior at the fracture intersection. In 2014, Moinfar et al.  simplified this problem. They approximated the mass transfer at the fracture intersection by assigning a transmissibility factor between intersecting fracture segments. The transmissibility factor is given as:
where Lint is the length of the intersection line and df1 and df2 are the weighted average of the normal distances from the centroids of the subsegments (on both sides) to the intersection line.
NNC Type IV: Connection between a Fracture Segment and the Matrix Cell in the Dual-Porosity Model
The matrix permeability and fracture geometry determine the NNC transmissibility factor between matrix and fracture segment. When a matrix cell is fully penetrated by a fracture segment, the matrix-fracture transmissibility factor is:
where Amfe is the area of the fracture segment on one side and kmfe is the permeability tensor, which is calculated as:
kfe is the fracture permeability tensor, km is the matrix permeability tensor, and df-m is the average normal distance from matrix to fracture, which is calculated as:
where V is the volume of the matrix cell, dV is the volume element of matrix, and xn is the distance from the volume element to the fracture plane.
If the fracture does not fully penetrate the matrix cell, the pressure distribution in the matrix cell may deviate from previous assumptions. As a result, the calculation of the transmissibility factor is complex. To make the method non-intrusive, it was assumed the transmissibility factor is proportional to the area of the fracture segment inside of the matrix cell.