## Simulation Study of CO2-EOR in Tight Oil Reservoirs with Complex Fracture Geometries

#### Pavel Zuloaga-Molero^{1}, Wei Yu^{2}, Yifei Xu^{1}, Kamy Sepehrnoori^{1} & Baozhen Li^{3}

^{1}Department of Petroleum and Geosystems Engineering, University of Texas at Austin, Austin, TX, 78712, USA. ^{2}Department of Petroleum Engineering, Texas A&M University, College Station, TX, 77843, USA. ^{3}CNOOC Research Institute, Beijing, 100027, China.

## Abstract

The recent development of tight oil reservoirs has led to an increase in oil production in the past several years due to the progress in horizontal drilling and hydraulic fracturing. However, the expected oil recovery factor from these reservoirs is still very low. CO_{2}-based enhanced oil recovery is a suitable solution to improve the recovery. One challenge of the estimation of the recovery is to properly model complex hydraulic fracture geometries which are often assumed to be planar due to the limitation of local grid refinement approach. More flexible methods like the use of unstructured grids can significantly increase the computational demand. In this study, we introduce an efficient methodology of the embedded discrete fracture model to explicitly model complex fracture geometries. We build a compositional reservoir model to investigate the effects of complex fracture geometries on performance of CO_{2} Huff-n-Puff and CO_{2} continuous injection.

The results confirm that the appropriate modelling of the fracture geometry plays a critical role in the estimation of the incremental oil recovery. This study also provides new insights into the understanding of the impacts of CO_{2} molecular diffusion, reservoir permeability, and natural fractures on the performance of CO_{2}-EOR processes in tight oil reservoirs.

# Introduction

Tight oil reservoirs are typically characterized by low porosity (<10%) and low permeability (<0.1 mD)1. The successful economic development of tight oil reservoirs in recent years hinges on two advanced technologies: horizontal drilling and multi-stage hydraulic fracturing. U.S. Energy Information Administration2 reported that tight oil production will increase from 33% of total lower 48 onshore oil production to 51% in 2040. However, the decline curves of primary production are steep due to low permeability.

The Bakken Field located in North Dakota, is one the most productive tight oil reservoirs in North America. The Bakken formation is composed of shale and dolomite layers. The Middle Bakken and Three Forks layers are the most important because of the reservoir quality and hydrocarbon saturation. The primary oil recovery factor in Bakken Formation is very low and typically less than 10%, resulting in tremendous oil resource remaining in place. Hawthorne et al.^{3} presented that a minor improvement in oil recovery factor such as 1% could yield several billion barrels of additional oil. Hence, it is important to investigate enhanced oil recovery (EOR) methods to increase oil recovery in tight oil reservoirs.

It is challenging to apply water flooding in tight oil reservoir due to low injectivity and poor sweep efficiency. Also, the oil-wet nature of some tight oil formations such as Bakken minimizes the effectiveness of water flooding. Both recent experimental and simulation studies have shown that carbon dioxide (CO2) injection could be a feasible EOR method to improve the oil recovery and carbon storage and sequestration in tight oil reservoirs^{4–9}.

CO_{2} has a considerably lower minimum miscibility pressure (MMP) than other gases such as N_{2} and CH_{4}^{10,11}. CO_{2}-EOR has two common operation scenarios: continuous CO_{2} injection, which is referred to as CO_{2} flooding in this study, and CO_{2} Huff-n-Puff. Although CO_{2}-EOR in conventional reservoir is well understood, it is relatively a new concept in tight oil reservoirs^{1}. Hawthorne *et al*.^{3} proposed five conceptual steps for CO_{2}-EOR process in tight oil formation: “(1) CO_{2} flows into and through the fractures, (2) unfractured rock matrix is exposed to CO_{2} at fracture surfaces, (3) CO_{2} permeates the rock driven by pressure, carrying some hydrocarbon inward;

however, the oil is also swelling and extruding some oil out of the pores, (4) oil migrates to the bulk CO_{2} in the fractures via swelling and reduced viscosity, and (5) as the CO_{2} pressure gradient gets smaller, oil production is slowly driven by concentration gradient diffusion from pores into the bulk CO_{2} in the fractures”. The CO_{2} molecular diffusivity is a key physical mechanism for CO_{2}-EOR process in tight oil reservoirs, which must be taken into account correctly when building a numerical compositional model^{12,13}.

An accurate modeling CO_{2}-EOR process in tight oil reservoirs is challenging. The actual hydraulic fracturing treatment often creates complex fracture networks by opening and interconnecting the pre-existing natural fractures. The fracture diagnostic technologies like microseismic^{14,15,16} and the recent developed fracture propagation models^{17,18,19} indicate that the fracture geometry is complex and non-planar. Although many efforts in the literature are dedicated to modeling CO_{2}-EOR in tight oil reservoirs^{20,21,22,23}, they made an assumption of simple bi-wing planar fracture geometry without considering the more-realistic complex fracture geometry. Even though the complex fracture networks provide the primary conduits for oil production, they might be detrimental to CO_{2} flooding effectiveness^{24}. This is because of early CO_{2} breakthrough and poor sweep efficiency during CO_{2} flooding. CO_{2} Huff-n-Puff operation might overcome this issue. Hence, more studies are needed to understand this operation in detail. Based on our knowledge, there are no published studies to date that focused on modeling complex fracture networks and natural fractures explicitly and investigating their effects on the CO_{2}-EOR effectiveness in multiple horizontal wells.

