Compositional Simulation of Geological and Engineering Controls on Gas Huff-n-Puff in Duvernay Shale Volatile Oil Reservoirs, Canada

Duvernay shale is a world class shale deposit with a total resource of 440 billion barrels oil equivalent in the Western Canada Sedimentary Basin (WCSB). The volatile oil recovery factors achieved from primary production are much lower than those from the gas-condensate window, typically 5–10% of original oil in place (OOIP). The previous study has indicated that huff-n-puff gas injection is one of the most promising enhanced oil recovery (EOR) methods in shale oil reservoirs. In this paper, we built a comprehensive numerical compositional model in combination with the embedded discrete fracture model (EDFM) method to evaluate geological and engineering controls on gas huff-n-puff in Duvernay shale volatile oil reservoirs. Multiple scenarios of compositional simulations of huff-n-puff gas injection for the proposed twelve parameters have been conducted and effects of reservoir, completion and depletion development parameters on huff-n-puff are evaluated. We concluded that fracture conductivity, natural fracture density, period of primary depletion, and natural fracture permeability are the most sensitive parameters for incremental oil recovery from gas huff-n-puff. Low fracture conductivity and a short period of primary depletion could significantly increase the gas usage ratio and result in poor economical efficiency of the gas huff-n-puff process. Sensitivity analysis indicates that due to the increase of the matrix-surface area during gas huff-n-puff process, natural fractures associated with hydraulic fractures are the key controlling factors for gas huff-n-puff in Duvernay shale oil reservoirs. The range for the oil recovery increase over the primary recovery for one gas huff-n-puff cycle (nearly 2300 days of production) in Duvernay shale volatile oil reservoir is between 0.23 and 0.87%. Finally, we proposed screening criteria for gas huff-n-puff potential areas in volatile oil reservoirs from Duvernay shale. This study is highly meaningful and can give valuable reference to practical works conducting the huff-n-puff gas injection in both Duvernay and other shale oil reservoirs.


Introduction
Horizontal drilling and multi-stage hydraulic fracturing have facilitated to a large increase in shale oil production in North America through continuous drilling investment in recent years. However, the ultimate oil recovery of shale oil reservoirs is usually less than 10%, which is much lower than that of conventional reservoirs [1][2][3][4][5][6].
Although the research on enhanced oil recovery methods for conventional reservoirs has been in-depth, research on enhanced oil recovery (EOR) for unconventional reservoirs is still at the conceptual stage. In the past ten years, oil companies and related scientific represent the fracture segment in the physical domain. Each individual fracture could be discretized into several fracture segments by the matrix cell boundaries. The EDFM method will check the intersection between the complex fracture geometries and the original structured matrix grids. The complex fractures will be discretized into several small fracture segments by the boundaries of matrix grids; thus, the EDFM method will generate additional virtual grids for all small fracture segments. The non-neighboring connections (NNC) are used to represent the physical connections between the original matrix grids and the additional fracture grids. Correspondingly, the fluid transport between complex fractures, matrix, and wellbore can be simulated inside reservoir simulator using transmissibility and effective wellbore index. The equations used to calculate the transmissibility factor between different NNCs and wellbore index can be found in the work by Xu et al. [19,20]. In short, the volumetric flow rate of phase l between cells in an NNC pair is mainly based on two-point flux approximation and calculated by where q l is the volumetric flow rate of phase l, λ l is the relative mobility of phase l, T NNC is the NNC transmissibility factor, and ∆P is the potential difference between the two points.
For the matrix-fracture connection, the NNC transmissibility is calculated by where T NNC is the transmissibility, k NNC is the matrix permeability in the direction perpendicular to the fracture plane, A NNC is the contact area of the plane inside the matrix block, and d NNC is the average normal distance from the matrix block to the fracture plane. For the connections between fracture segments, T NNC is the harmonic average between the transmissibility factors of the two segments.
For the well-fracture intersection, the effective well index is calculated by Peaceman (1983) equation, as shown in Equation (3).
, r e = 0.14 L 2 where W I f is the well index representing well-fracture intersection, k f is the fracture permeability, w f is the fracture aperture, L S is the length of the fracture segment, H S is the height of the fracture segment, r e is the effective radius, and r w is the wellbore radius. The EDFM model is then incorporated into a compositional model. The workflow to combine the compositional model with the EDFM model is as follows: • Build the compositional model.

•
Get new dimensions using EDFM method. The EDFM method will generate additional virtual grids to represent fractures. Extra grids will be appended in X direction. Thus the number of total grids will increase. Number and positions of fractures are defined in this step. • Add non-neighboring connections to represent the physical connections between the original matrix grids and the additional fracture grids. Fracture permeability, aperture and dip angle are defined in this step. • Assign pore volume for fracture cells. Fracture porosity is defined in this step. • Define fracture well intersection to represent the fluid transport between complex fractures, matrix, and wellbore. • Eventually, a comprehensive compositional model in combination with the EDFM method is built.

Nanopore Confinement on Shale Reservoir
Zarragoicoechea and Kuz [26] obtained the variation formula of critical properties with pore radius based on actual experimental data fitting. Singh et al. [27] further verified the correlation of critical property deviations using Monte Carlo simulation. In this method, Monte Carlo simulation is performed in a standard grand-canonical ensemble, where the volume, chemical potential, and temperature remain constant, while the number of particles and energy fluctuate.
Critical temperature and pressure shift in nanopores can be calculated by the following formulas: σ LJ = 0.244 3 T c /P c (6) where ∆T c is the relative critical temperature shift caused by nanopore confinement; ∆P c is the relative critical pressure shift caused by nanopore confinement; r p is the pore radius; T c is the bulk critical temperature; T cp is the pore critical temperature; P c is the bulk critical pressure; P cp is the pore critical pressure; and σ LJ is the Lennard-Jones size parameter.

