Numerical simulation of CO2 huff-n-puff in complex fracture networks of unconventional liquid reservoirs

Accepted Manuscript Numerical Simulation of CO2 Huff-n-puff in Complex Fracture Networks of Unconventional Liquid Reservoirs Jianlei Sun, Amy Zou, Edith Sotelo, David Schechter PII:

S1875-5100(16)30138-X

DOI:

10.1016/j.jngse.2016.03.032

Reference:

JNGSE 1344

To appear in:

Journal of Natural Gas Science and Engineering

Received Date: 30 November 2015 Revised Date:

15 February 2016

Accepted Date: 9 March 2016

Please cite this article as: Sun, J., Zou, A., Sotelo, E., Schechter, D., Numerical Simulation of CO2 Huff-n-puff in Complex Fracture Networks of Unconventional Liquid Reservoirs, Journal of Natural Gas Science & Engineering (2016), doi: 10.1016/j.jngse.2016.03.032. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Numerical Simulation of CO2 Huff-n-puff in Complex Fracture Networks of Unconventional Liquid Reservoirs Jianlei Sun a*, Amy Zou a, Edith Sotelo a, David Schechter a a

Harold Vance Department of Petroleum Engineering, Texas A&M University, College Station, TX 77840, USA

RI PT

* Corresponding author, [email protected]

ABSTRACT

M AN U

SC

Horizontal drilling and hydraulic fracturing produce Unconventional Liquid Reservoirs (ULRs) economically. One of the characteristics of the ULRs is high production declining rate, and thus enhanced oil recovery techniques such as Huff-n-Puff will ensue. In addition, the existence of complex fracture network has been proven by both field and simulation work. This study investigated the potential of CO2 huff-n-puff in ULRs with complex fracture networks. Two fracture characterization techniques for complex fracture network such as fractal-based, and microseismic-based were applied to reduce uncertainties of fracture networks and thus improve the accuracy of simulation input. In order to study complex CO2 recovery mechanism, lab-measured data was upscaled to the field-scale, and then used for evaluating production performance of a simple fracture network. Detailed sensitivity analysis showed that diffusion in the field scale is negligible compared with other recovery mechanisms such as oil swelling, gas expansion, viscosity reduction, etc. The following comparison between dual continuum and the explicit fracture modeling approaches highlights the importance of heterogeneity due to the explicit fractures. Dual continuum models yields less accurate pressure behavior than explicit discrete fracture models even though production behavior could be history matched.

1. Introduction

TE D

Keywords: Unconventional Liquid Reservoirs (ULRs); CO2 huff-n-puff; Production performance; Dual continuum model; discrete fracture network

AC C

EP

Unconventional Liquid Reservoirs (ULRs) present great challenges to develop because of low permeability. However, with the improvement of technologies in petroleum industry, hydraulic fractured horizontal wells have been successfully applied to economically produce ULRs. During the process of hydraulic fracturing, high-pressurized fluids together with high-strength proppants are injected into the reservoir, creating high-conductivity flow paths, or in other words, induced hydraulic fractures to the horizontal wellbore. Because of the complex interactions between hydraulic fractures and natural fractures, fluid flow paths become complex fracture networks rather than simple planar fracture. Researchers (Weng et al. 2011) have proposed to apply fracture geomechanics to simulate hydraulic fracture propagation in complex fracture networks. One of the important applications of fracture propagation models was to investigate the “Stress Shadow effect” (Kresse et al. 2012), which provides useful guidance on hydraulic fracturing design and optimization. One of the main assumptions of current fracture modeling approaches is that reservoir fields are produced through the fracture networks including both hydraulic and natural fractures, and matrix rock contribution to fluid flow is negligible. Such assumption is reasonable since evidence (Gale et al. 2014, Gale et al. 2007) has demonstrated the existence of natural fractures in ULRs as well as its crucial impact on production performance. When it comes to characterization of natural fractures, Kim (2007) proposed stochastic algorithms to model density, length and orientation distributions of natural fractures by resorting to various data resources such as core, log, and outcrop data. Another important information that we can use to characterize natural fractures comes from microseismic data during hydraulic fracturing monitoring process. Even

ACCEPTED MANUSCRIPT 2

AC C

EP

TE D

M AN U

SC

RI PT