Numerical modeling of the well performance from tight oil reservoirs based on explicitly handling the complex fracture networks remains a challenging topic. Although significant attempts have been focused on developing unstructured grids method to model the complex fracture networks^{25,26}, it is still limited in field-scale application due to its complexity in gridding and large computational demand, especially in the simulation of CO_{2}-EOR process using compositional numerical models with multiple components. Accordingly, an efficient approach to simulate CO_{2}-EOR process from tight oil reservoir with complex fracture geometry is still lacking in the petroleum industry.

In this study, we introduced an embedded discrete fracture model (EDFM), which was originally proposed by Li and Lee^{27}. Further extensions have been done by Moinfar *et al*.^{28} and Cavalcante Filho *et al*.^{29}. Based on the EDFM, we can modify the compositional reservoir simulator in a non-intrusive manner to accurately and efficiently handle the complex fracture geometry^{30}. We verified the EDFM methodology for simple hydraulic fractures against the traditional local grid refinement (LGR) approach^{31}.

In this study, we first built a compositional numerical model, including four hydraulic fractures within a single stage based on the actual fluid and reservoir properties of the Bakken Formation. Four case studies with different fracture geometries from simple to complex were modelled using the EDFM. In addition, we evaluated the impacts of different fracture geometries on the well performance of CO_{2} Huff-n-Puff after three years of primary production. Subsequently, we extended the model to include two wells, each well has four fractures within single stage, to simulate CO_{2} flooding. Finally, we built a field-scale reservoir model, including two horizontal wells to study the effects of complex non-planar hydraulic fractures and natural fractures on the well performance of both CO_{2} injection scenarios under two different reservoir permeabilities of 0.01 md and 0.1 md. This work provides critical insights about the effect of fracture complexity on the well performance of CO_{2} Huff-n-Puff and CO_{2} flooding in tight oil reservoirs.

## Results

### 2

We performed four case studies for CO2 Huff-n-Puff simulation by considering different fracture geometries within a single stage including four hydraulic fractures.

- Case 1: Planar fractures (Note that case 1 is considered as reference case)
- Case 2: Diagonal fractures
- Case 3: Reoriented fractures
- Case 4: Fracture networks

These cases consider different fracture geometries with different degrees of complexity, such as the striking angle between horizontal well and fractures, irregular fracture length of individual fracture, and the creation of fracture networks. These cases were used to quantify the effect of the increasing fracture complexity on the well performance for primary production and Huff-n-Puff scenario.

#### Case 1: Planar fractures.

This case presents four simple planar fractures orthogonal to the horizontal wellbore. All fractures have the same properties (fracture half length: 210 ft, fracture height: 40 ft and fracture conductivity: 50 md-ft).

Figure 1(a) shows the projection of the aforementioned fractures. Figure 2(a) shows the comparison of oil recovery factor between the primary production and Huff-n-Puff scenario. It can be seen that the difference of the oil recovery factor at the end of simulation period is around 3.9%. As a reference case, this model was used to evaluate the impact of various cases on the well performance.

**Figure 1.** Illustration of fractures geometries evaluated for the Huff-n-Puff scenario. (a) Case 1: Planar fractures. (b) Case 2: Diagonal fractures. (c) Case 3: Reoriented fractures. (d) Case 4: Fracture networks. The yellow planes represent the fractures and the red line represents horizontal wellbore.

#### Case 2: Diagonal fractures.

The first degree of complexity considered in this study was the orientation of the fractures with respect to the horizontal wellbore. In case 2, four fractures from an angle of 45° with the wellbore, as shown in Fig. 1(b). The other fracture properties are kept the same as case 1. Although the distance between perforation clusters is the same as case 1, the orthogonal fracture distance between two neighboring fractures was reduced to from 140 ft to 99 ft due to the change of fracture orientation.

Figure 2(b) compares the difference of the oil recovery factor between cases 1 and 2. As shown, the primary oil recovery without CO_{2} injection is almost the same. However, the incremental oil recovery factor of case 2, which is around 2.9%, is less than case 1 with 3.9%. This difference is due to the changes in CO_{2} molecule distribution around the fractures. There is a lower production interference between the fractures in case 1 than in case 2, resulting in a less effective CO_{2} Huff-n-Puff for case 2. The increase in the interference between the fractures of case 2 is expected since the orthogonal distance between the fractures decreases.