Overview of the Studied Reservoir
The Duvernay Formation is an Upper Devonian source rock in the Western Canadian Sedimentary Basin (WCSB). The Duvernay shale was deposited during a basin-fill sequence with relatively high organic-matter productivity and preservation during episodes of Late Devonian reef growth. Since 2011, the liquids-rich or volatile oil window has been targeted by most of the development wells in the Duvernay [23].
The studied reservoir is located at the volatile oil window of the Kaybob area. The reservoir is overpressured and the pressure coefficient is approximately 1.8. The reservoir thickness range is 30-45 m. Total organic carbon (TOC) of the reservoir ranges between 2% to 6%. The thermal maturity (Ro) range is 0.9-1.2%. The porosity averages around 5%, and matrix permeability is of hundreds nano Darcy [28]. The GOR of the volatile oil window in this region ranges from 2000 to 3500 scf/STB (stock tank barrel). Figure 1 shows the GOR of the Duvernay shale and the location of the subject well.
The Diagnostic Fracture Injection Test (DFIT) of the subject well in Duvernay shale indicates complex fractures mixed with hydraulic fractures and natural fractures near the wellbore. The density of natural fractures is 0.1 fracture per meter along the horizontal wellbore according to OBMI (Oil-Base Micro-Imager) image log interpretation. The strike of natural fractures is dominated by NW-SE, which is parallel with the trajectory of the horizontal well. The estimated formation pressure and permeability for the subject well are 63.84 MPa and 0.0028 mD, respectively.

Fluid Characterization
Eight pseudocomponents were assigned in this study: N 2 , CO 2 , CH 4 , C 2 H 6 , C 3 H 8 , C 4 , C 5 -C 6 , and C 7+ . The corresponding molar fractions were 0.0026, 0.0053, 0.5403, 0.1260, 0.0764, 0.0480, 0.0522, and 0.1492, respectively. We tune the Peng-Robinson equation of state (PR EOS) by the bulk fluid PVT data from laboratory reports such as Constant Composition Expansion (CCE) and Differential Liberation (DL) tests. The oil properties were determined for this well using the Peng-Robinson EOS and flash calculations under a reservoir temperature of 124 • C. Oil gravity was 45 • API, gas/oil ratio was 4000 scf/STB, formation volume factor was 1.65 RB/STB, the bubble point pressure was 3492 psi, and the minimum miscible pressure (MMP) was 22.48 MPa. Detailed data for the PR EOS and the binary interaction parameters are listed in Tables 1 and 2. The Diagnostic Fracture Injection Test (DFIT) of the subject well in Duvernay shale indicates complex fractures mixed with hydraulic fractures and natural fractures near the wellbore. The density of natural fractures is 0.1 fracture per meter along the horizontal wellbore according to OBMI (Oil-Base Micro-Imager) image log interpretation. The strike of natural fractures is dominated by NW-SE, which is parallel with the trajectory of the horizontal well. The estimated formation pressure and permeability for the subject well are 63.84 MPa and 0.0028 mD, respectively.

Fluid Characterization
Eight pseudocomponents were assigned in this study: N2, CO2, CH4, C2H6, C3H8, C4, C5-C6, and C7+. The corresponding molar fractions were 0.0026, 0.0053, 0.5403, 0.1260, 0.0764, 0.0480, 0.0522, and 0.1492, respectively. We tune the Peng-Robinson equation of state (PR EOS) by the bulk fluid PVT data from laboratory reports such as Constant Composition Expansion (CCE) and Differential Liberation (DL) tests. The oil properties were determined for this well using the Peng-Robinson EOS and flash calculations under a reservoir temperature of 124 °C. Oil gravity was 45°API, gas/oil ratio was 4000 scf/STB, formation volume factor was 1.65 RB/STB, the bubble point pressure was 3492 psi, and the minimum miscible pressure (MMP) was 22.48 MPa. Detailed data for the PR EOS and the binary interaction parameters are listed in Tables 1 and 2.

Numerical Modeling Set-Up with Complex Fractures
The subject well was hydraulically-fractured in 21 stages, with 3 clusters per stage along a lateral length of 1860 m. The hydraulic fracture pattern is assumed to be multiple bi-wing fractures with evenly fracture spacing, the number of fractures is 63. The initial pressure of this volatile oil reservoir is about 63.8 MPa. The primary production of the well consists of oil, gas, and water. The model size is set based on the lateral length, well spacing and reservoir thickness of the subject well. The lateral length, well spacing and reservoir thickness are 1860 m, 300 m, and 45 m, respectively. In order to make sure that the model can fully cover the well and avoid boundary effect, we set the model size to 2600 m × 610 m × 45 m, which corresponds to length, width, and height, respectively. Since the spacing of hydraulic fractures is 30 m, to ensure that there are at least two sets of grid blocks among each hydraulic fractures, the grid size along the x-axis (length) is set to 10 m. We set the same grid size for the y-axis (width). The grid size for the vertical axis is set to 5 m, which we believe it is enough for running fracture models. Thus, the number of grid blocks is 260 × 61 × 9, and the grid cell dimension is 10 m × 10 m × 5 m.
Both natural fractures and hydraulic fractures were considered in numerical modeling, and EDFM was performed to establish fracture models [29,30]. Based on the equivalent conductivity, the width of hydraulic fractures is set to 0.4 m and the permeability is 2.5 mD, that is, the conductivity of hydraulic fracture is 100 mD·cm. Considering the fracture density of the natural fractures, assuming that the density of natural fractures is 0.1/m, the permeability of the unstimulated area is 0.001 mD, and the permeability of the stimulated area is equal to the permeability of hydraulic fractures to 2.5 mD. According to core data from the subject well, the width of natural fractures is set to 0.004 m. Then, the conductivity of natural fracture in the unstimulated area is 0.0004 mD·cm. Figure 2 shows the numerical model with 63 hydraulic fractures (in red) and 1000 natural fractures (in grey). Considering the simulation time, the model mainly considers two scales of pores, nano-scale matrix systems, and micro-centimeter-scale fracture systems (including natural fractures and hydraulic fractures) [31,32]. By adjusting the critical properties of different components, two sets of EOS are set for different pore systems. The critical properties of the components corresponding to each pore radius were calculated according to Equation (4) to Equation (6). Table 3 lists the critical pressure and temperature of the eight pseudocomponents corresponding to the pore radius of 10 nm, 30 nm, 50 nm, 2000 nm and bulk phase. The base case of the pore radius for the matrix is set to 10 nm.  Considering the simulation time, the model mainly considers two scales of pores, nano-scale matrix systems, and micro-centimeter-scale fracture systems (including natural fractures and hydraulic fractures) [31,32]. By adjusting the critical properties of different components, two sets of EOS are set for different pore systems. The critical properties of the components corresponding to each pore radius were calculated according to Equation (4) to Equation (6). Table 3 lists the critical pressure and temperature of the eight pseudocomponents corresponding to the pore radius of 10 nm, 30 nm, 50 nm, 2000 nm and bulk phase. The base case of the pore radius for the matrix is set to 10 nm. Molecular diffusion can be described by Fick's law, which can be used to describe the characteristics of pore-scale mass transfer and diffusion caused by concentration differences [33]. It has a certain impact on the effect of gas injection and shale oil production. Molecular diffusion parameters can usually be obtained by fitting the results of laboratory experiments. The molecular diffusion characteristic value of the methane-propane system is about 2 × 10 −4 -10 −3 cm 2 /s (70 • C, 10-20 MPa). Usually other mechanisms will also be considered in the molecular diffusion mechanism, so this parameter is usually larger than the actual value and is usually used as tuning parameters for calibrating numerical models. The gas diffusion coefficients are usually an order of magnitude larger than the oil diffusion coefficients. Hamdi et al. [6] used diffusion coefficients of 0.01 cm 2 /s [16] and 0.001 cm 2 /s for the light components in the gas and oil phases, respectively. In this paper, we took Hamdi's assumption as reference, and set the diffusion coefficient of CH 4 as 0.01 cm 2 /s, then the diffusion coefficient of the other components were set according to their molecular weight. The molecular diffusion parameters of different components are shown in Table 4. It should be pointed out that the model we built is a single well model. Well spacing and well interference cannot be investigated using the model. This is mainly because we aim to discuss the geological and engineering control on huff and puff, however, a multiple well model could complicate the effect of these parameters.