though there are still lots of uncertainties about the source mechanism of microseismic events (Cipolla et al. 2012, Warpinski et al. 2013a), several approaches (Cipolla et al. 2012, Sotelo Gamboa 2014) have already showed results in microseismic-constrained fracture network generation. The biggest concern for fracture characterization lies in its uncertainties, and we usually end up with more than one fracture realizations matching with the production history but yielding very limited production prediction. It is crucial to perform sensitivity studies not only on fracture properties but also fracture discretization related parameters. Sun and Schechter (2014) investigated stochastically generated natural fractures on production performance of shale gas reservoirs. On the other hand, unstructured grids are flexible to model complex fracture geometries, and Sun et al. (2015) showed that fracture discretization related parameters are less important than fracture characterization properties. The CPU performance can be improved by coarsening the reservoir background, reducing unnecessary refinement around fractures, and using hybrid-gridding approach where there are unstructured grids surrounding fractures and structured grids in the reservoir background.. So far, many aspects have been investigated such as vertical extruded PEBI grid (Prévost et al. 2005, Verma and Aziz 1997), intersection gridding technique (Branets et al. 2008), high-resolution gridding application (Olorode et al. 2012), optimization-based approach (Romain et al. 2011), flow-based grid generation (Moog 2013), discrete fracture simulation (Karimi-Fard et al. 2003), and workflow integration (Caumon et al. 2004), CPU performance analysis of unstructured grids (Sun et al. 2015). Furthermore, from both field observations and reservoir simulation, URLs present very steep decline. Hydraulically fractured horizontal wells usually die out in a few years once the fluids within the fractures are produced, which is followed by a very long production history with extremely low production rate. Such characteristic requires Enhanced Oil Recovery (EOR) technique to restimulate the depleted wells. So far, the potential of CO2 huff-n-puff in ULRs has been studied for simple fracture geometry or dual continuum models. Hawthorne et al. (2013) proposed five conceptual steps for CO2 huff-n-puff based EOR approach in unconventional reservoirs. In order to verify the proposed five-step mechanisms, they performed a series of laboratory CO2 exposure experiments and showed that CO2 effectively recovered oil from the Bakken source rock and reservoir rock. According to Hawthorne et al. (2013), CO2 first flows into and through the fractures getting in touch with the matrix surface. When the pressure gradient is high, CO2 permeates from the fracture into the matrix rock driven by the pressure difference. As a result, the oil gets swelled and extruded back into the fracture. When the pressure gradient becomes small, the diffusion rather than convection dominates the flow. With the force of concentration gradient, CO2 diffuses into the matrix pores, reducing oil viscosity and bring oil back into the fracture. In order to evaluate the CO2 EOR potential in simple fracture geometries of ULRs, several studies (Mohanty et al. 2013, Wan et al. 2013) have been carried out. However, there is not too much work which investigates CO2 huff-n-puff in complex fracture networks. First, previous work (Ghaderi 2013) simplified the fracture geometry into either planar or orthogonal, which was not common in the field. Second, commercial simulators with structured grids are not capable to discretize the domain with arbitrary fracture orientation and non-uniform fracture aperture, and new discretization scheme has to be developed and tailored for fractures. Third, CO2 huffn-puff involves complicated phase behavior, a thorough core-scale to field-scale study is required to reveal discrepancies of recovery mechanisms in different scales. Complex fracture geometry might complicate the recovery mechanism further due to high heterogeneous permeability paths. Fourth, CO2 huff-n-puff in dual continuum approaches with history-matched parameters might introduce inaccuracies due to over-simplified fracture distribution, which should be evaluated against solutions of explicit fracture models. To these ends, we will first present model development, which includes two approaches to generate complex fracture networks. The fracture discretization scheme will be devoted to new algorithms and techniques to handle complex fracture networks with unstructured grids. Then, core-scale CO2 huff-npuff simulation will be performed to history-match lab results of two preserved Barnett core plugs. The same rock and fluid properties will be used for a simple fracture network to investigate complex CO2

ACCEPTED MANUSCRIPT 3

phase behavior during field-scale CO2 huff-n-puff. Moreover, the developed fracture discretization scheme will be validated against well-known solutions, and then applied for complex fracture networks. The microseismic-constrained fracture generation model will be compared with the single-stage fracture network model. The fractal-based fracture generation models will be compared with dual continuum models.

RI PT

2. Model Development 2.1 Complex fracture network generation based on fractal theory

M AN U

SC

Unlike simple planar, orthogonal, or wire-mesh fractures, complex fracture network or discrete fracture network (DFN) consists of non-orthogonal plane with low-angle intersections and non-uniform aperture distributions. Complex fracture network usually results from complicated interactions with insitu natural fractures. In addition, it has been shown that natural fractures present scale-independent fractal features in several studies (Fardin et al. 2001, Kulatilake and Um 1999), and fractal theory was thus applied to combine different data sources such as outcrop studies, image logs and core data for generating more realistic stochastic fracture networks. Kim and Schechter (2009) developed fractal discrete fracture network (FDFN) generation code for natural fractures without performing reservoir simulation to evaluate the impact of hydraulic fractures on well production performance. In this paper, we not only introduce hydraulic fractures and horizontal wells, but also apply a new gridding scheme to discretize the fractured domain for CO2 huff-n-puff simulation.

2.2 Complex fracture network generation based on micro-seismic information

AC C

EP

TE D

Williams-Stroud et al. (2013) suggests that we use microseismic data recorded during hydraulic fracturing treatments to constrain the modeling of in-situ natural fractures. In this study (Sotelo Gamboa 2014), two assumptions were used to generate stochastic complex fracture models by incooperating micro-seismic data, core data, and logging data. The first assumption is that microseismic event locations constrain the location of natural fractures, i.e., only one fracture pass through a microseismic event location, because reactivation of plane of weaknesses is considered as one of the main mechanisms of induced microseismicity (Warpinski et al. 2013b). The other assumption is that there will be unconnected fractures to the main hydraulic-natural fracture network, which is a reasonable outcome that agrees with suggested remote reactivation mechanisms of pre-existing discontinuities (Dahi Taleghani and Olson 2014, Warpinski et al. 2013b). In addition, fractures are assumed to be planes that are perpendicular and fully penetrating into a horizontal pay layer. Fracture length is assumed to follow a Power Law distribution, and fracture strike to follow a Fisher distribution. For a natural fracture network with two fracture sets, the ratio between the number of microseismic events for each set is assumed the same as the ratio between the number of fractures per set. During fracture generation process, first of all, the locations of microseismic event will be perturbed to honor source uncertainties. Then, the microseismic event locations will be categorized into different fracture sets. Fracture length and strike will be assigned based on distribution functions. Finally, natural fractures are interconnected with hydraulic fracture branches to form complex fracture networks based on material balance equations.

2.3 Development of fracture discretization scheme Fig. 1 shows a workflow that we proposed in this study to discretize and simulate complex fracture networks. In total, there are seven components—fracture characterization, reservoir and fracture properties, fracture discretization, reservoir simulation, visualization, sensitivity analysis and history

ACCEPTED MANUSCRIPT 4

Reservoir and Fracture Characterization Reservoir Properties

Fracture Properties

Meshing Parameters

Reservoir Simulation

SC

Fracture Discretization

RI PT