**Figure 2.** Comparison of the oil recovery factors. (a) Oil recovery factor for primary production and Huff-n- Puff scenario for planar fractures. (b) Comparison between Case 1 and Case 2. (c) Comparison between Case 1 and Case 3. (d) Comparison between Case 1 and Case 4.

### Case 3: Reoriented fractures.

Case 3 represents four non-planar fractures with outer fractures longer than inner fractures, as shown in Fig. 1(c). Each fracture is composed of two fracture segments with an orientation of 45° and 135° or 135° and 45° for the upper and lower segment, respectively. The fracture half-length for outer and inner fractures is 295 ft and 125 ft, respectively. It should be noted that the total fracture length, i.e. the summation of the fracture segments, remains the same as cases 1 and 2. Figure 2(c) compares the difference in the oil recovery factor between cases 1 and 3. As shown, there is a small difference of primary production. Nevertheless, when compared to case 1, the incremental oil recovery factor of case 3 increases from 3.9% to around 4.5%. Again, this improvement is associated with the interference between the fractures. In this case, the Huff-n-Puff is more effective than case 1 because the outer fractures have a longer length and less interference with the inner fractures.

### Case 4: Fracture networks.

The last case of fractures analyzed represents four systems of fracture networks created around each perforation cluster, as shown in Fig. 1(d).

The fracture networks are composed of several segments that intersect each other. Similar to previous cases, the total fracture length of all segments is equivalent to the planar fractures of case 1. Case 4 represents a more realistic fracture geometry. Figure 2(d) compares the difference of the oil recovery factor between cases 1 and 4.

_{2}_{2}_{2}_{2}

**Table 1.** Summary of oil recovery factor (RF) of primary production and incremental RF after CO2 Huff-n-Puff and flooding for the different cases and CO2-contacted area for each case.

In addition, Table 1 lists the oil recovery factor (RF) for primary production and the incremental RF after Huff-n-Puff for each case. The results show small differences of primary production, but significant changes after Huff-n-Puff for different fracture geometries. We can notice that there is a consistent reduction in the incremental RF of CO_{2} Huff-n-Puff as the fractures become shorter and closer to each other, and the CO_{2}-contacted area decreases as a result of the fracture interference. Case 4, which has the highest fracture complexity, has the lowest incremental oil recovery. This shows that characterizing the actual fracture geometry and accurately modeling the well performance from this geometry play an important role in the estimation of the additional oil recovery after the Huff-n-Puff stimulation. The use of simple planar fractures to simplify the complex fractures might overestimate the CO_{2}-EOR effectiveness.

#### Case Studies of CO_{2} Flooding

For the evaluation of CO_{2} flooding scenario, the same single-stage fracture geometries as CO_{2} Huff-n-Puff scenario were considered, as shown in Fig. 3. For these case studies, the simulation model was extended to include two horizontal wells with identical fracture geometries. The reservoir and fracture properties remain the same as those mentioned in the previous section. The distance between two wells is fixed at 1,020 ft for each case. The CO_{2} flooding scenario considers an initial period of primary production of 3 years for both wells.

After that, one of the producing wells is converted to injector and used for CO_{2} injection until the end of the production time (18 years in total). The cumulative CO_{2} injection is comparable to the amount of the Huff-n-Puff scenario. The primary production simulated without CO_{2} injection is used to measure the incremental oil recovery after CO_{2} flooding.

**Figure 3. **Illustration of fractures geometries evaluated for the flooding scenario. (a) Case 1: Planar fractures. (b) Case 2: Diagonal fractures. (c) Case 3: Reoriented fractures. (d) Case 4: Fracture networks. The yellow planes represent the fractures and the red line represents horizontal wellbores.

Figure 4 compares the oil recovery factor curves for these four cases. Table 1 summarizes the oil RF of primary production and incremental RF after CO_{2} flooding and CO_{2}-contacted area in each case at the end of production. As it can be seen, the incremental RF of cases 1 and 4 is lower than that of cases 2 and 3. It can be noticed from the results that there is not a direct relationship between the increase in fracture complexity and the CO_{2} flooding effectiveness. However, the incremental RFs are still in agreement with the area contacted by the CO_{2} injected, even though there are some small divergences as it can be observed in cases 2 and 3. Case 2 shows a higher recovery than case 3 but it has a slightly lower CO_{2} contacted area. Nonetheless, this can be expected since the average CO_{2} concentration in the contacted area in case 3 is higher than case 2 (0.60 and 0.55, respectively). For the flooding scenario a higher contacted area is due to the location of the fractures. More particularly, it is because of the dimensions of the cross sectional area covered by the fractures between the injector and producer. It can be clearly observed in the Fig. 3 that cases 2 and 3 have a larger extension in the x-axis direction (717 ft and 572 ft respectively), whereas cases 1 and 4 have a smaller extension (420 ft and 529 ft respectively). A longer extension of the fractures in the x-axis direction allows a higher CO_{2}-contacted area and therefore a higher incremental oil recovery. For the CO_{2} flooding, the fracture geometry is also a key factor affecting the estimation of the additional oil recovery. However, the inclusion of the complex fractures for the flooding scenario, unlike the Huff-n-Puff, does not necessarily implies a negative effect in the incremental oil recovery.