History-Matching Parameters
Based on the numerical model constructed in this paper, history matching was conducted. Although the production data is noisy, a good history match of oil rate was achieved as shown in Figure 3. The history matched parameters are listed in Table 5, which are the basic reservoir and fracture inputs for production simulations. We use EDFM to efficiently and easily handle hydraulic fractures in the model. The flowing bottom-hole pressure (BHP) was applied for the reservoir simulation constraint.

History-Matching Parameters
Based on the numerical model constructed in this paper, history matching was conducted. Although the production data is noisy, a good history match of oil rate was achieved as shown in Figure 3. The history matched parameters are listed in Table 5, which are the basic reservoir and fracture inputs for production simulations. We use EDFM to efficiently and easily handle hydraulic fractures in the model. The flowing bottom-hole pressure (BHP) was applied for the reservoir simulation constraint.

Effect of Reservoir, Completion, and Depletion Parameters on Huff-n-Puff in Duvernay Shale
Oil recovery during gas huff-n-puff in shale reservoirs is a complex function of mass transfer (including miscibility, interfacial-tension, swelling, and molecular diffusion) and other controlling factors such as hydraulic fracturing and geological heterogeneities [3,6,21].
Based on the numerical model constructed and the history-matched parameters, we proposed a set of parameters that can represent hydraulic fracturing and geological heterogeneities. These parameters can be divided into three types: reservoir, hydraulic fracturing and primary depletion. Table 6 summarizes the type, name, unit, range and base case of the geological and engineering parameters for the simulation. To examine the effects of these parameters on gas huff-n-puff in Duvernay shale, 37 cases were simulated. We set fixed operational parameters for all the cases. The gas huff-n-puff process began after five years of primary depletion. We used rich gas as the injected gas, which is composed of 70% of CH 4 , 20% of C 2 H 6 , and 10% of C 3 -C 4 . The injection rate is set to 3.5 MMscf/d and the cumulative gas injection is 176.5 MMscf. Then, the well is brought back on production after soaking for 30 days. We shut the well when the oil rate declined to the rate before the gas huff-n-puff process (the oil rate after 5 years of primary depletion). The period from the end of soaking to the shut-in date is defined as the duration of the first huff-n-puff cycle. Cumulative Oil Production for Primary Depletion, incremental oil recovery, and gas usage ratio was calculated for each case.
The incremental oil recovery equals to incremental oil production with huff-n-puff divided by OOIP of the subject well. The OOIP of the subject well is calculated under the condition of 300 m well spacing, 1860 m lateral length, and 45 m reservoir thickness. The calculated OOIP of the subject well is 2682 Mbbl. Gas usage ratio is defined as total volume of injected gas divided by total barrels of improved oil production by gas injection.  Table 7 summarizes results for huff-n-puff with different reservoir thickness. What needs to be clarified is that we built five models for these five cases of reservoir thickness. The grid size of all of the cases is 10 m × 10 m × 5 m, and the difference among each case is the vertical number of grid blocks. In particular, we adjusted the cumulative volume of injected gas with respect to different reservoir thickness. The results of cumulative oil production, average reservoir pressure, incremental oil recovery, duration of the first huff-n-puff cycle and gas usage ratio were plotted in Figures 4 and 5. As can be seen from these figures, thick reservoirs produced more oil during primary depletion, whereas the lowest incremental oil recovery and highest gas usage ratio were observed in the thickest reservoir. We also noticed the duration of the first huff-n-puff cycle increased with the increase of reservoir thickness. Even if the incremental oil recovery was lower in thicker reservoirs, a minimum thickness of 30 m was proposed for gas huff-n-puff due to the higher primary production in thicker reservoirs.

TOC
The effect of TOC was quantified by adsorbed gas content which increased with the increase of TOC. As shown in Table 8, TOC of 2%, 4%, 6%, and 10% were examined in this

TOC
The effect of TOC was quantified by adsorbed gas content which increased with the increase of TOC. As shown in Table 8, TOC of 2%, 4%, 6%, and 10% were examined in this study. Figures 6 and 7 indicate that cases with higher TOC shows higher cumulative oil production due to the desorption of the adsorbed gas. The incremental oil recovery shows a negative correlation with TOC, but the variance among different TOC is minor. Due to the desorption effect, the duration shows a positive correlation with TOC.