matching, and EOR application. The preprocessor reads in reservoir properties, fracture properties, and user-defined meshing parameters to generate simulation input files for connection-list based simulators. The postprocessor reads in simulation results and performs visualization for both time-dependent and time-independent parameters. Sensitivity analysis generates numerous fracture realizations, and investigates the effect of both natural fractures as well as simulation grids on production performance. Among those realizations, the most possible fracture network can be chosen by history matching with production data, and then being used for production forecast. Finally, different EOR application scenarios can be investigated in the history-matched fracture network model.

M AN U

Visualization(post-processor) Sensitivity Analysis and History Matching EOR Application

Fig. 1—Integrated workflow for characterization and simulation of unconventional reservoirs. Input fracture geometry

TE D

Project 3D to 2D

Fractures

AC C

EP

Nonuniform fracture aperture

Refinement around fractures

Fracture intersection and clustering

Fixed points around fractures

Reservoir background Optimization for background mesh density and orientation Flexible points for reservoir background

All Voronoi Cell All Voronoi Centers Cell Centers 2D Voronoi grid 2.5D Voronoi grid Compute simulator input

Fig. 2—Workflow for discretization of complex fracture networks with Voronoi or PEBI grids.

The developed fracture discretization scheme can handle non-orthogonal, low-angle intersections of extensively clustered fractures with non-uniform aperture distributions. As seen in Fig. 2, the detailed workflow within the fracture discretization component consists of five steps: projection of 3D to 2D; calculation of fixed and flexible points; construction of 2D Voronoi; construction of 2.5D Voronoi grid where the grid boundary is vertical; and computation of simulator input. Firstly, by projecting input 3D fracture geometry onto the horizontal plane, the 3D meshing problem is reduced to a 2D meshing problem. This projecting process is justified when the fractures are almost vertical, and thus complicated 3D Voronoi cells can be avoided. Secondly, fractures and reservoir background are treated separately.

ACCEPTED MANUSCRIPT 5

RI PT

For the fractures, fixed Voronoi cell centers are computed to resolve fracture refinement, fracture intersection, and clustering, as well as variable fracture aperture. Next, flexible Voronoi cell centers are computed for reservoir background from force-based optimization algorithms, where different reservoir background mesh density, orientation, and type can be obtained by iteratively updating the flexible Voronoi cell centers. Thirdly, once all of the Voronoi cell centers are calculated, the 2D PEBI grid can be constructed. Fourthly, 2D grid is extruded vertically to each geological layer to obtain the 2.5 PEBI grids. Finally, simulator inputs are prepared such as the pore volume of each PEBI cell, transmissibility between two PEBI cells, cell center depth, as well as well-related properties.

3. Results and Discussions

3.1 Simple fracture geometry

M AN U

SC

This study deals with two issues for modeling CO2 huff-n-puff in complex fracture networks, which are complex fracture geometries and complex CO2 huff-n-puff phase behaviors. In the previous section, we introduced how we computed complex fracture geometries using either fractal or microseismicconstraint models as well as the new developed fracture discretization scheme. Regarding the complex CO2 huff-n-puff phase behaviors, we performed a systematic study from core-scale lab experiment, to core-scale simulation, to field-scale simulation. In this paper, we will only present the field-scale simulation results in the first part of this section. Then, we will show the effect of both complex fracture geometry and CO2 phase behaviors.

AC C

EP

TE D

A synthetic field case was developed to study oil recovery mechanisms on the field scale. To save computational time, a single stage with one fracture on a horizontal well was modeled. The main focus is to use the properties from the core model (Zou 2015) to study complex CO2 phase behaviors. The stimulated reservoir volume (SRV) of one stage can be viewed in the schematic shown in Fig. 3.

Fig. 3—Schematic of the SRV around the fracture highlighted by the red rectangle.

The field-scale input parameters such as rock and PVT properties were obtained from the cores-scale simulation of experimental work (Zou 2015). In short, we will fully explore the effect of production pressure, permeability and porosity of both matrix and fractures, capillary pressure, diffusion coefficient, CO2 injection beginning time, injection rate, injection pressure, injection cycle length, soaking length, and number of cycles on well performance.

3.1.1 Rock and PVT model The fluid model for the core-scale simulation was built based on the dead oil data from the reservoir. To model the fluid in the reservoir with little knowledge of the live oil properties, methane is combined

ACCEPTED MANUSCRIPT 6

with the dead oil components to create a live oil model to be used in the field model. Table 1 contains fluid properties of the recombined live oil model. The bubble point pressure at 150 ºF occurs at 1200 psi. Component

Mole Fraction, %

Critical Pressure, psi

Critical Temperature, o R

Acentric Factor

Molecular Weight, g/mol

CO2

1.00E-06

1055.88

547.56

0.23

44.01

C1

0.34

667.20

343.08

0.01

16.04

COMP1

0.22

497.47

1007.532

0.25

COMP2

0.20

347.60

1068.05

0.40

124.12

COMP3

0.15

248.97

1212.01

0.62

181.06

COMP4

0.09

136.43

700.45

1.00

259.99

Pseudo-Fracture Width, ft. Pseudo-Fracture Permeability, mD

150 1

30

45

M AN U

Fracture Porosity, %

63 x 13 x 10

SC

Grid Dimension Fracture Half-length, ft.

RI PT

Table 1—Live oil composition and properties.

95.33

Matrix Permeability, mD

1.00E-04

Matrix Porosity, %

6

Initial Water Saturation, %

10

Initial Reservoir Pressure, psi

3000

Reservoir Temperature, F

150

Pay Thickness, ft.

100

HCPV, MSTB 40.72 Table 2— Simulated reservoir properties.

TE D

The live oil has lower Minimum Miscibility Pressure (MMP) than the dead oil due to the change in composition. Using a commercial MMP module, the MMP is estimated to be 1417 psi, and First Contact Miscibility (FMC) is 1789 psi at 150 °F. The reservoir properties are listed in Table 2. The diffusion coefficient in Table 3 which was obtained from core-scale simulation will be used for the field-scale simulation.