**Figure 4.** Comparison of oil recovery factor curves for four cases under CO_{2} flooding scenario.

### Field-Scale Case Study.

This section is devoted to analyze the influences of complex fracture geometries on the well performance of CO_{2} Huff-n-Puff and CO_{2} flooding by building a realistic field-scale reservoir model. Two permeabilities such as 0.01 md and 0.1 md are considered to represent low and high permeability, which are within the range of permeability in the Bakken Formation. In addition, we also examine the effect of natural fractures on the production performance. For these purposes, the original model was extended to 5,240 ft × 2,680 ft × 40 ft (262 × 134 × 1 cells), which is able to model two horizontal wells with lateral length of 4,640 ft for each one and well spacing of 1,340 ft.

Each well has 15 stages and each stage is assumed to have a single effective hydraulic fracture. The fracture spacing is 290 ft. We performed both CO_{2} Huff-n-Puff and CO2 Flooding simulations by considering the following cases:

- Non-planar fractures
- Non-planar fractures with one set of natural fractures

In the first case, shown in Fig. 5(a), all the fractures are non-planar and have different dimensions. The fracture half-length varies from 194 to 445 ft.The fracture height and conductivity used are 40 ft and 50 md-ft, respectively. Both fracture geometries have the same fracture spacing of 290 ft.

One additional cases were set in order to evaluate the impact of the natural fractures on the performance of the CO_{2} Huff-n-Puff and CO_{2} flooding, as presented in Fig. 5(b). This case includes a set of 300 natural fractures, which are randomly distributed with the assumption that their orientations are parallel to the horizontal wellbore. These natural fractures have a length ranging from 100 to 200 ft, a height of 40 ft, and a fracture conductivity of 5 md-ft. The combination of modelling complex non-planar hydraulic fractures and natural fractures permits to model more realistic fracture networks.

**Figure 5. ** Illustration of non-planar fracture geometry used for the field-scale simulation, taking into account the presence of natural fractures. (a) Non-planar fractures with irregular dimensions (b) Non-planar fractures with a set of natural fractures parallel to the horizontal wellbore.

_{2}

We performed the simulations for CO_{2} Huff-n-Puff scenario and the results for both high and low permeabilities are shown in Fig. 6. For the cases evaluated there was no significant differences in the oil recovery factor during primary production regardless of high or low permeability. For the high permeability of 0.1 md, the maximum difference between the cases is less than 0.05%, and for the low permeability of 0.01 md the difference is about 0.4%. Similarly, the enhanced oil recovery at the high permeability, shown in Fig. 6(a), reflects a small impact. The incremental oil recovery factor is about 7.7% for the non-planar fractures when compared to the primary production. The previous results observed in the single stage study showed a bigger difference because the fracture spacing used was 140 ft, which is lower than the fracture spacing of 290 ft used for the field case. Then, the fracture interference decreases with increasing fracture spacing and there is only a small difference of 0.6% in the incremental oil recovery factor between the planar and non-planar fracture geometries. In addition, the effect of natural fractures is also small. The incremental oil recovery factor is 8.2% for the additional set of natural fractures.

On the other hand, the results at low permeability of 0.01 md, presented in Fig. 6(b), show a significant difference in the incremental oil recovery factor when the natural fractures are included. The incremental recovery factor in the case with one set of natural fracture increases from 5.2% to 6.8% when compared to the non-planar fractures without natural fractures. Hence, the presence of natural fractures significantly impacts the well performance of CO_{2} Huff-n-Puff under the low permeability of 0.01 md. It can be implied that the natural fractures should be correctly included in the numerical model for simulating CO_{2} Huff-n-Puff in some tight formations with low permeability and high density of natural fractures.

**Figure 6.** Oil recovery factors for CO2 Huff-n-Puff and CO2 flooding under matrix permeability values of 0.1 md and 0.01 md. (a) Oil recovery factors of CO2 Huff-n-Puff for permeability of 0.1 md. (b) Oil recovery factors of CO2 Huff-n-Puff for permeability of 0.01 md. (c) Oil recovery factors of CO2 flooding for permeability of 0.1 md. (d) Oil recovery factors of CO2 flooding for permeability of 0.01 md.

The previous observations can be verified by comparing CO_{2} molecule distribution maps, as shown in Fig. 7. The molecular fraction distribution maps of CO_{2} presented in Fig. 7(a,b) show the concentrations of CO_{2} for the non-planar fractures with and without the set of natural fractures under low permeability of 0.01 md. It can be clearly noticed that under low permeability, the non-planar fracture geometry with one set of natural fractures has a higher concentration of CO_{2} compared to without natural fractures especially for the values of around 0.6 displayed in cyan. Hence, the CO_{2} distribution maps show that the performance of CO_{2} Huff-n-Puff under the low permeability is more sensitive to the presence of natural fractures.