TOC
The effect of TOC was quantified by adsorbed gas content which increased with the increase of TOC. As shown in Table 8, TOC of 2%, 4%, 6%, and 10% were examined in this study. Figures 6 and 7 indicate that cases with higher TOC shows higher cumulative oil production due to the desorption of the adsorbed gas. The incremental oil recovery shows  a negative correlation with TOC, but the variance among different TOC is minor. Due to the desorption effect, the duration shows a positive correlation with TOC.

Pore Radius of Shale Matrix
The nanopore confinement in shale reservoirs would lead to an increase in capillary pressure, which has a significant impact on the phase characteristics, critical properties,

Pore Radius of Shale Matrix
The nanopore confinement in shale reservoirs would lead to an increase in capillary pressure, which has a significant impact on the phase characteristics, critical properties, and physical properties of the fluid [34]. In numerical simulation, the effect of nanopore confinement can be calculated by modifying the equation of state.
The pores in Duvernay shale are mainly organic pores and the pore throats are of nano-scale, their radius is mainly 2-50 nm, and the radius of the organic pores is about 10 nm. We set the base case of pore radius to 10 nm. Additionally, the pore radius of 30 nm, 50 nm and 2000 nm were examined for gas huff-n-puff. It can be seen from Table 9, Figures 8 and 9 that the smaller the pore radius, the higher the cumulative oil production, but the difference is slight. This is due to the low density and low viscosity of the oil phase under the low pore radius, which is conducive to oil phase flow. The results imply that the effect of pore radius on gas huff-n-puff was marginal. Four cases of different matrix permeability were simulated and results were shown in Table 10, Figures 10 and 11. Lower matrix permeability could result in lower cumulative oil production and higher average reservoir pressure. Higher matrix permeability is beneficial for incremental oil recovery and gas usage ratio, but the difference is minor. However, higher matrix permeability could shorten the duration of the first huff-n-puff cycle due to its high flowing rate in the reservoir. The recommended matrix permeability for huff-n-puff is ≥5·10 −5 mD. The higher permeability is conducive to injecting gas to displace the fluid in the matrix.

Matrix Permeability
Four cases of different matrix permeability were simulated and results were shown in Table 10, Figures 10 and 11. Lower matrix permeability could result in lower cumulative oil production and higher average reservoir pressure. Higher matrix permeability is beneficial for incremental oil recovery and gas usage ratio, but the difference is minor. However, higher matrix permeability could shorten the duration of the first huff-n-puff cycle due to its high flowing rate in the reservoir. The recommended matrix permeability for huff-n-puff is ≥5·10 −5 mD. The higher permeability is conducive to injecting gas to displace the fluid in the matrix.  (a) (b) Figure 9. Correlation between pore radius of shale matrix and huff-n-puff results: (a) Cumulative oil production and incremental oil recovery for different pore radius of shale matrix; (b) Duration of the first huff-n-puff cycle and gas usage ratio for different pore radius of shale matrix.

Matrix Permeability
Four cases of different matrix permeability were simulated and results were shown in Table 10, Figures 10 and 11. Lower matrix permeability could result in lower cumulative oil production and higher average reservoir pressure. Higher matrix permeability is beneficial for incremental oil recovery and gas usage ratio, but the difference is minor. However, higher matrix permeability could shorten the duration of the first huff-n-puff cycle due to its high flowing rate in the reservoir. The recommended matrix permeability for huff-n-puff is ≥5·10 −5 mD. The higher permeability is conducive to injecting gas to displace the fluid in the matrix.  Figure 10. Effect of different matrix permeability with one huff-n-puff cycle on oil production and reservoir pressure: (a) Cumulative oil production for different matrix permeability; (b) Average reservoir pressure for different matrix permeability.  (a) (b) Figure 11. Correlation between matrix permeability and huff-n-puff results: (a) Cumulative oil production and incremental oil recovery for different matrix permeability; (b) Duration of the first huff-n-puff cycle and gas usage ratio for different matrix permeability.

Natural Fracture Density
The effect of gas huff-n-puff under different natural fracture densities was analyzed. Five scenarios are considered. Figure 12 shows the numerical models for different scenarios. The results are shown in Table 11, and Figures 13 and 14. It can be seen from the figures: The existence of natural fractures has a greater impact on the cumulative oil production for primary depletion which is 1.7 times compared with that without natural fractures.
When natural fractures are not considered, the oil rate after gas injection is low, but the duration of the first huff-n-puff cycle is much longer, which leads to a higher incremental oil recovery. When natural fractures exist, the density of natural fractures has a minor effect on the cumulative oil production for primary depletion. This is mainly due to the ultralow conductivity of natural fractures. According to core data from the subject well, the conductivity of natural fracture is 0.0004 mD·cm which is ultralow compared with that of hydraulic fractures (100 mD·cm). The simulation results indicate reducing fracture conductivity will significantly lowering primary depletion. The ultralow conductivity of natural fractures could result in the low impact of natural fracture density on primary depletion. Furthermore, the orientation of natural fractures is parallel to the horizontal well, so the hydraulic fractures can only connect limited natural fractures.
However, the incremental oil recovery increased significantly with the increase of natural fracture density. This can be attributed to a larger amount of injected gas available to contact the matrix area and extract hydrocarbons through natural fractures. The development of natural fractures is a favorable condition for gas huff-n-puff. The effect of gas huff-n-puff under different natural fracture densities was analyzed. Five scenarios are considered. Figure 12 shows the numerical models for different scenarios. The results are shown in Table 11, and Figures 13 and 14. It can be seen from the figures:        The effect of gas huff-n-puff under different natural fracture permeability was analyzed, and the results are shown in Table 12, and Figures 15 and 16. It can be seen from the figures: Natural fracture permeability shows a distinct positive correlation with cumulative oil production, and high oil rate can be achieved with high natural fracture permeability.
The incremental oil recovery decreases with the increase of natural fracture permeability. This is attributed to the high oil recovery for primary depletion and the oil abundance near wellbore has been reduced by primary depletion. Meanwhile, natural fractures The existence of natural fractures has a greater impact on the cumulative oil production for primary depletion which is 1.7 times compared with that without natural fractures.
When natural fractures are not considered, the oil rate after gas injection is low, but the duration of the first huff-n-puff cycle is much longer, which leads to a higher incremental oil recovery. When natural fractures exist, the density of natural fractures has a minor effect on the cumulative oil production for primary depletion. This is mainly due to the ultralow conductivity of natural fractures. According to core data from the subject well, the conductivity of natural fracture is 0.0004 mD·cm which is ultralow compared with that of hydraulic fractures (100 mD·cm). The simulation results indicate reducing fracture conductivity will significantly lowering primary depletion. The ultralow conductivity of natural fractures could result in the low impact of natural fracture density on primary depletion. Furthermore, the orientation of natural fractures is parallel to the horizontal well, so the hydraulic fractures can only connect limited natural fractures.
However, the incremental oil recovery increased significantly with the increase of natural fracture density. This can be attributed to a larger amount of injected gas available to contact the matrix area and extract hydrocarbons through natural fractures. The development of natural fractures is a favorable condition for gas huff-n-puff.

