Abstract
Carbonate reservoirs in the South China Sea mostly contain natural fractures with various length scales and different intensities, which causes great challenges in efficient reservoir modeling and flow simulation. Existing efforts based on dual-porosity and dual-permeability models could not reflect the characteristics of production data in certain wells. To accurately and efficiently characterize multiscale fractures, a hybrid fracture characterization method is proposed. Firstly, fractures are divided into two types according to the geometrical size and interpretation approach. Then, small-scale fractures, characterized mainly by image log interpretations, are modeled by the traditional dual-porosity/dual-permeability (DP) method. And large-scale fractures, which are characterized by seismic interpretations and dominate the flow regime, are modeled by the embedded discrete fracture method (EDFM) to achieve both accuracy and efficiency. Lastly, transmissibilities among these three types of grid mediums are calculated to generate the hybrid DP+EDFM model for flow simulation. The proposed approach is applied to a carbonate, fractured reservoir in the South China Sea. The overall procedure is fast and reliable, and water cut matches of both field and specific wells are dramatically improved. Comparing the simulation results with the conventional DP model, the proposed approach yields much more accurate predictions on rapid water breakthrough and high water cut in fractured reservoirs.
1. Introduction
Effective modeling of fractured reservoir remains a challenge since realistic characterization of fractures is difficult, and highly recognized algorithms to derive field-scale fracture properties from seismic processing results are not yet available. Advances in seismic and image log interpretations for fracture mapping have been made for the last two decades, which make quantitative modeling of naturally fractured systems possible. As a result, it is viable, though with significant uncertainty, to represent actual fracture distribution in reservoir modeling and run detailed flow simulation via a discrete approach [1].
The discrete fracture model (DFM) was proposed to precisely characterize each fracture based on its actual geometry ([2–8]; Matthai et al., 2005). But it is widely recognized that DFMs are difficult to be applied to field-scale problems with complex geological structures and large number of fractures in multiple scales. Hybrid models that combine different scales of fractures and model them in different types and scales of grids offer an effective approach to solve this problem. Moinfar et al. [9] proposed an embedded discrete fracture model (EDFM) to model different scales of fractures all together into a dual-porosity type of media. The most significant advantage of EDFM is that it avoids unstructured mesh generation, which is typically very challenging given complex fracture distribution and interactions. However, in the original EDFM methodology, the transmissibilities between matrix cells and fractures are computed analytically assuming a linear pressure distribution along the flow path [10, 11]. Using DP models for the dense microfractures and DFM models for the sparse macroscale fractures, respectively, Wang [9] proposed a hybrid DP-DFM model to simulate SRV regions resulted from hydraulic fracturing in unconventional reservoirs. As an improvement, the transmissibilities between DP and DFM media are calculated according to their grid geometry to some extent. However, the shape factor computed in this work is based on sugar-cube distributed fracture assumptions and with specific boundary conditions, thus limits its application to more realistic and arbitrarily distributed fracture systems.
In this work, we propose a systematic approach to characterize the actual fracture distribution and, based on which, to generate a hybrid grid for discrete fracture model construction, and furthermore to perform flow simulation. In this hybrid model, the large-scale fractures are modeled as embedded discrete fractures and the small-scale fractures are modeled as dual-porosity fractures, respectively. And the transmissibilities between matrix and small-scale fractures, matrix and large-scale fractures, and different scales of fractures can be accurately computed from their actual geometries. The final hybrid model could be simulated through connection-list-based flow simulators.
This paper has proceeded as follows: Section 2 presents the methodology and workflow of the proposed approach. Section 3 demonstrates the application of the new modeling approach to an actual fractured reservoir in the South China Sea. Flow simulation results are then presented to demonstrate the advantages of the hybrid model over existing models.
2. Methodology and Workflow
The workflow to build the proposed hybrid model by applying different numerical models to multiscale fractures is shown in Figure 1. Firstly, the small-scale fracture network is generated following a standard discrete fracture network (DFN) modeling process based on image log interpretations. Secondly, a traditional DP grid system is constructed as a base grid model. Thirdly, the large-scale fractures are identified through a fracture reconstruction procedure based on poststack seismic interpretations. Fourthly, the discrete large-scale fractures are embedded into the base DP grid model, and the transmissibilities between the discrete large-scale fractures and both matrix and fracture in the base DP grid system are calculated. Finally, the hybrid model can be simulated in a connection-list-based flow simulator.
2.1. Large-Scale Fracture Reconstruction
Previous studies showed that seismic properties such as amplitude and velocity vary ellipsoidally with respect to azimuth in anisotropic mediums [12, 13]. So azimuthal analysis of P-wave seismic attributes offers an effective way for fracture predictions by extracting P-wave seismic attributes at different azimuths and then fitting these attributes to an ellipse (Figure 2). In this analysis, the ellipticity indicates the density or strength of the fracture, and the directions of major and minor axes indicate the orientation of the fracture.
Given the fracture density and azimuth interpreted from seismic anisotropy analysis, it is viable to reconstruct a discrete fracture distribution explicitly from the seismic grid [14]. Firstly, the fracture density and azimuth interpretations are mapped from seismic grid onto stratigraphic reservoir grid to form a fracture element network. Then, in the reservoir grid, by merging adjacent fracture elements that are similar in directions (azimuth) and density, the discrete fractures are reconstructed. The reconstruction results may vary due to different searching directions and routes applied. But it only introduces slight translation smaller than the length of grid cell, since the size of the background stratigraphic grid is significantly smaller compared to the fracture length (Figure 3). Similarly, in vertical directions, if the distance between the two fracture elements is small enough, and their dip angles are adequately consistent, we can merge the two fracture elements from neighboring layers (Figure 4).
2.2. Transmissibility Calculation for the Hybrid DP+EDFM Model
Using the connection-list-based simulator, the hybrid EDFM+DP model is specified by a connection list of the control volumes. This requires the input of cell volume, porosity, permeability, saturation, and depth information for each of the grid cells. Meanwhile, the two-point flux approximation is applied to correlate transmissibility with pressure difference and flow rate as in
Here, denotes the phase, and denote adjacent cells, is the flow rate, is the density of phase , and are the pressures of cell and , and is the transmissibility between cell and depending on the permeability and geometry of the control volumes.
The standard dual-porosity model assumes that the fractured media is considered as the overlapping continuum of a matrix continuum with well-defined porosity and permeability, in which they are locally connected to each other. The governing equations can be written for each continuum based on Darcy’s law and mass conservation: where the subscript and denote matrix and fracture, for the permeability of the subsurface formation, for fluid pressure, and and for fluid density and viscosity, respectively.
In EDFM, as the fractures are defined with individual cell quantities, the overall grid system turns into unstructured, although the background reservoir grid is typically a structured corner-point grid system. Thus, the key is to calculate transmissibility between nonneighboring matrix and fracture cells as in
where is the area open to flow, is the harmonic average of permeability, and is the characteristic distance between two control volumes associated with a nonneighbor connection (NNC).
We need to properly calculate transmissibilities for the three types of NNCs in Figure 5. The transmissibility for the three types of NNCs is calculated as follows: (1)For a pair of nonneighboring connected matrix and fracture cells, is denoted as the fracture surface area in the grid block, and is denoted as the harmonic average of the fracture and matrix permeabilities and is close to the lower-end matrix permeability typically. In addition, Li and Lee [10] and Hajibeygi et al. [11] assumed that in a grid block, the pressure variation is linear in the normal direction of each fracture surface. The following equation is used to calculate the average normal distance, :where is the volume element for integral, is the normal distance of the element to the fracture, and is the volume of a grid block (2)For a pair of nonneighboring connection between two intersecting fracture cells, we extend the approach proposed by Karimi-Fard et al. [5] to obtain the transmissibilities as inwhere is the length of the intersection line illustrated as the black solid line in Figure 5(b), and are fracture aperture and permeability, and is the average distance between the center of the fracture subsurface and the intersection line. In addition, for scenarios with more than two fractures intersect in one grid block, each pair of intersecting fracture volumes form a NNC and their transmissibility needs to be calculated. Finally, when two fractures penetrate the same grid block but they do not intersect with each other in the grid block, there is no NNC of the fractures in that grid block (3)For a pair of nonneighboring connection between two individual fracture cells, is equal to the fracture permeability, and is the distance between the centers of two fracture surfaces. The intersection line of the fracture surface and the common face of two neighboring grid blocks is shown in Figure 5(c), and is the intersection area calculated as the product of intersection line length and fracture aperture
(a)
(b)
(c)
3. Field Application
3.1. Background
The target reservoir is an offshore carbonate formation in the South China Sea that has been in production since 1996. It has complex geology and developed faults and fractures. The reservoir depth ranges from 1,198 to 1,273 m, its average porosity is 21.49%, and its average permeability is 363.14 mD. The reserves are evaluated to be 190.8 million sm3 with heavy oil (46.5-162.1 mPa·s). Reservoir pressure maintains well due to the large aquifer, and it is mainly developed by horizontal wells. The total production declined rapidly, which is a result of the rich fractures as they serve as water channels. The water cut rapidly increased above 90% by the end of 2001.
3.2. Hybrid Modeling of Multiscale Fractures
In this work, we have investigated different seismic attributes and picked the instantaneous azimuthal amplitude from P-wave reflection signals to map the fracture density and orientation. Figure 6 shows the fracture density map interpreted from poststack seismic analysis. This is usually a reflection of large-scale fractures.
Mapped fracture elements are shown in Figure 7. The azimuth of each fracture element is obtained through a local fitting process. Calibrated with image logs, the fracture density is then used to yield fracture apertures based on which to calculate fracture volumes and permeabilities.
Following the proposed methodology of fracture reconstruction, the large-scale fracture elements are merged laterally and vertically, and the discrete large-scale fracture network of the entire reservoir is reconstructed (Figure 8).
Through a regular DFN modeling process, the small-scale fractures are characterized by image log interpretations as shown in Figure 9. Small-scale fractures are modeled by the DP method, based on which the reconstructed large-scale fractures are embedded into following the proposed methodology. Finally, we generated a hybrid grid system containing 701,125 cells for the matrix and fracture cells, and the complete transmissibility and connection list is obtained for flow simulation.
3.3. Simulation Results and Discussion
To highlight the effectiveness of this approach, flow simulation results are presented together with the actual well history data, as well as the previous simulation results based on the standard DP model. Four horizontal wells with a relative long production history are of the greatest interest. And it is notable that this study only involves a few field-level history matching adjustments.
The simulation results are shown in the following, including water cut for the entire field (Figure 10), and oil rates (Figure 11, control condition), water cut (Figure 12), and water rates (Figure 13) of the four major wells. In terms of the overall water cut, results from our proposed hybrid model have a very good match with the historical data as shown in Figure 12. For the early stage of rapid water breakthrough, the DP+EDFM model simulates the trend much better than the DP model, due to the detailly represented discrete fractures. In terms of simulation time cost, the DP+EDFM model slightly increases to 11,162 s, which is 8,640 s for the DP model running with Eclipse. It is reasonable since the discrete large-scale fractures cause more nonlinearity.
Comparison of single well water cut results are similar. Our proposed hybrid model is able to capture the drastic water breakthrough period for four major wells; in contrast, the DP model fails to match the high water cut and rapid climbing for wells 1, 3, and 4. Other than well 2, all wells are intersected with large-scale fractures, resulting in a poor match for the DP model. In the DP model, water production rates are relatively low and changing more uniformly. While in our proposed model, both the trend and the level of water production, especially for the transient periods, are accurately matched.
The effects of different types of fractures can be seen between well 2 and other 3 wells. Due to their strong transmissibilities of large-scale fractures, large amount of water is transported from the aquifer into wells in a short period and rapid oil rate declines are seen in wells 1, 3, and 4. For well 2, the oil rate declined less rapidly from the peak rate. Bottom water invades more evenly through small-scale fractures compared to the large-scale ones. Both processes can be perfectly captured by the hybrid EDFM+DP model.
4. Conclusions
An efficient hybrid fractured reservoir modeling method is proposed and successfully applied to an actual field study. (1)A hybrid grid system is generated to precisely and efficiently model both large-scale and small-scale fractures through an integrated approach. It offers a convenient way of discretely characterizing the key fractures under the background structured grid(2)According to the field case study results, the hybrid EDFM+DP model can help achieve highly accurate simulation results compared to conventional DP models. The DP model fails to reproduce the rapid water breakthrough in early period caused by large-scale fractures. In terms of the computational cost, the proposed approach is slightly higher compared to conventional DP models and is significantly lower than a full DFM approach which requires unstructured gridding and introduces significantly more connections. It achieves a good balance between accuracy and efficiency in modeling and simulating naturally fractured reservoirs(3)To our best knowledge, the proposed approach is stable and robust when the large-scale fractures become even more complex. The hybrid grid generation process can be done in a few minutes with a fracture set of over 5,000 elements. Although there are no limits of fracture number and complexity on the grid generation part, it might be very challenging for the flow simulation part. Thus, it is recommended to keep the discrete fractures at a proper level for larger-scale field applications
Data Availability
Because the oilfield involved in the paper is still under development and the oilfield data is still in a confidential state, the oilfield data used in the paper cannot be made public now.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge the project named “Application of an Accurate and Efficient Modeling Approach to a Multiscale Fractured Reservoir in South China Sea” which is provided by the Science and Technology Major Projects of CNOOC China Limited (No. CNOOC-KJ 135 ZDXM 37 SZ 04 SZ).