**Figure 7.** Global CO_{2} molecule distribution for Huff-n-Puff and flooding scenarios with and without natural fractures under matrix permeability values of 0.1 md and 0.01 md. (a) Global CO_{2} molecule distribution for Huff-n-Puff after one year of CO_{2} injection under permeability of 0.01 md without natural fractures (b) Global CO_{2 } molecule distribution for Huff-n-Puff after one year of CO_{2} injection under permeability of 0.01 md with natural fractures (c) Global CO_{2} molecule distribution for flooding after 15 years of CO_{2} injection under permeability of 0.1 md. (d) Global CO_{2} molecule distribution for flooding after 15 years of CO_{2} injection under permeability of 0.01 md.

### CO_{2} flooding for field-scale study.

Similarly, we performed the simulations for CO_{2} flooding scenario and the results for both high and low permeabilities are shown in Fig. 6. The effect of the natural fractures, which are parallel to the horizontal wellbore, is not significant and the incremental oil recovery factor is similar to that obtained from the non-planar fractures without natural fractures (21.5% and 21.2%, respectively). Regarding the comparison for the low permeability of 0.01 md and high permeability of 0.01 md, Fig. 7(d) show that the flooding scenario for low permeability is not favorable because the injected CO_{2} cannot reach the production well (Well 1) due to the poor connectivity between two wells, unlike the case of high permeability shown in Fig. 7(c). The final recovery factors are lower than the primary production, which includes the production of both wells. Hence, we can notice that the presence of natural fractures has a higher importance in the production response of the CO_{2}-EOR in the cases with low permeability, especially for the Huff-n-Puff scenario.

## Discussion

The reservoir model built allowed to evaluate the effects of the complex fracture geometries on the well performance of CO_{2}-EOR using the EDFM approach. From the presented results we observed the fracture geometry has an important effect on the Huff-n-Puff treatment. The CO_{2} Huff-n-Puff is quite sensitive to the fracture geometry because its effectiveness is mainly determined by the fracture interference.

As we have mentioned previously, the interference between fractures increases as the fracture segments become closer and shorter, which are the cases of complex fractures and fracture networks. Therefore, the expected incremental recovery decreases. Most of the models reported in the scientific literature use planar fractures, but if this model is used for the evaluation of CO_{2} Huff-n-Puff, the incremental recovery will be certainly overestimated.

For the CO_{2} flooding scenario the non-planar fractures do not show a negative effect on the EOR performance. It can improve the effectiveness of the flooding if the location of the fractures turns into a higher contacted area, which depends on the configuration of the fractures. The CO_{2} flooding is less sensitive to the fracture geometry itself, but it is affected by the extension of the injection front developed by such fracture geometry.

Although there has been an improvement in the techniques of fracture diagnostics, the actual fracture geometry is still not well understood. The fracture geometry obtained from such techniques is non-unique and the uncertainty in the fracture geometry should be taken into account for the evaluation of EOR treatments. For an adequate estimation, it would be useful to make a sensitivity analysis of all the possible fracture geometries to obtain a range of values for the incremental recovery, especially for the CO_{2} Huff-n-Puff. For the CO_{2} flooding the focus should be on the extension of the cross-sectional area of the gas front.

The case studies performed with natural fractures present in the reservoir show that the improvement of the CO_{2} Huff-n-Puff is only significant for the low permeability of 0.01 md. This increment in the oil recovery responds to the increase in the CO_{2} contacted area. In the cases of high permeability the contacted area was already high and the natural fractures did not affect it significantly. However, The CO_{2} flooding is not favorable for low permeability of 0.01 md due to poor connectivity between two wells, while it is suitable for high permeability of 0.1 md. The CO_{2} Huff-n-Puff performs better than the CO_{2} flooding for low permeability of 0.01 md.

## Methods

### Embedded Discrete Fracture Model.

EDFM is an efficient model to handle complex fracture geometry in reservoir simulators. Using this method, fractures can be discretized into small fracture segments with matrix cell boundaries. Virtual cells are created to represent these fracture segments. This method can be applied in traditional reservoir simulators in a non-intrusive manner by appending cells in the grid domain and adding non-neighboring connections (NNCs) for these cells to account for the mass transport associated with fractures, including the flow between matrix and fractures, flow inside an individual fracture, and flow between intersecting fractures. Through transmissibility factors between corresponding cells, the volume flow rate of phase l between cells in a NNC pair is:

where λ_{l} is the relative mobility of phase l, TNNC is the NNC transmissibility factor, and ΔP is the potential difference between the cells. For the calculation of transmissibility factors, generally, TNNC can be expressed as