Natural Fracture Permeability
The effect of gas huff-n-puff under different natural fracture permeability was analyzed, and the results are shown in Table 12, and Figures 15 and 16. It can be seen from the figures:

Half-length of Hydraulic Fractures
Fracture half-length of 35 m, 60 m, and 85 m were simulated to examine its effect on gas huff-n-puff. The results are shown in Table 13, and Figures 17 and 18. It can be seen from the figures: Long fracture half-length has a significant contribution in shale reservoirs by increasing oil production during primary depletion. When the same volume of gas was injected, fracturing could weaken the effect of gas injection, and with the increase of fracture halflength, the incremental oil recovery could be decreased. This is because the injected gas needs to seep along the long fractures to the far well area, reducing the gas injection effect. Compared with the long fractures, if there are micro-fractures in the near well area, the effect of gas huff-n-puff is better.
The duration of the gas huff-n-puff cycle increases with the increase of fracture halflength. For long fracture half-length, the total volume of injected gas should be increased to improve the effect of gas huff-n-puff.  Natural fracture permeability shows a distinct positive correlation with cumulative oil production, and high oil rate can be achieved with high natural fracture permeability.
The incremental oil recovery decreases with the increase of natural fracture permeability. This is attributed to the high oil recovery for primary depletion and the oil abundance near wellbore has been reduced by primary depletion. Meanwhile, natural fractures improve the pressure drop near wellbore, which makes it difficult to form a high-pressure zone near the wellbore and leads to a shorter duration of the gas huff-n-puff cycle. Increasing the total volume of injected gas could improve the effect of gas huff-n-puff and yield more additional oil.

Half-Length of Hydraulic Fractures
Fracture half-length of 35 m, 60 m, and 85 m were simulated to examine its effect on gas huff-n-puff. The results are shown in Table 13, and Figures 17 and 18. It can be seen from the figures: Long fracture half-length has a significant contribution in shale reservoirs by increasing oil production during primary depletion. When the same volume of gas was injected, fracturing could weaken the effect of gas injection, and with the increase of fracture halflength, the incremental oil recovery could be decreased. This is because the injected gas needs to seep along the long fractures to the far well area, reducing the gas injection effect. Compared with the long fractures, if there are micro-fractures in the near well area, the effect of gas huff-n-puff is better.

Fracture Conductivity
"System leakage" is a very practical and interesting topic for huff and puff. Fracture conductivity is an indicator of "leakage" or sealing capability. We analyzed the effect of gas huff-n-puff under different hydraulic fracture conductivity, results are shown in Table  14, and Figures 19 and 20. It can be seen from the results that with the increase of fracture conductivity, the cumulative oil production for primary depletion and the incremental oil recovery increase, but the increase gradually flattens. When the fracture conductivity reaches 1000 mD·cm, there is a slight increase in cumulative oil production and oil recovery. So, it can be inferred that there is a limitation of conductivity. However, we believe the sealing capability of the reservoir more likely depends on the extension of hydraulic fractures and natural fractures instead of conductivity.

Fracture Conductivity
"System leakage" is a very practical and interesting topic for huff and puff. Fracture conductivity is an indicator of "leakage" or sealing capability. We analyzed the effect of gas huff-n-puff under different hydraulic fracture conductivity, results are shown in Table  14, and Figures 19 and 20. It can be seen from the results that with the increase of fracture conductivity, the cumulative oil production for primary depletion and the incremental oil recovery increase, but the increase gradually flattens. When the fracture conductivity reaches 1000 mD·cm, there is a slight increase in cumulative oil production and oil recovery. So, it can be inferred that there is a limitation of conductivity. However, we believe the sealing capability of the reservoir more likely depends on the extension of hydraulic fractures and natural fractures instead of conductivity. The duration of the gas huff-n-puff cycle increases with the increase of fracture halflength. For long fracture half-length, the total volume of injected gas should be increased to improve the effect of gas huff-n-puff.

Fracture Conductivity
"System leakage" is a very practical and interesting topic for huff and puff. Fracture conductivity is an indicator of "leakage" or sealing capability. We analyzed the effect of gas huff-n-puff under different hydraulic fracture conductivity, results are shown in Table 14, and Figures 19 and 20. It can be seen from the results that with the increase of fracture conductivity, the cumulative oil production for primary depletion and the incremental oil recovery increase, but the increase gradually flattens. When the fracture conductivity reaches 1000 mD·cm, there is a slight increase in cumulative oil production and oil recovery. So, it can be inferred that there is a limitation of conductivity. However, we believe the sealing capability of the reservoir more likely depends on the extension of hydraulic fractures and natural fractures instead of conductivity.
Wells with higher fracture conductivity is preferable for gas huff-n-puff. Fractures with high conductivity are conducive to gas injection into the formation to interact with crude oil. The recommended conductivity is greater than 100 mD·cm. Wells with higher fracture conductivity is preferable for gas huff-n-puff. Fractures with high conductivity are conducive to gas injection into the formation to interact with crude oil. The recommended conductivity is greater than 100 mD·cm.  Figure 19. Effect of different fracture conductivity with one huff-n-puff cycle on oil production and reservoir pressure: (a) Cumulative oil production for different fracture conductivity; (b) Average reservoir pressure for different fracture conductivity.
(a) (b) Figure 20. Correlation between fracture conductivity and huff-n-puff results: (a) Cumulative oil production and incremental oil recovery for different fracture conductivity; (b) Duration of the first huff-n-puff cycle and gas usage ratio for different fracture conductivity. Figure 19. Effect of different fracture conductivity with one huff-n-puff cycle on oil production and reservoir pressure: (a) Cumulative oil production for different fracture conductivity; (b) Average reservoir pressure for different fracture conductivity. Wells with higher fracture conductivity is preferable for gas huff-n-puff. Fractures with high conductivity are conducive to gas injection into the formation to interact with crude oil. The recommended conductivity is greater than 100 mD·cm.  Figure 19. Effect of different fracture conductivity with one huff-n-puff cycle on oil production and reservoir pressure: (a) Cumulative oil production for different fracture conductivity; (b) Average reservoir pressure for different fracture conductivity.