Gas, ft2/day Oil, ft2/day

CO2

COMP1

EP

Component

COMP2

COMP3

COMP4

2.69E-05

1.33E-06

1.20E-06

8.43E-07

6.32E-07

9.18E-06

1.09E-06

1.01E-06

8.12E-07

6.35E-07

AC C

Table 3—Diffusion coefficients used to match the ultimate recovery volume in field units.

3.1.2 Base case setup

The bubble-point pressure for the recombined live oil is at 1200 psi, which is below the MMP value. In this study, since the reservoir is slightly over pressured, the initial pressure of the reservoir is 3000 psi. Cases with different production pressures were considered for the primary depletion as well as the CO2 huff-n-puff. The bottom-hole flowing pressure for different cases are 2000 psi, which is above the FMC pressure; 1550 psi, which is slightly above the MMP; 1300 psi, which is below the MMP but above the bubble-point pressure; and 1000 psi, which is below the bubble point pressure. Besides, 1 cycle of CO2 was injected following the primary depletion cases. The injection time was on day 1000 for all cases. The injectors were constrained by a surface rate of 500 Mscf/day with BHP constraint at 3500 psi, which was assumed to be below reservoir fracture pressure at the corresponding depth. After 30 days of injection and 15 days of soaking, the production took place at the same pressure as the corresponding primary depletion BHP. Table 4 shows the producer BHP, cumulative oil before and after

ACCEPTED MANUSCRIPT 7

huff-n-puff, and incremental percentage of the recoverable oil up to 5000 days. Fig. 4 shows cumulative oil production of all eight cases under pressure depletion and CO2 huff-n-puff up to 5000 days. Cum Oil, STB

Cum Oil after Huff-n-Puff, STB

Case 1

1000

1674.12

1841.53

10.0

Case 2

1300

910.13

942.56

3.56

Case 3

1550

758.02

769.95

1.57

2000

503.60

512.08

Case 4

Oil Incremental Percentage, %

RI PT

Producer BHP, psi

1.68

Table 4—Cases with different producer BHP and the corresponding recovery factors. 2500

SC

Case 1 Depletion Case 2 Depletion

1500

M AN U

Cumulative Oil Production, STB

2000

1000

500

0 1000

2000

3000

TE D

0

4000

Case 3 Depletion Case 4 Depletion Case 1 1 Cycle Case 2 1 Cycle Case 3 1 Cycle Case 4 1 Cycle

5000

Time, day

Fig. 4—Cumulative oil production for the primary depletion cases with different producer BHP.

AC C

EP

The oil incremental percentage for the cases in Table 4 is modest as is expected for huff-n-puff. In reality, natural fractures may exist in these reservoirs so that production should be higher than the calculated value in this model depending on the intensity of the natural fractures. From Fig. 4, one can see that Case 1 Depletion and Case 1 Cycle huff-n-puff show the largest production, as well as the largest recovery increment among all cases. This is because two-phase flow occurs when producing at pressures below the bubble point. The pressurizing effect of injected gas is more pronounced. When production occurs at pressure below the bubble point following the injection, gas comes out of the solution and provides extra energy for the flow. The other cases have very little increment in recovery even when CO2 is mixed with the oil phase. This is because the fractures were saturated with oil instead of fully saturated with pure CO2. This greatly reduced the chemical potential of CO2 existing in the fractures and very little CO2 was able to reach further into the nano-Darcy permeability matrix with the small diffusion coefficient values obtained from laboratory experiment simulations. Therefore, the production after CO2 injection showed very minimal improvement because the primary change of oil properties still exists in the fracture regions. The improvement of production was mainly due to the viscosity reduction of the oil existed in the fractures. In the next section, Case 1 with 1 Cycle of huff-npuff was chosen to be the base case for sensitivity study of other parameters.

3.1.3 Sensitivity analysis

ACCEPTED MANUSCRIPT 8

Parameters

Default

1

Matrix Porosity

0.06

2

Matrix Permeability

3

Fracture Porosity

4

Fracture Permeability, mD

30

5

Capillary Pressure

w/

6

Diffusion

w/

7

Injection Time, Days

1000

9

Injection Rate, MSCF

10

Injection Pressure, psi

11

Injection Length, days

12

15 16 17 18

Soaking Time, days

AC C

14

Number of Cycles

Oil Incremental Percentage, %

2395.67

30.09

1.00E-04

1.00E-03

3880.47

110.72

0.45

0.2

1591.93

-13.55

200

1961.55

6.52

w/t

1841.53

0

w/t

1841.53

0

2000

1886.96

2.47

500

1804.35

-2.02

500

5000

1927.55

4.67

3500

1400

1734.22

-5.83

30

EP

13

Cum Oil after Huff-n-puff, STB

0.12

TE D

8

Sensitivity

M AN U

Case #

SC

RI PT

Similarly, as the previous section, the effect of other parameters was investigated on both cumulative oil production and oil incremental percentage. Table 5 shows comparison of incremental percentage of the recoverable oil for all the cases. The third column represents the default input parameters from the previous Case 1, for which the cumulative oil production after huff-n-puff equals to 1841.53 bbl. The incremental percentage for recovery oil was computed by column 5 minus 1841.53, which is divided by 1841.53. In addition, among all the parameters in Table 5, we notice that the matrix/fracture porosity and permeability are the most significant parameters, and the higher the matrix/fracture porosity and permeability, the higher the incremental oil. On the other hand, diffusion, capillary pressure, and soaking time has almost zero effect on the field model. For the field model, convection due to pressure drawdown can be considered to be the dominating mechanism instead of diffusion, which is different from the core-scale simulations where the diffusion is dominating mechanism for oil recovery. Furthermore, regarding injection operating parameters, the later the initial injection time, the larger the injection rate, the higher the injection pressure, the longer the injection length, and the more the cycles, the more chances to repressurize the reservoir and provide gas drive to assist production under immiscible condition. However, economic analysis needs to be investigated in order to estimate the optimal CO2 operating parameters.