where *k*_{NNC}, *A*_{NNC}, and *d*_{NNC} are the permeability, contact area, and distance associated with this connection, respectively. For matrix-fracture connection, in Eq. 2, *k*_{NNC} is the matrix permeability in the direction perpendicular to the fracture plane, *A*_{NNC} is the area of the fracture plane inside the matrix block, *d*_{NNC} is the average normal distance from matrix block to fracture plane.

For connections between fracture segments, *k*_{NNC} is calculated as an average of fracture permeability, *A*_{NNC} is the common area between fracture segments, and *d*_{NNC} is the distance between centroids of the fracture segments. More details of the calculation can be found in Xu^{30}. In addition, the fracture-wellbore intersections are modeled as effective well indices in the EDFM. The effective well index in the EDFM can be calculated as^{30}

where *k*_{f} is the fracture permeability, *w*_{f} is the fracture aperture, *L* and *W* are the length and height of the fracture segment, respectively. Eqs 3 and 4 can be derived by replacing the dimensions and permeability of the gridblock in Peaceman’s model with those of the fracture segments.

By discretizing large fractures into interconnected small fractures, the EDFM has been proven to be effective in modeling complex hydraulic fracture geometries such as non-planar fractures and fractures with variable width^{30}. Furthermore, the computational performance of the EDFM in traditional simulators was also verified through detailed comparison with LGR model. In this study, we first present a case study to verify the EDFM with LGR model, then for other cases, we apply the EDFM method in the simulator to simulate complex fracture geometries.

### Reservoir Simulation Model

A compositional numerical reservoir model with single fracture stage was built based on the typical reservoir and fracture properties of the Middle Bakken Formation using a reservoir simulator^{32}. A Cartesian grid system was used, which consists of 70 grids in *x* direction, 51 grids in *y* direction, and 1 grid in *z* direction. The dimension of the reservoir model is 1,400 ft × 1,020 ft × 40 ft, which corresponds to length, width, and thickness, respectively. A porosity of 7% and a permeability of 0.1 md were considered which are the typical reservoir and fracture properties from the Middle Bakken Formation^{33,34}. The study considers homogeneous and isotropic reservoir. The relative permeability curves including the water-oil relative permeability and liquid-gas relative permeability are from our previous study^{12}, which were generated based on history matching with a field production well from the Middle Bakken Formation. More detailed information can be found in the Supplementary Information document. The molecular diffusion is one the main mechanisms of EOR in tight reservoirs, since the gravity and viscous forces are not dominant because of the low permeability. The Sigmund method was used to model the physical mechanism of CO_{2} molecular diffusion^{35,36}. The CO_{2} molecular diffusion coefficient in the oil phase was set at 0.001 cm^{2}/s in this study. More details about the CO_{2} molecular diffusion effect can be found in the work by Yu *et al*.^{13}.

The oil composition used in the model and the main oil properties were taken from the study of crude oil composition for the Middle Bakken Formation by Yu *et al*.^{13}, which were determined based on matching the key fluid properties reported for the Bakken formation in the scientific literature^{4,33}. The Peng-Robinson equation of state^{37} was used and the model considers seven pseudo-components, i.e., CO_{2}, N_{2}, CH_{4}, C_{2}-C_{4}, C_{5}-C_{7}, C_{8}-C_{9}, C_{10+}, with the corresponding molar fractions of 0.02%, 0.04%, 25%, 22%, 20%, 13%, and 19.94%, respectively. The phase behavior of the crude sample was modelled using CMG-WinProp^{38}. The final determined oil properties are: oil gravity is 42° API, GOR is 1,000 SCF/bbl, bubble point is 2,000 psi, oil formation factor is 1.6 STB/bbl. It is important to point out that these oil properties are within the reasonable range of typical oil properties of the Middle Bakken Formation. In addition, the MMP was calculated as 3,334 psi, which was also in the range reported by other experiments^{6}.

The initial conditions considered for the simulation were an initial pressure of 8,000 psi for a reservoir depth of 11,000 ft, an initial oil saturation of 80%, an initial water saturation of 20%, and a reservoir temperature of 240 °F. The initial pressure does reflect an overpressured reservoir, as it is usually observed in tight and shale oil fields^{39}. The overpressure is related to the cracking of the organic material and the increase in fluid volume. Also, it is an indicator of the maturity of the shale formation and correlated to a higher oil production.

The model uses a no-flow boundary condition. The calculated original oil in place was 3.746 MMSTB. The bottomhole pressure (BHP) was used as the primary constraint for the reservoir simulation and a value of 1,800 psi was set for the production wells. For the injection wells, two constraints were used. The first one is a constant injection rate. 1,000 MSCF/D and 2,500 MSCF/D were used for the flooding case and Huff-n-Puff case, respectively. Both cases have the same cumulative injection volume. The second one is a maximum BHP of 8,000 psi. The BHP constraint and injection rates were set in such a way that the MMP could be reached after a short time of injection and the pressure is maintained then above the MMP to achieve the miscible condition during the whole simulation period. If the miscible condition is not achieved, the CO_{2} Huff-n-Puff scenario will show a poor performance because the presence of multiphase flow results in a decrease in relative permeability of oil. The optimized production schemes reported in the literature performed also the CO_{2} injection under miscible conditions^{40,41}.