Type of SRV Region
(a) (b) Figure 20. Correlation between fracture conductivity and huff-n-puff results: (a) Cumulative oil production and incremental oil recovery for different fracture conductivity; (b) Duration of the first huff-n-puff cycle and gas usage ratio for different fracture conductivity.

Type of SRV Region
The SRV region is defined by modification of permeability in these grid cells. We multiply the primary permeability of matrix and fractures by a number of times to represent the effect of stimulation. Considering the stimulated volume around the main hydraulic fractures, we assume four types of SRV region to represent no stimulation around main hydraulic fractures (Figure 21A), stimulation within 10 m of the single wing ( Figure 21B), stimulation within 10 m of the bi-wing ( Figure 21C), and stimulation within 20 m of the single wing ( Figure 21D), respectively. The effect of the types of SRV region on gas huffn-puff was analyzed. The results are summarized in Table 15, and Figures 22 and 23, it can be seen from the figures that the types of SRV region show a very slight effect on gas huff-n-puff.
The SRV region is defined by modification of permeability in these grid cells. We multiply the primary permeability of matrix and fractures by a number of times to represent the effect of stimulation. Considering the stimulated volume around the main hydraulic fractures, we assume four types of SRV region to represent no stimulation around main hydraulic fractures ( Figure 21A), stimulation within 10 m of the single wing ( Figure  21B), stimulation within 10 m of the bi-wing ( Figure 21C), and stimulation within 20 m of the single wing ( Figure 21D), respectively. The effect of the types of SRV region on gas huff-n-puff was analyzed. The results are summarized in Table 15, and Figures 22 and 23, it can be seen from the figures that the types of SRV region show a very slight effect on gas huff-n-puff.

Magnification of Permeability in SRV Region
To further examine the effect of hydraulic fracturing, we set different permeability magnifications in the SRV region, which are 1×, 2×, 5×, and 8×. We increased the permeability in the SRV region (Orange grids in Figure 21C) by one time, two times, five times, and eight times. The results are summarized in Table 16, and Figures 24 and 25, from the figures it can be seen:

Magnification of Permeability in SRV Region
To further examine the effect of hydraulic fracturing, we set different permeability magnifications in the SRV region, which are 1×, 2×, 5×, and 8×. We increased the permeability in the SRV region (Orange grids in Figure 21C) by one time, two times, five times, and eight times. The results are summarized in Table 16, and Figures 24 and 25, from the figures it can be seen: As the permeability of natural fractures in the SRV region increases, the cumulative oil production during primary depletion increases, but because the area of the SRV region is small, the impact on the cumulative oil production is minor.
The increase of permeability in the SRV region has little effect on the increase in oil recovery and gas usage ratio during gas huff-n-puff.

Magnification of Permeability in SRV Region
To further examine the effect of hydraulic fracturing, we set different permeability magnifications in the SRV region, which are 1×, 2×, 5×, and 8×. We increased the permeability in the SRV region (Orange grids in Figure 21C) by one time, two times, five times, and eight times. The results are summarized in Table 16, and Figures 24 and 25, from the figures it can be seen: As the permeability of natural fractures in the SRV region increases, the cumulative oil production during primary depletion increases, but because the area of the SRV region is small, the impact on the cumulative oil production is minor.
The increase of permeability in the SRV region has little effect on the increase in oil recovery and gas usage ratio during gas huff-n-puff.   As the permeability of natural fractures in the SRV region increases, the cumulative oil production during primary depletion increases, but because the area of the SRV region is small, the impact on the cumulative oil production is minor.
The increase of permeability in the SRV region has little effect on the increase in oil recovery and gas usage ratio during gas huff-n-puff.
Cumulative oil production with one huff-n-puff cycle

Period of Primary Depletion
Primary depletion has been examined for gas huff-n-puff in the period of one year, three years, five years, and seven years. The results are summarized in Table 17,  The reservoir pressure drop increases with the increase in the period of primary depletion. The incremental oil recovery also increases with the period of primary depletion.
If the gas huff-n-puff starts with a short period of primary depletion, it could result in a poor huff-n-puff effect. This is because the high reservoir pressure is not conducive to gas injection and migration to the deep part of the reservoir. At the same time, it will lead to a too high gas-oil ratio during the production stage.
The proposed period of primary depletion is five years, and the reservoir pressure should be above saturation pressure to prevent condensation in the wellbore which could affect the oil recovery.

Period of Primary Depletion
Primary depletion has been examined for gas huff-n-puff in the period of one year, three years, five years, and seven years. The results are summarized in Table 17 The reservoir pressure drop increases with the increase in the period of primary depletion. The incremental oil recovery also increases with the period of primary depletion.
If the gas huff-n-puff starts with a short period of primary depletion, it could result in a poor huff-n-puff effect. This is because the high reservoir pressure is not conducive to gas injection and migration to the deep part of the reservoir. At the same time, it will lead to a too high gas-oil ratio during the production stage.
The proposed period of primary depletion is five years, and the reservoir pressure should be above saturation pressure to prevent condensation in the wellbore which could affect the oil recovery.

Period of Primary Depletion
Primary depletion has been examined for gas huff-n-puff in the period of one year, three years, five years, and seven years. The results are summarized in Table 17, and Figures 26 and 27, from the figures it can be seen: Cumulative oil production with one huff-n-puff cycle