15

1

10

1831.27

-0.56

60

1866.86

1.38

200

1896.77

3.00

1

1841.53

0

15

1841.53

0

100

1841.53

0

2

1950.442

5.91

3

2051.36

11.39

Table 5—A brief summary of sensitivity analysis of input parameters.

3.2 Complex fracture geometry 3.2.1 Model validation In order to validate the proposed fracture discretization scheme for modeling complex fracture networks, a detailed CO2 huff-n-puff comparison between developed facture discretization approach (unstructured) and conventional discretization approach (structured) will be performed. Table 6 to Table 9 summarize simulation input parameters, which include reservoir and fracture properties, relative permeability curves, EOS input parameters, and binary interaction coefficients. The reservoir consists of 40 vertical and orthogonal hydraulic fractures, which will be discretized by a commercial simulator

ACCEPTED MANUSCRIPT 9

using conventional Tar-tan grid, as well as by the developed unstructured grid discretization approach. Fig. 5a is the top view of the unstructured PEBI grid, and Fig. 5b is the top view of structured Cartesian grids by the commercial simulator. In total there are 290,000 cells for the structured grid, and 108, 805 cells for the unstructured grid.

0.06

Reservoir Length, ft.

5000

1.00E-04

Reservoir Width, ft.

1800

Rock Compressibility, psi

3.00E-06

Reservoir Height, ft.

100

Hydraulic Fracture Permeability, md

100

Reference Datum Depth, ft.

12050

Hydraulic Fracture Half Length, ft.

240

Reference Pressure, psi

6000

Matrix Permeability, md -1

Hydraulic Fracture Width, ft.

RI PT

Matrix Porosity

3

0.01

62.4

-1

0.3

Water Compressibility, psi

SC

Hydraulic Fracture Porosity

Water Density, lb./ft

3.36E-06

Reservoir Top Depth, ft.

12000

Water Viscosity, cp

0.2053

Three-phase Kro Calculation Method

STONE 2

Water Formation Volume Factor

1.0803

Sw

Krw

M AN U

Table 6—Reservoir and Fracture Properties Krow

Sg

Krg

Krog

0

1

0

0

1

0.3

0

0.656

0.05

0.001

0.826

0.356

0.002

0.385

0.113

0.009

0.633

0.413

0.009

0.208

0.175

0.038

0.465

0.469

0.021

0.1

0.238

0.084

0.323

0.525

0.038

0.041

0.3

0.15

0.207

TE D

0.25

0.059

0.013

0.363

0.234

0.116

0.079

0.005

0.375

0.255

0.103

0.638

0.084

0.003

0.425

0.338

0.052

0.694

0.115

0

0.488

0.459

0.013

0.75

0.15

0

0.55

0.6

0

1

1

0

0.75

0.8

0

EP

0.581 0.625

Table 7—Input oil and gas relative permeability curves

Molar Weight, g/gmol

Critical Temperature, K

Critical Pressure, psi

Critical Gas Compressibility Factor

Acentric Factor

CO2

44.01

304.2

1070.16

0.2736

0.2250

N2-C1

16.21

189.67

665.028

0.2898

0.0084

C2-C4

44.79

412.17

639.303

0.2561

0.1481

C5-C7

83.46

556.92

554.043

0.2618

0.2486

C8-C12

120.52

667.52

456.288

0.2488

0.3279

C13-C19

220.34

673.76

283.563

0.2830

0.5672

C20-C30

321.52

792.4

226.086

0.2701

0.9422

AC C

Component

Table 8—Compositional data for the Peng-Robinson EOS equation

COMPONENT

CO2

N2-C1

0.1013

C2-C4

0.1317

0

C5-C7

0.142

0

N2-C1

C2-C4

0

C5-C7

C8-C12

C13-C19

ACCEPTED MANUSCRIPT 10 C8-C12

0.1501

0

0

0

C13-C19

0.1502

0

0

0

0

C20-C30

0.1503

0

0

0

0

0

Table 9—Binary interaction coefficient of each component

TE D

M AN U

SC

RI PT

The well first produces for 7665 days at a constant bottom-hole pressure of 3000 psi, which is followed by ten huff-n-puff cycles. Each cycle consists of 6 months’ CO2 injection at a rate of 2 MMscf/d, three months’ soaking, and one year’s production. Total simulated production time is around 15,000 days.

EP

Fig. 5—Comparison of top view between structured Cartesian grids (290,000) and unstructured 2.5D PEBI grids (108,805 cells).

AC C

Fig. 6 shows cumulative oil production up to 15,000 days. The black curve from the developed unstructured discretization approach shows a reasonable match with the red curve of the structured grids by the commercial simulator. The blue curve represents the case where there is no huff-n-puff. In addition, a good match of oil production rate between two discretization approaches is illustrated in Fig. 7. During each production cycle after CO2 soaking, we see good matches between unstructured grids and structured grids. Therefore, the developed discretization scheme is accurate enough to reproduce results of commercial simulators.

ACCEPTED MANUSCRIPT

RI PT

11

TE D

M AN U

SC

Fig. 6—Comparison of cumulative oil production between structured and unstructured discretization approaches.

Fig. 7—Comparison of oil rate between structured and unstructured discretization approaches.

EP

3.2.2 Comparison between explicit and dual continuum models

AC C