### Acknowledgements

We gratefully acknowledge Computer Modeling Group Ltd. for providing the CMG software for this study.

### Author Contributions

P.Z.M. and W.Y. created the simulation model, designed the cases and performed the simulation studies. Y.X. designed the preprocessor of EDFM for fracture modelling. All authors wrote the paper.

### Additional Information

- Competing financial interests: The authors declare no competing financial interests.
- doi: 10.1038/srep33445 (2016).
- Wei Yu (email: yuwei127@tamu.edu)

## References

- Sorensen, J. A. et al. Characterization and evaluation of the Bakken petroleum system for CO2 enhanced oil recovery. Presented at the Unconventional Resources Technology Conference, San Antonio, Texas, USA, July 20–22. SPE 178659-MS (2015).
- U.S. Energy Information Administration. 2013 Early Release Overview. (2013).
- Hawthorne, S. B. et al. Hydrocarbon mobilization mechanisms from upper, middle, and lower Bakken reservoir rocks exposed to CO2. Presented at the SPE Unconventional Resources Conference, Calgary, Canada, November 5–7. SPE 167200-MS (2013).
- Kurtoglu, B., Sorensen, J. A., Braunberger, J., Smith, S. & Kazemi, H. Geologic characterization of a Bakken reservoir for potential CO2 EOR. Presented at the Unconventional Resources Technology Conference, Denver, Colorado, USA, August 12–14. SPE 168915- MS (2013).
- Song, C. & Yang, D. Performance evaluation of CO₂ huff-n-puff processes in tight oil formations. Presented at the SPE Unconventional Resources Conference Canada, Calgary, Alberta, Canada, November 5–7. SPE 167217-MS (2013).
- Adekunle, O. O. & Hoffman, B. T. Minimum miscibility pressure studies in the Bakken. Presented at the SPE Improved Oil Recovery Symposium, Tulsa, Oklahoma, USA, April 12–16. SPE 169077-MS (2014).
- Dai, Z. et al. An integrated framework for optimizing CO2 sequestration and enhanced oil recovery. Environmental Science & Technology Letters 1, 49–54 (2014).
- Ren, B. et al. Performance evaluation and mechanisms study of near-miscible CO2 flooding in a tight oil reservoir of Ji Lin oilfield China. Journal of Natural Gas Science and Engineering 27, 1796–1805 (2015).
- Ren, B., Ren, S., Zhang, L., Chen, G. & Zhang, H. Monitoring on CO2 migration in a tight oil reservoir during CCS-EOR in Jilin oilfield China. Energy 98, 108–121 (2016).
- Stalkup, F. I. Displacement behavior of the condensing/vaporizing gas drive process. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, USA, September 27–30. SPE 16715-MS (1987).
- Holm, L. W. Evolution of the carbon dioxide flooding processes. Journal of Petroleum Technology 39, 1337–1342 (1987).
- Yu, W., Lashgari, H. & Sepehrnoori, K. Simulation study of CO₂ huff-n-puff process in Bakken tight oil reservoirs. Presented at the SPE Western North American and Rocky Mountain Joint Meeting, Denver, Colorado, USA, April 17–18. SPE 169575-MS (2014).
- Yu, W., Lashgari, H. R., Wu, K. & Sepehrnoori, K. CO₂ injection for enhanced oil recovery in Bakken tight oil reservoirs. Fuel 159, 354–363 (2015).
- Fisher, M. K. et al. Optimizing horizontal completion techniques in the Barnett shale using microseismic fracture mapping. Presented at the SPE Annual Technical Conference and Exhibition, Houston, Texas, USA, September 26–29. SPE 90051-MS (2004).
- Warpinski, N. R., Mayerhofer, M. J., Vincent, M. C., Cipolla, C. L. & Lolon, E. Stimulating unconventional reservoirs: maximizing network growth while optimizing fracture conductivity. Journal of Canadian Petroleum Technology 48, 39–51 (2009).
- Maxwell, S. C., Weng, X., Kresse, O. & Rutledge, J. Modeling microseismic hydraulic fracture deformation. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, USA, September 30–October 2. SPE 166312-MS (2013).
- Wu, R., Kresse, O., Weng, X., Cohen, C. & Gu, H. Modeling of interaction of hydraulic fractures in complex fracture networks. Presented at the SPE Hydraulic Fracturing Technology Conference, The Woodlands, Texas, USA, February 6–8. SPE 152052-MS (2012).
- Xu, G. & Wong, S. W. Interaction of multiple non-planar hydraulic fractures in horizontal wells. Presented at the International Petroleum Technology Conference, Beijing, China, March 26–28. IPTC-17043 (2013).
- Wu, K. & Olson, J. E. Simultaneous multi-frac treatments: fully coupled fluid flow and fracture mechanics for horizontal wells. SPE Journal 20, 337–346 (2015).
- Hoffman, B. T. & Shoaib, S. CO₂ flooding to increase recovery for unconventional liquids-rich reservoirs. Journal of Energy Resources Technology 136, 022801.1–022801.10. (2013).
- Wan, T., Yang, Y. & Sheng, J. J. Comparative study of enhanced oil recovery efficiency by CO₂ injection and CO₂ huff-n-puff in stimulated shale oil reservoirs. Presented at AICHE annual meeting, Atlanta, Georgia, USA, November 17 (2014).
- Sanchez Rivera, D. Reservoir simulation and optimization of CO₂ huff-and-puff operations in the Bakken shale. MS Thesis, University of Texas at Austin, Austin, Texas (2014).
- Alharthy, N. et al. Enhanced oil recovery in liquid-rich shale reservoirs: laboratory to field. Presented at the SPE Annual Technicaxl Conference and Exhibition, Houston, Texas, USA, September 28–30. SPE 175034-MS (2015).
- Jarrell, P. M., Fox, C. E. & Stein, M. H. Practical aspects of CO₂ flooding (SPE monograph series, Richardson, Texas, 2002).
- Hoteit, H. & Firoozabadi, A. Compositional modeling of discrete-fractured media without transfer dunctions by the discontinuous Galerkin and mixed methods. SPE Journal 11, 341–352 (2006).
- Hui, M. H., Mallison, B. T., Fyrozjaee, M. H. & Narr, W. The upscaling of discrete fracture models for faster, coarse-scale simulations of IOR and EOR processes for fractured reservoirs. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, USA. September 30-October2. SPE 166075-MS (2013).
- Li, L. & Lee, S. Efficient field-scale simulation of black oil in naturally fractured reservoir through discrete fracture networks and homogenized media. SPE Reservoir Evaluation & Engineering 11, 750–758 (2008).
- Moinfar, A., Varavei, A., Sepehrnoori, K. & Johns, R. T. Development of an efficient embedded discrete fracture model for 3D compositional reservoir simulation in fractured reservoirs. SPE Journal 19, 289–303 (2014).
- Cavalcante Filho, J. S. A., Shakiba, M., Moinfar, A. & Sepehrnoori, K. Implementation of a preprocessor for embedded discrete fracture modeling in an IMPEC compositional reservoir simulator. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, USA, February 23–25. SPE-173289-MS (2015).
- Xu, Y. Implementation and application of the embedded discrete fracture model (EDFM) for reservoir simulation in fractured reservoirs. MS Thesis, University of Texas at Austin, Austin, Texas (2015).
- Rubin, B. Accurate simulation of non Darcy flow in stimulated fractured shale reservoirs. Presented at the SPE Western Regional Meeting, Anaheim, California, USA, May 27–29. SPE-132093-MS (2010).
- CMG-GEM. GEM User’s Guide, Computer Modeling Group Ltd (2015).
- Nojabaei, B., Johns, R. T. & Chu, L. Effect of capillary pressure on phase behavior in tight rocks and shales. SPE Reservoir Evaluation & Engineering 16, 281–289 (2013).
- Pilcher, R. S., Ciosek, J. M. McArthur, K., Hohman, J. & Schmitz, P. J. Ranking production potential based on key geological drivers-Bakken case study. International Petroleum Technology Conference, Bangkok, Thailand, Febreuary 7–9 (2012).
- Sigmund, P. M. Prediction of molecular diffusion at reservoir conditions. Part II - estimating the effects of molecular diffusion and convective mixing in multicomponent systems. Journal of Canadian Petroleum Technology 15, 53–62 (1976).
- Sigmund, P. M., Kerr, W. & MacPherson, R. E. A laboratory and computer model evaluation of immiscible CO2 flooding in a low-temperature reservoir. Presented at the SPE Enhanced Oil Recovery Symposium, Tulsa, Oklahoma, USA, April 15–18. SPE-12703- MS (1984).
- Peng, D. Y. & Robinson, D. B. A new two-constant equation of state. Industrial & Engineering Chemistry Fundamentals 15, 59–64 (1976).
- CMG-WinProp. WinProp User’s Guide, Computer Modeling Group Ltd (2015).
- Meissner, F. F. Petroleum geology of the Bakken formation, Williston Basin, North Dakota and Montana. The economic geology of the Williston Basin: Proceedings of the Montana Geological Society, 24th Annual Conference, 207–227 (1978).
- Ampomah, W. et al. Farnsworth field CO2-EOR project: performance case history. Presented at the SPE Improved Oil Recovery Conference. Tulsa, Oklahoma, USA. April 11–13. SPE-179528-MS (2016).
- Gunda, D., Ampomah, W., Grigg, R. B. & Balch, R. S. Reservoir fluid characterization for miscible enhanced oil recovery. Presented at the Carbon Management Technology Conference, Sugarland, Texas USA, November 16–19 (2015).

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/