BHP of Primary Depletion
The bottom hole pressure (BHP) controls the well production. The lower the BHP, the greater the production pressure difference, the higher the well production, and the faster the reservoir pressure drops. In shale reservoirs, multiple studies indicated that the production characteristics of tight reservoirs with high initial production and sharp decline in production will cause irreversible damage to fractures near the wellbore and affect conductivity.

BHP of Primary Depletion
The bottom hole pressure (BHP) controls the well production. The lower the BHP, the greater the production pressure difference, the higher the well production, and the faster the reservoir pressure drops. In shale reservoirs, multiple studies indicated that the production characteristics of tight reservoirs with high initial production and sharp decline in production will cause irreversible damage to fractures near the wellbore and affect conductivity. The reservoir pressure drop increases with the increase in the period of primary depletion. The incremental oil recovery also increases with the period of primary depletion.
If the gas huff-n-puff starts with a short period of primary depletion, it could result in a poor huff-n-puff effect. This is because the high reservoir pressure is not conducive to gas injection and migration to the deep part of the reservoir. At the same time, it will lead to a too high gas-oil ratio during the production stage.
The proposed period of primary depletion is five years, and the reservoir pressure should be above saturation pressure to prevent condensation in the wellbore which could affect the oil recovery.

BHP of Primary Depletion
The bottom hole pressure (BHP) controls the well production. The lower the BHP, the greater the production pressure difference, the higher the well production, and the faster the reservoir pressure drops. In shale reservoirs, multiple studies indicated that the production characteristics of tight reservoirs with high initial production and sharp decline in production will cause irreversible damage to fractures near the wellbore and affect conductivity.
We analyzed four cases of BHP of primary depletion, which are 35 MPa, 40 MPa, 45 MPa, and 50 MPa. The base case is 45 MPa. The effect of gas huff-n-puff under different BHP of primary depletion is examined. The results are summarized in Table 18, and Figures 28 and 29, it can be seen from the figures: ing primary depletion and the production period of gas huff-n-puff, the lower the BHP, the greater the pressure difference, and the higher the cumulative oil production. The BHP should be controlled within a reasonable range. High BHP is not conducive in improving oil recovery. If BHP is too low, it will cause gas precipitation and affect the liquid phase recovery; therefore, it is recommended that the reservoir pressure should be controlled within a suitable range above the saturation pressure.

Cumulative Oil Production for Primary Depletion
Prior to discussing the results of geological and engineering controls on gas huff-npuff, the sensitivity of the reservoir, completion, and depletion parameters on primary production should be discussed. Figure 30 is a Tornado plot that indicates the sensitivity of the 12 parameters on cumulative oil production for primary production. The sensitivity values are measured with reference to a base model which is based on the geological fea- The incremental oil recovery increases with the decrease of BHP. This is because during primary depletion and the production period of gas huff-n-puff, the lower the BHP, the greater the pressure difference, and the higher the cumulative oil production. The BHP should be controlled within a reasonable range. High BHP is not conducive in improving oil recovery. If BHP is too low, it will cause gas precipitation and affect the liquid phase recovery; therefore, it is recommended that the reservoir pressure should be controlled within a suitable range above the saturation pressure. Prior to discussing the results of geological and engineering controls on gas huff-npuff, the sensitivity of the reservoir, completion, and depletion parameters on primary production should be discussed. Figure 30 is a Tornado plot that indicates the sensitivity of the 12 parameters on cumulative oil production for primary production. The sensitivity values are measured with reference to a base model which is based on the geological features, hydraulic fracturing parameters and a good history match of primary production data. These values are obtained by setting each parameter to its minimum and maximum values one-at-a-time as shown in Table 6, while the rest of the parameters are kept at their base values. The green bars in Figure 30 correspond to the sensitivity measure (compared to the base case) when a selected parameter is set to its minimum and the red bars correspond to the sensitivity measure when a selected parameter is set to its maximum value. It can be seen from Figure 30 that reservoir thickness, BHP of primary depletion, and natural fracture permeability are recognized as the most sensitive parameters. It is followed by the period of primary depletion, half-length of hydraulic fractures, natural fracture density, TOC, fracture conductivity, and matrix permeability in order. Pore radius of shale matrix, type of SRV region, and magnification of permeability in the SRV region have the least effect on cumulative oil production for primary production. An implication of this is that these least sensitive parameters could have an insignificant impact on the gas huff-n-puff process.

Incremental Oil Recovery
Multiple scenarios of compositional simulations of huff-n-puff gas injection for the proposed 12 parameters have been conducted and the incremental oil recovery for each scenario was obtained.
To evaluate the relative impact of the parameters, a Tornado plot (Figure 31) was Figure 30. Tornado chart indicating the most and least sensitive parameters for cumulative oil production for primary depletion. Each green bar corresponds to a sensitivity measure (compared to a base case) when a parameter is set to its lower bound. The red bars are for the upper bounds. The vertical blue line denotes the base case.

Incremental Oil Recovery
Multiple scenarios of compositional simulations of huff-n-puff gas injection for the proposed 12 parameters have been conducted and the incremental oil recovery for each scenario was obtained.
To evaluate the relative impact of the parameters, a Tornado plot ( Figure 31) was created for sensitivity analysis on the incremental oil recovery. It is clear that fracture conductivity, natural fracture density, period of primary depletion and natural fracture permeability are the most sensitive parameters. It is then followed by reservoir thickness, TOC, BHP of primary depletion, matrix permeability, and half-length of hydraulic fractures in order. The effect of pore radius of shale matrix, type of SRV region and magnification of permeability in SRV region on incremental oil recovery is marginal. This rank is based on the range of each parameter listed in Table 6. The range for the percentage of recovery increase over the primary recovery for one huff-n-puff cycle (nearly 2300 days of production) is about 0.23-0.87%.