Dual continuum models can be used to represent simple sugar-cube concepts for naturally fractured reservoirs. For modeling complex fracture networks, dual continuum approaches are able to match production history by tuning parameters such as fracture spacing and fracture permeability without explicit honoring fracture geometry information, resulting in inaccurate pressure graphs and oil recovery. In this section, we will see if we can apply the dual continuum models to represent complex fracture networks or discrete fracture networks (DFNs) during CO2 huff-n-puff. The fractal-based approach was applied to generate stochastically distributed natural fractures. Table 7 summarizes input parameters for generating three DFN realizations. Fig. 8 shows reservoir and fracture geometry as well as fracture discretization. Each DFN model consists of 40 vertical hydraulic fractures, and stochastically generated natural fractures. Number of Natural Fracture Sets Minimum Fracture Length (ft.) Fracture Fractal Density Fracture Fractal Length Constant Density Term Scale Ratio Multifractal Dimension

2 0.03 * reservoir domain size 1.9 1.8 1.5 3 2

ACCEPTED MANUSCRIPT 12

Table 10—Input parameters for generating stochastic discrete fracture networks.

TE D

M AN U

SC

RI PT

With the same simulator input parameters as seen in the previous section, we performed CO2 huff-npuff for all three DFN models. Fig. 9 shows the comparison of cumulative oil production among three DFN models, the previous DFN model without natural fractures (Fig. 5), and a corresponding dual porosity model. The blue curve represents cumulative oil production for the DFN model without natural fractures. The pink curve, green curve, red curve represent the three DFN models with natural fractures, respectively. The black curve is the matched corresponding dual porosity model with parameters listed in Table 8. Before huff-n-puff, the DFN model without natural fractures show lower cumulative oil production than the other models. However, after the huff-n-puff treatment starts, we see a significant difference between the DFN model without natural fractures (blue curve) and the three DFN models with natural fractures. The existence of natural fractures increases the sweeping efficiency of injected CO2. The regions far away from hydraulic fractures can be swept by the injected CO2, and therefore more oil is produced from the DFN model with natural fractures. The three DFN models with natural fractures show very similar cumulative oil production performance, because they were generated with the same input parameters. Besides, by varying natural fracture parameters of the dual porosity model, we found that natural fracture permeability of 0.00003 and fracture porosity of 0.000005 provides a reasonable match between three DFN models and the dual porosity model. Similarly, from oil production rate plot in Fig. 10, we reach the same conclusions. We were able to find a corresponding dual porosity model to match production performance of the complicated DFN models for CO2 huff-npuff studies. Fig. 11 illustrates pressure graphs of the three DFN models and the matched dual porosity model. The dual porosity model shows very smooth pressure diffusion front. However, for all three DFN models, the pressure front depends on how the natural fractures are distributed. We observe pressure depletion in the regions far away from the horizontal wellbore because of the existence of natural fractures. Even though we could roughly match cumulative oil production of the DFN models with the dual porosity model, the pressure diffusion graph is different, and thus inevitably introduces inaccuracies for production forecast. Therefore, when it comes to modeling complex fracture networks, explicit models not only accurately predict production performance, but also capture pressure front behavior. Dual continuum models such as dual porosity tends to smooth the pressure diffusion front, which is not accurate enough to model complex fracture networks. 0.06

Natural Fracture Spacing Lx, ft.

25

Matrix Permeability, md

0.0001

Natural Fracture Spacing Ly, ft.

25

Rock compressibility

3.00E-06

Natural Fracture Spacing Lz, ft.

20

Hydraulic Fracture permeability, md

100

Natural Fracture Porosity

0.00005

Hydraulic Fracture half length, ft.

240

Natural Fracture Permeability Kx, md

0.00003

Hydraulic Fracture width, ft

0.01

Natural Fracture Permeability Ky, md

0.00003

Hydraulic Fracture porosity

0.3

Natural Fracture Permeability Kz, md

0.00003

AC C

EP

Matrix Porosity

Table 11—Input parameters for the matched dual porosity model

ACCEPTED MANUSCRIPT

TE D

M AN U

SC

RI PT

13

AC C

EP

Fig. 8—Three realizations of stochastically generated fractured network. From the top to the bottom, DFN_1, DFN_2, and DFN_3.

ACCEPTED MANUSCRIPT 14

M AN U

SC

RI PT

Fig. 9—Comparison of cumulative oil production among three DFN models (DFN_1, DFN_2, and DFN_3), the matched dual porosity model (DP_NF), and the DFN model without natural fractures (DP) from Fig. 5.

AC C

EP

TE D

Fig. 10—Comparison of oil production rate among three DFN models (DFN_1, DFN_2, and DFN_3), the matched dual porosity model (DP_NF), and the DFN model without natural fractures (DP) from Fig. 5.

ACCEPTED MANUSCRIPT

RI PT

15

Fig. 11—Comparison of pressure graphs among three DFN models and the matched dual porosity model.

SC

3.2.3 CO2 huff-n-puff for the microseismic-constraint characterization approach

AC C

EP

TE D

M AN U

In section 3.1, we see an orthogonal hydraulic fracture in a single stage, which is not usually the case in the field. Because reservoir heterogeneities and intersections between induced and in-situ natural fractures usually result in complex fracture networks. It is necessary to investigate complex CO2 huff-npuff behavior in complex fracture networks. The compositional PVT fluid model is the same as seen in section 3.2.1. In addition, the fracture geometry of a single stage will be generated based on microseismic and core data using the algorithms in section 2.2. Detailed algorithms and assumptions are discussed in Sotelo Gamboa (2014). In Fig. 12, the green dots represent the locations of microseismic events. Red line segments indicate hydraulic fractures, and the blue ones are the connected natural fractures. Solid black line segment indicates the well trajectory. The rest black line segments are the disconnected natural fractures, which were assumed to make no contribution to well production. The well produces from both hydraulic fractures and connected natural fractures.

Fig. 12—Complex fracture geometry which is generated from microseismic and core data.

Table 9 shows input reservoir and fracture properties as well as CO2 injection related parameters. The total pore volume from the simulator shows around 1.48 MM bbl. reservoir barrels of fluids. Considering CO2 is injected for 30 days at a rate of 500 Mscf/day, the total injected CO2 volume is around 2.4 reservoir pore volume. Fig. 13 shows the unstructured mesh after applying the proposed fracture discretization algorithms in section 2.3, which are designed to honor and model complex fracture geometry. Fig. 14 to Fig. 16 show production performance as well as the pressure graph of the complex fracture network. From Fig. 15, the changes in oil saturation are below 0.01, which is due to the ultra-low reservoir permeability. Correspondingly, the reservoir pressure in Fig. 14 depletes significantly and reaches bottom-hole flowing pressure of 3000 psi around the fractures. Fig. 16 shows

ACCEPTED MANUSCRIPT 16

comparison of cumulative oil production (COP), oil production rate (QOP) between two cases with and without CO2 huff-n-puff. With CO2 huff-n-puff, a spike increase in oil production rate is observed, which is followed by a steep decline. After soaking, it takes 181 days for the production rate goes back to the case without huff-n-puff. Cumulative oil production shows better performance (26.7 M bbl. after huff-n-puff) with an incremental oil around 438.61 STB at the end of 5000 days’ production, which represents 1.64 % increase in oil production. 0.06

Reservoir Length, ft.

Matrix Permeability, md

0.0001

Reservoir Width, ft.

1011

RI PT

Matrix Porosity

451

1.48

Reservoir Height, ft.

300

3984

Reservoir Top, ft.

10000

Hydraulic Fracture width, ft.

0.01

Initial Reservoir Pressure, psi

6000

Hydraulic Fracture porosity

0.3

Initial Oil Saturation

0.75

CO2 Soaking length, days

15

CO2 Injection length, days

30

CO2 injection rate, Mscf/day

500

Production length before huff-n-puff, days

1000

SC

Total Pore Volume, MM bbl. Hydraulic Fracture permeability, md

TE D

M AN U

Table 12—Reservoir and Fracture Properties

AC C

EP

Fig. 13—Discretization of the complex fracture network with unstructured PEBI grids.

Fig. 14—Pressure graph at the end of the huff-n-puff simulation.

ACCEPTED MANUSCRIPT

SC

RI PT

17

TE D

M AN U

Fig. 15—Oil Saturation at the end of the huff-n-puff simulation.

AC C

4. Conclusions

EP

Fig. 16—Cumulative oil production and oil production rate at the end of the huff-n-puff simulation.

From sensitivity analysis of complex CO2 huff-n-puff behaviors, it was observed that gas expansion is an important mechanism during depletion, thus huff-n-puff is mostly effective below the bubble-point pressure. In addition, field-scale simulation shows that the incremental oil is sensitive to matrix porosity, matrix permeability, fracture porosity, fracture permeability, time of first injection, cycling length, injection rate, injection pressure, and number of cycles, and it is not sensitive to capillary pressure, diffusion coefficient, and length of soaking time. Due to the short duration of huff-n-puff, diffusion is negligible, whereas we expect diffusion to play a greater role during continuous injection. From numerical simulation of CO2 huff-n-puff in complex fracture geometries, which are generated from two fracture characterization approaches, we found that dual continuum models introduce inaccurate results because it cannot really resolve the impact of stochastically distributed fracture networks. Pressure graphs are completely off the trend even though production rates could be matched with DFN models by varying natural fracture related parameters. On the other hand, CO2 huff-n-puff in complex fracture networks was proven to be feasible by the proposed fracture discretization approach. Pressure depletion and residual oil saturation can be accurately resolved.

ACCEPTED MANUSCRIPT 18

5. References

AC C

EP

TE D

M AN U

SC

RI PT