Gas Usage Ratio
The gas usage ratio is a key parameter to evaluate the efficiency of gas injection, which is an indicator of output-to-input ratio or economical efficiency. It is defined as the total volume of injected gas divided by total barrels of improved oil production by gas injection. For a given gas huff-n-puff scheme with a fixed composition of injected gas, a lower gas usage ratio usually implies higher economical efficiency. Figure 32 is a Tornado plot that indicates the sensitivity of the 12 parameters on gas usage ratio for the first huff-n-puff cycle. A direct implication of Figure 32 is that low fracture conductivity and short period of primary depletion could significantly increase the gas usage ratio and result in poor economical efficiency of the gas huff-n-puff process. Nevertheless, increasing the natural fracture density could reduce the gas usage ratio. Therefore, we believe natural fractures associated with secondary hydraulic fractures (natural fractures activated by hydraulic fracturing) are the most sensitive parameters for gas huff-n-puff in Duvernay shale oil reservoirs. This is because higher fracture density (both natural fractures and hydraulic fractures) provides a much larger matrix-surface area. It can be seen from Figure 31 that the incremental oil recovery is positive related to natural fracture density, fracture conductivity, matrix permeability, and period of primary depletion, while it is inversely related to reservoir thickness, TOC, BHP of primary depletion, natural fracture permeability, and half-length of hydraulic fractures.
It is indicated from Figure 31 that natural fracture is key to gas huff-n-puff effectiveness. The increase of natural fracture density from 0.1 fractures/m to 0.5 fractures/m has resulted in a 0.27% increase of incremental oil recovery, which is 45% higher than the base case value (0.6%). This is attributed to a larger amount of injected gas available to contact the matrix area and extract hydrocarbons through natural fractures.

Gas Usage Ratio
The gas usage ratio is a key parameter to evaluate the efficiency of gas injection, which is an indicator of output-to-input ratio or economical efficiency. It is defined as the total volume of injected gas divided by total barrels of improved oil production by gas injection. For a given gas huff-n-puff scheme with a fixed composition of injected gas, a lower gas usage ratio usually implies higher economical efficiency. Figure 32 is a Tornado plot that indicates the sensitivity of the 12 parameters on gas usage ratio for the first huff-n-puff cycle. A direct implication of Figure 32 is that low fracture conductivity and short period of primary depletion could significantly increase the gas usage ratio and result in poor economical efficiency of the gas huff-n-puff process. Nevertheless, increasing the natural fracture density could reduce the gas usage ratio. Therefore, we believe natural fractures associated with secondary hydraulic fractures (natural fractures activated by hydraulic fracturing) are the most sensitive parameters for gas huff-n-puff in Duvernay shale oil reservoirs. This is because higher fracture density (both natural fractures and hydraulic fractures) provides a much larger matrix-surface area.

Suggestions for Gas Huff-n-Puff Potential Area Screening
Based on compositional simulation of geological and engineering controls on gas huff-n-puff in Duvernay shale volatile oil reservoirs, screening criteria for gas huff-n-puff potential area are proposed. We believe that areas with higher natural fracture density, thicker reservoir, higher fracture conductivity, moderate fracture half-length, and moderate period of primary depletion are preferable for gas huff-n-puff in Duvernay shale volatile oil reservoirs to enhance oil recovery. Table 19 lists the screening criteria in the order of sensitivity for incremental oil recovery. Table 19. Screening criteria for gas huff-n-puff potential area in Duvernay shale volatile oil reservoirs.

Sensitivity Order
Parameter Screening Criteria 1 Fracture conductivity ≥100 mD·cm 2 Natural fracture density ≥0.1 fracture/m 3 Period of primary depletion ≥5 years Figure 32. Tornado chart indicating the most and least sensitive parameters for gas usage ratio. Each green bar corresponds to a sensitivity measure (compared to a base case) when a parameter is set to its lower bound. The red bars are for the upper bounds. The vertical blue line denotes the base case.

Suggestions for Gas Huff-n-Puff Potential Area Screening
Based on compositional simulation of geological and engineering controls on gas huff-n-puff in Duvernay shale volatile oil reservoirs, screening criteria for gas huff-n-puff potential area are proposed. We believe that areas with higher natural fracture density, thicker reservoir, higher fracture conductivity, moderate fracture half-length, and moderate period of primary depletion are preferable for gas huff-n-puff in Duvernay shale volatile oil reservoirs to enhance oil recovery. Table 19 lists the screening criteria in the order of sensitivity for incremental oil recovery.

Conclusions
In this study, we built a comprehensive numerical compositional model in combination with the EDFM method that simulates the huff-n-puff gas injection process in an actual Duvernay horizontal well based on a good history match. Multiple scenarios of compositional simulations of huff-n-puff gas injection for the proposed twelve parameters have been conducted and effects of reservoir, completion, and depletion development parameters on huff-n-puff in Duvernay Shale volatile oil reservoirs are evaluated. The sensitivity of the proposed twelve parameters on cumulative oil production for primary depletion, incremental oil recovery, and gas usage ratio were examined. Moreover, we proposed a set of screening criteria for gas huff-n-puff potential areas in Duvernay shale volatile oil reservoirs. The following conclusions can be made from this study: (1) A complex numerical model consisting of hydraulic fractures and complex natural fractures can be easily built using the EDFM method. (2) Base on the actual field production data and laboratory test results, a good history match has been achieved using a comprehensive numerical compositional model in combination with the EDFM method. (3) The range for the recovery increase over the primary recovery for one gas huff-n-puff cycle (nearly 2300 days of production) in the Duvernay shale volatile oil reservoir is about 0.23-0.87%. (4) The sensitivity of the reservoir, completion, and depletion parameters on incremental oil recovery by huff-n-puff can be ranked as: fracture conductivity, natural fracture density, period of primary depletion, natural fracture permeability, reservoir thickness, TOC, BHP of primary depletion, matrix permeability and half-length of hydraulic fractures. (5) Sensitivity analysis indicates that natural fractures associated with secondary hydraulic fractures are the most sensitive parameters for gas huff-n-puff in Duvernay shale oil reservoirs which can be attributed to that both natural fractures and hydraulic fractures provide much larger matrix-surface area during the gas huff-n-puff process. (6) Screening criteria for gas huff-n-puff potential area in Duvernay shale volatile oil reservoirs are proposed. Areas with higher natural fracture density, thicker reservoir, higher fracture conductivity, moderate fracture half-length, and moderate period of primary depletion are preferable for gas huff-n-puff EOR. (7) A good numerical tool (EDFM) is proposed in this study that can simulate the gas Huff-n-Puff process in the shale oil reservoir with multiple hydraulic fractures and complex natural fractures more effectively and faster.