Branets, L. V., Ghai, S. S., Lyon, S. L. et al. 2008. Challenges and Technologies in Reservoir Modeling. Communications in Computational Physics 6:23. Caumon, G., Lepage, F., Sword, C. et al. 2004. Building and Editing a Sealed Geological Model. Mathematical Geology 36 (4):405-424. 10.1023/B:MATG.0000029297.18098.8a. Cipolla, C. L., Maxwell, S. C., Mack, M. G. 2012. Engineering Guide to the Application of Microseismic Interpretations. Presented at the SPE Hydraulic Fracturing Technology Conference, The Woodlands, Texas, USA, 01/01/2012. 152165. 10.2118/152165-ms. Dahi Taleghani, A., Olson, J. E. 2014. How Natural Fractures Could Affect Hydraulic-Fracture Geometry. SPE Journal 19 (01):161-171. 10.2118/167608-PA. Fardin, N., Stephansson, O., Jing, L. 2001. The Scale Dependence of Rock Joint Surface Roughness. International Journal of Rock Mechanics and Mining Sciences 38 (5):659–669. http://dx.doi.org/10.1016/S1365-1609(01)00028-4. Gale, J. F. W., Laubach, S. E., Olson, J. E. et al. 2014. Natural Fractures in Shale: A Review and New Observations. AAPG Bulletin 98 (11):2165-2216. Gale, J. F. W., Reed, R. M., Holder, J. 2007. Natural Fractures in the Barnett Shale and Their Importance for Hydraulic Fracture Treatments. AAPG bulletin 91 (4):603-622. Ghaderi, S. M. 2013. Simulation Study of Co2 Injection and Storage in Alberta. PhD, University of Calgary, Calgary, Alberta. Hawthorne, S. B., Gorecki, C. D., Sorensen, J. A. et al. 2013. Hydrocarbon Mobilization Mechanisms from Upper, Middle, and Lower Bakken Reservoir Rocks Exposed to Co, 2013/11/5/. 10.2118/167200-MS. Karimi-Fard, M., Durlofsky, L. J., Aziz, K. 2003. An Efficient Discrete Fracture Model Applicable for General Purpose Reservoir Simulators. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 3-5 February. 79699. http://dx.doi.org/10.2118/79699-ms. Kim, T. H. 2007. Fracture Characterization and Estimation of Fracture Porosity of Naturally Fractured Reservoirs with No Matrix Porosity Using Stochastic Fractal Models. Doctor of Philosophy, Texas A&M University. Kim, T. H., Schechter, D. 2009. Estimation of Fracture Porosity of Naturally Fractured Reservoirs with No Matrix Porosity Using Fractal Discrete Fracture Networks. SPE Reservoir Evaluation & Engineering 12 (2):232– 242. SPE-110720-PA. http://dx.doi.org/10.2118/110720-PA. Kresse, O., Weng, X., Wu, R. et al. 2012. Numerical Modeling of Hydraulic Fractures Interaction in Complex Naturally Fractured Formations, 2012/1/1/. Kulatilake, P. H. S. W., Um, J. 1999. Requirements for Accurate Quantification of Self-Affine Roughness Using the Roughness–Length Method. International Journal of Rock Mechanics and Mining Sciences 36 (1):5– 18. http://dx.doi.org/10.1016/S0148-9062(98)00170-3. Mohanty, K. K., Chen, C., Balhoff, M. T. 2013. Effect of Reservoir Heterogeneity on Improved Shale Oil Recovery by Co Huff-N-Puff, 2013/4/10/. 10.2118/164553-MS. Moog, G. J. E. A. 2013. Advanced Discretization Methods for Flow Simulation Using Unstructured Grids. Olorode, O. M., Freeman, C. M., Moridis, G. J. et al. 2012. High-Resolution Numerical Modeling of Complex and Irregular Fracture Patterns in Shale Gas and Tight Gas Reservoirs. Presented at the SPE Latin America and Caribbean Petroleum Engineering Conference, Mexico City, Mexico, 16–18 April. http://dx.doi.org/10.2118/152482-ms. Prévost, M., Lepage, F., Durlofsky, L. J. et al. 2005. Unstructured 3d Gridding and Upscaling for Coarse Modelling of Geometrically Complex Reservoirs. Petroleum Geoscience 11 (4):339-345. 10.1144/1354079304-657. Romain, M., Caumon, G., Levy, B. et al. 2011. Building Centroidal Voronoi Tessellations for Flow Simulation in Reservoirs Using Flow Information. Presented at the SPE Reservoir Simulation Symposium, The Woodlands, Texas, USA, 21–23 February. http://dx.doi.org/10.2118/141018-MS. Sotelo Gamboa, E. 2014. Natural Fracture Characterization by Source Mechanism Estimation and SemiStochastic Generation of Discrete Fracture Networks Using Microseismic and Core Data. Master of Science, Texas A&M University, College Station (December 2014). Sun, J., Schechter, D. S. 2014. Optimization-Based Unstructured Meshing Algorithms for Simulation of

ACCEPTED MANUSCRIPT 19

AC C

EP

TE D

M AN U

SC

RI PT

Hydraulically and Naturally Fractured Reservoirs with Variable Distribution of Fracture Aperture, Spacing, Length and Strike. Presented at the SPE Annual Technical Conference and Exhibition, Amsterdam, The Netherlands, 27–29 October. http://dx.doi.org/10.2118/170703-MS. Sun, J., Schechter, D. S., Huang, C.-K. 2015. Sensitivity Analysis of Unstructured Meshing Parameters on Production Forecast of Hydraulically Fractured Horizontal Wells. Presented at the Abu Dhabi International Petroleum Exhibition and Conference, Abu Dhabi, UAE, 9-12 November. SPE-177480-MS. Verma, S., Aziz, K. 1997. A Control Volume Scheme for Flexible Grids in Reservoir Simulation. Presented at the SPE Reservoir Simulation Symposium, Dallas, Texas, 8–11 June. http://dx.doi.org/10.2118/37999-MS. Wan, T., Sheng, J. J., Soliman, M. Y. 2013. Evaluate Eor Potential in Fractured Shale Oil Reservoirs by Cyclic Gas Injection, 2013/8/12/. 10.1190/URTEC2013-187. Warpinski, N. R., Mayerhofer, M., Agarwal, K. et al. 2013a. Hydraulic-Fracture Geomechanics and Microseismic-Source Mechanisms. 10.2118/158935-PA. Warpinski, N. R., Mayerhofer, M., Agarwal, K. et al. 2013b. Hydraulic-Fracture Geomechanics and Microseismic-Source Mechanisms. SPE Journal 18 (4): 766–780. SPE-158935-PA. http://dx.doi.org/10.2118/158935-PA. Weng, X., Kresse, O., Cohen, C. et al. 2011. Modeling of Hydraulic-Fracture-Network Propagation in a Naturally Fractured Formation. Spe Production & Operations 26 (4):368-380. Williams-Stroud, S., Ozgen, C., Billingsley, R. L. 2013. Microseismicity-Constrained Discrete Fracture Network Models for Stimulated Reservoir Simulation. Geophysics 78 (1):B37-B47. 10.1190/geo2011-0061.1. Zou, A. 2015. Compositional Simulation of Co2 Enhanced Oil Recovery in Unconventional Liquid Reservoirs. M.S., Texas A&M University, College Station.

ACCEPTED MANUSCRIPT

• Numerical simulation shows that CO2 Huff-n-puff improves oil production from unconventional liquid reservoirs with complex fracture networks. • Two different fracture characterization techniques such as fractal-based and microseismic-constrained reduce uncertainties of fracture networks and thus

RI PT

improve accuracy of the simulation results. • Dual continuum models yields less accurate pressure behavior than explicit discrete fracture models even though production behavior could be history

AC C

EP

TE D

M AN U

SC

matched.