Effects of Diffusion, Adsorption, and Hysteresis on Huff-n-Puff Performance in Ultratight Reservoirs with Different Fluid Types and Injection Gases

: Cyclic solvent injection, known as solvent huff-n-puff, is one of the promising techniques for enhancing oil recovery from shale reservoirs. This study investigates the huff-n-puff performance in ultratight shale reservoirs by conducting large-scale numerical simulations for a wide range of reservoir ﬂuid types (retrograde condensate, volatile oil, and black oil) and different injection gases (CO 2 , C 2 H 6 , and C 3 H 8 ). A dual-porosity compositional model is utilized to comprehensively evaluate the impact of multicomponent diffusion, adsorption, and hysteresis on the production performance of each reservoir ﬂuid and the retention capacity of the injection gases. The results show that the huff-n-puff process improves oil recovery by 4–6% when injected with 10% PV of gas. Huff-n-puff efﬁciency increases with decreasing gas-oil ratio (GOR). C 2 H 6 provides the highest recovery for the black oil and volatile oil systems, and CO 2 provides the highest recovery for retrograde condensate ﬂuid type. Diffusion and adsorption are essential mechanisms to be considered when modeling gas injection in shale reservoirs. However, the relative permeability hysteresis effect is not signiﬁcant. Diffusion impact increases with GOR, while adsorption impact decreases with increasing GOR. Oil density reduction caused by diffusion is observed more during the soaking period considering that the diffusion of the injected gas caused a low prediction error, while adsorption for the injected gas showed a noticeable error.


Introduction
Tight and ultratight (shale) hydrocarbon reservoirs have noticeably contributed to oil and condensate productions in the last decade, and it is projected to be doubled in the next decade [1]. The increase in oil production from tight formations is credited to long horizontal wells and the hydraulic fracturing technique. Despite higher production rates due to the advanced drilling techniques, fast decline production rates are always associated with producing tight formations due to the low permeability nature of these formations. This rapid decline in production rates motivated the implementation of different enhanced oil recovery (EOR) techniques to multistage hydraulic fractured horizontal wells in order to increase hydrocarbon production. Among various EOR techniques, solvent (gas) cyclic injection, known as huff-n-puff, has been one of the promising techniques for shale oil formations [2][3][4]. The cyclic gas injection has been shown to improve oil recovery factors (RF) by 33% to 85% based on the reservoir's characteristics and the operation parameters experimentally [5,6]. Note that RF is defined as the ratio of produced hydrocarbon to total hydrocarbon originally in the reservoir. The common solvents that have been considered include CO 2 , N 2 , CH 4 , C 2 H 6 , and C 3 H 8 . Due to the experimental success of cyclic injection, several field tests have been conducted. However, there is a gap between field recovery factors and the ones predicted from modeling [7]. Therefore, the validity of some common assumptions in shale reservoirs should be investigated, including the importance of diffusion, adsorption, and permeability hysteresis.
Introducing new fluid to a reservoir causes compositional changes due to component interphase transfer. Tracking these compositional changes and the associated fluid properties modifications due to the introduction of new fluid is crucial in accurately predicting fluid flow in porous media. The hydrocarbon recovery factor from shale reservoirs increases with higher Gas-Oil-Ratio (GOR): typically less than 10% for black oil reservoirs, 30-70% for gas reservoirs, and 10-30% for condensate [8][9][10][11]. Cronin et al. [11] investigated the impact of the reservoir fluids and injection gas and concluded that cyclic injection performed better for lower GOR than condensate. They highlighted the important role of density reduction and composition dilution in explaining the low recovery factor from ultralight oil reservoirs compared to gas reservoirs. Moreover, they concluded that the mass concentration reduction, which can be achieved by fluid expansion or composition dilution, is the main factor for improving the recovery factor from ultratight formations.
During the huff-n-puff process in ultratight formations, convection dominates the fluid flow in the fracture network, while diffusion is comparable to advection in the shale matrix due to the nanoscale pore sizes and nano-Darcy permeability [11][12][13]. Consequently, when studying multiphase, multicomponent fluid transport in fractured shale reservoirs, diffusion and advection have to be considered. Alharthy et al. [3] studied the impact of diffusion experimentally and proved that neglecting diffusion could result in a 20% underestimation in the oil recovery factor. Cronin et al. [12] developed a semi-analytical model to address the need for a compositional diffusion-based approach for simulating huff-n-puff in ultratight oil reservoirs.
Diffusion coefficients that define the molecular-level transport of fluids should be included in the modeling of cyclic gas injection processes. The diffusion coefficients are functions of pressure, temperature, density, and composition. For example, the CO 2 diffusion coefficient value had decreased by an average of 8% when CO 2 was injected into the black oil reservoir. This decrease in the diffusion coefficient value is more pronounced closer to the well and is less further away. Empirical coefficients are usually used to define the diffusion coefficient because of the complex process of measuring the diffusion coefficient. Different correlations have been introduced to determine the multicomponent diffusion coefficients, but the widely used ones in shales include [14][15][16]. Alharthy et al. [3] investigated the reliability of four correlations to determine the diffusion coefficients for shale oil reservoirs. They concluded that the Hayduk and Minhas correlation [16] agrees well with the experimental data. Du and Nojabaei [17] investigated different diffusion coefficient correlations and concluded that the Sigmund correlation [14] shows the closest match to the experimental data.
Shale reservoirs have a large storage capacity because of their high organic matter content [18], where hydrocarbon molecules can be adsorbed. Adsorption is a function of the adsorbent surface area, the effective molecular size of the sorbed gas, and the combined attraction energy between adsorbate and adsorbent. Due to the significant internal surface area in shale organic matters and the high attraction between the gas molecules and the organic matter in shales, the adsorption effect is substantial [19]. Langmuir adsorption isotherm has been commonly used for shales [20]. Adsorption of the light hydrocarbon components causes a reduction in total pressure, resists pressure buildup, and mitigates solvent injection [21]. Therefore, adsorption is an essential mechanism when modeling hydrocarbon production or gas injection in organic-rich shales in order to accurately predict the reserves, recovery factors, and reservoir storage capacity of the injected gas. While many studies investigated the adsorption impact on gas shale reservoirs, the studies of the effect of adsorption on the EOR/EGR in tight shales are limited [21][22][23]. Adsorption affects the volume of gas in the reservoir and the shale retention capacity of the injected solvent, which is important in assessing the storage capacity of injection gases. Neglecting adsorption when a solvent is injected into shale reservoirs may result in an overestimation of production [24,25]. It would also affect the evaluation of different solvent capacities relative to shale formations [23]. Wang et al. [22] investigated the impact of adsorption on CO 2 huff-n-puff in shales and concluded that oil recovery was overestimated by 13% when adsorption is neglected.
The hysteresis effect describes a property that relies on the current and historical values of another property. The hysteresis of gas relative permeability and capillary pressure should be considered because their values are dependent on fluid saturation history. While the gas relative permeability hysteresis cannot be neglected when gas is injected into oil reservoirs, the capillary pressure hysteresis can be ignored in tight shales because of the small capillary length ratio to the grid resolution [26][27][28].
This study investigates diffusion, adsorption, and hysteresis effects when different gases are injected into ultratight reservoirs with different fluid types (black oil to condensate). First, we investigate hydrocarbon primary production in the presence and absence of diffusion and adsorption. Then, we examine their impact during huff-n-puff for different types of fluid and injection gas. Finally, we compare the performance of each injection gas for different fluid types and make recommendations on the injection gas based on the reservoir fluid type.

Theory and Method
We selected the Eagle Ford shale (South Texas, USA) to investigate the impact of adsorption, diffusion, and relative permeability hysteresis on oil production for different reservoir fluids. Two main reasons for this selection are as follows: first, the Eagle Ford spans the entire PVT window from dry-gas, wet-gas, retrograde condensate, and volatile oil to black oil; second, the Eagle Ford is one of the most prolific/actively developed shale reservoirs in North America with an extensive literature.

Numerical Model
We utilized a compositional and unconventional simulator (GEM) by Computer Modelling Group (CMG) in order to solve spatial-temporal mass conservation equations and predict multiphase flow and transport phenomena in shale reservoirs. A dual-porosity model was used to include natural fractures in the model and couple fracture and matrix flow mechanisms. The heterogeneity of the Dykstra-Parsons coefficient was assumed to be less than 0.55 in order to minimize the impact of heterogeneity on production [28]. The spatial reservoir characteristics were assumed to be homogeneous and isotropic. Figure 1 (top) depicts the Cartesian dual-porosity model built to represent the matrix and natural fracture network. Figure 1 (bottom) shows the element of symmetry representing the reservoir model discretized by different grid block sizes. The model dimensions are 1880 ft in length, 300 ft in width, and 50 ft in height. A 1680 ft horizontal well was assumed to be drilled in the middle of the reservoir with ten hydraulic fractures. The hydraulic fractures are equally spaced at 180 ft, the fracture half-length is 170 ft, and the slab length is 210 ft. The grid blocks near the hydraulic fractures were locally refined logarithmically in I and J directions. The domain under consideration was split into 40 symmetric slabs because of its symmetric nature, where each slab represents a quarter fracture connected with its respective matrix volume. Table 1 provides reservoir and fracture properties, which are based on Eagle Ford shale. Figure 2 shows relative permeability curves for an Eagle Ford reservoir obtained from history matching [23]. The PVT data of three different fluid compositions from Eagle Ford representing low GOR black oil, high GOR volatile oil, and gas condensate with low condensate gas ratio (CGR) were obtained from Orangi et al. [29]. Table 2 provides the fluid composition of three reservoir fluids, and the corresponding P-T diagrams are shown in Figure 3. It can be readily observed from these data that the selected fluid types cover a large PVT window from retrograde condensate to black oil, with CH 4 mole fractions ranging from 31% to 64%.

Hysteresis
We used the Killough method to include the relative permeability hysteresis effect [30]. During production, liquid saturation increases while gas saturation decreases. When the gas saturation proceeds below residual saturation, the injected gas becomes trapped and immobile. When gas injection starts, gas saturation increases and liquid saturation decreases, and gas mobility is a function of the history of gas saturation. Accordingly, oil-water and gas-liquid relative permeability scanning curves are calculated by using the following: where k im rg is the bounding drainage curve, k d rg is the bounding drainage curve, S gi is the gas saturation when the process is switched, S gi,max is the maximum trapped gas saturation, S gt is the trapped gas saturation, and S gt,max is the maximum trapped gas saturation. S gt,max = 0.4 was used, which is the highest reported value for carbonate reservoirs [28].

Diffusion
We selected Sigmund correlation [14] to predict the spatial-temporal multicomponent diffusion coefficients for all the components. The empirical correlation for predicting the diffusion coefficient of component i in the mixture is as follows: where D * i is the diffusion coefficient of component i in the mixture, D ij is the binary diffusion coefficient, v ci is the component critical molar volume, ρ is the molar density, R is the gas universal constant, and y i is the mole fraction. The collision diameter Ω ij and the collision integral σ ij of the Lennard-Jones potential are related to the critical component properties [31]. The diffusion coefficients are then corrected with D e f f = D * i /τ to account for the effects of tortuosity and porosity, where τ k = Fφ is the tortuosity, and F is the resistivity factor.

Adsorption
We employed the extended Langmuir isotherm for multicomponent to define the competitive adsorption behavior in the model [32,33]: where B i [1/psi] is the Langmuir isotherm pressure constant, ω i [SCF/ton] is the adsorbed volume per unit mass of rock, ω i,max is Langmuir volume constant, P is the pressure, and y i is the molar fraction of the adsorbed component i in the gas phase. Table 3 presents the Langmuir adsorption parameters for light hydrocarbon components and injection gases that were used in this study [34]. We only considered the adsorption of the light hydrocarbon components and gases as they have significant sorption potential.

Gas Injection Protocol
We considered miscible solvent injection by maintaining the reservoir pressures above the minimum miscible pressure (MMP) at all stages. Before starting the solvent injection process, the reservoir produces by natural depletion for two years, known as the primary production period. The primary depletion period is followed by five injection cycles. A total injected solvent volume comprising 10% of the reservoir pore volume is injected during the injection cycles. The total injected solvent is split equally throughout the five stages. Each injection cycle length is six months divided into three months injection followed by a three-month soaking period. Each injection cycle is followed by two years of the production period.

Grid Size Sensitivity Analysis
The grid size in the numerical model should be sufficiently small for capturing the molecular transfer between the injected gas and the reservoir fluid. However, such a model will be extensive due to the ultra-small pore sizes in shales and would be associated with numerical difficulties. Therefore, one should choose a grid block size that does not associate with a high numerical error. We tested different grid sizes with a Logarithmic refinement (LRF) gridding system in the stimulated area and compared it to a 1 × 1 ft 2 grid block size in order to optimize the grid block size. The considered grid block sizes are 5 × 5, 10 × 10, and 20 × 20 ft 2 with LRF. Figure 4 shows the oil recovery factor (RF) when different grid block sizes are used, revealing a slight difference in oil recovery prediction of 0.01%. The 20 × 20 ft 2 provides a slightly higher recovery factor, and the recovery factor decreases with the number of the grid block size. The grid block associated error is consistent through primary depletion and the huff-n-puff process. While there is no significant difference in oil recoveries for different gridding systems, a noticeable difference in the molar distribution of the injected solvent is pronounced in Figure 1a-d. This slight difference in the performance is expected because LRF maintained the gridding sizes to 1 × 1 ft 2 around the fracture location where most of the compositional changes occur. This slight difference in prediction is tolerated in favor of lower simulation time when larger grids with LRF are used. The simulation time is less in the LRF case due to the lower number of grid blocks (56,700 grid blocks for the 1 × 1 ft 2 model vs. 2328 for the 5 × 5 ft 2 LRF model). Using the 1 × 1 ft 2 gridding system required 87 times more CPU time compared to the 5 × 5 ft 2 gridding system. Therefore, we decided to use the 5 × 5 ft 2 model with the LRF gridding system.  Figure 5 shows the desorption effect on the primary oil recovery factor for the three reservoir fluids. The adsorption impact was noticeable because of the light component desorption from the organic matter. Diffusion and hysteresis were not present due to the absence of the concentration gradient and gas phase. The results show that neglecting the adsorption during the primary depletion underestimates oil recovery by 1.69%, 1.06%, and 1.16% for black oil, volatile oil, and condensate. While the underestimation percentages are close, the associated error for each fluid type is different. The associated errors were calculated by using Error = RF With Adsorption − RF Without Adsorption /RF With Adsorption , which are 31%, 13%, and 8% for black oil, volatile oil, and condensate, respectively. The associated error values show that the adsorption impact is higher for black oil reservoirs and decreases with an increase in GOR.  Figure 6 presents oil recovery for different injection gases and fluid types, where adsorption, diffusion, and hysteresis are included. The huff-n-puff process improved the oil recovery factor by 4.07% to 6.2% in black oil, volatile oil, and condensate reservoirs. While the increase in the recovery factor for the different fluids is relatively close, we can see a significant difference in the effectiveness of the injection in different reservoir fluids by comparing the percentage increase to primary production. Therefore, we evaluated the impact of the injection on each fluid type based on the efficiency of the injection using E f f iciency = RF Hu f f −n−pu f f − RF Primary /RF Primary . Accordingly, C 2 H 6 performed best for black and volatile oil by increasing the recovery factor by 5.98% and 6.02%, which shows an efficiency equivalent to 116.34% and 44.6%, respectively. CO 2 performed best for condensate reservoirs by increasing the recovery factor by 4.67%, which means 19.61% efficiency. A comparison of injection efficiency suggests that cyclic injection is more efficient in black oil reservoirs and less efficient in condensate reservoirs, i.e., cyclic injection efficiency is negatively correlated to the fluid light component content and GOR.

Diffusion and Adsorption Effects
When a solvent is injected into the reservoir, adsorption increases the solvent mobility to the matrix. Hence, adsorption increases the contact between the solvent and the oil. During the soaking period, diffusion facilitates the mobility of the light and intermediate molecules from the oil to the solvent, decreasing oil density. Hence, diffusion increases flowing light and intermediate components. Since intermediate components are characterized by higher adsorption/desorption, they depart the matrix faster. Hence diffusion increases the mobility of intermediate components during the soaking period, and adsorption increases the mobility during the depletion process. This process explains the increase in oil recovery when diffusion is included. Figure 7 shows the oil recovery profiles for CO 2 , C 2 H 6 , and C 3 H 8 , respectively, where the impact of the adsorption and diffusion is observed. Adsorption and diffusion provided higher production, credited to adsorption and enhanced by diffusion later as the diffusion impact became more pronounced after the second injection cycle. The recovery factor decreases at a later production time, which overturns the effect of neglecting adsorption and diffusion from underestimating to overestimating the recovery factor. This behavior is observed in C 3 H 8 injection but not CO 2 and C 2 H 6 because it happens at a later time beyond the studied period. The diffusion impact is highest for the condensate reservoirs because of the higher content of the light and intermediate components. Diffusion impact increases with GOR while adsorption impact decreases.
In order to better understand the impact of diffusion, we investigated the changes in oil density. Since the diffusion mechanism controls the molecular transfer from higher to lower concentrations, it enhances oil mobility by reducing oil density. Figure 8 shows the average oil density changes in the reservoir for different fluid types during the first two injection cycles. Since the concentration gradient is the driving force for diffusion, its impact is expected to be more presented during the soaking period. At the end of the injection period, a large concentration gradient exists between the solvent in the natural fractures and the oil in the matrix. Since the convective dispersion force is small during the soaking period, the molecular exchange is dominated by diffusion. As depicted in Figure 8, the reduction in oil density due to molecular transfer is evident during the soaking period. Oil density reduction during the soaking period is continuous and more efficient in each cycle, explaining the increase in diffusion impact with more cycles. These observations conclude that more extended soaking periods are expected to reduce oil density, which will increase oil recovery and the efficiency of the cyclic injection.

Retention Capacity
The available molecules in the reservoir compete for the same storage space in the formation. Hence, the retention capacity of the injected gases is sensitive to the changes in the fluid composition that occur during production. The retention capacity is evaluated by comparing the produced volume of the solvent relative to the injected volume. In this study, we investigated the retention capacity for CO 2 in order to assess potential CO 2 sequestration in ultratight formation. Figure 9 shows that black oil reservoirs provide higher retention capacity for the CO 2 followed by volatile oil and condensate reservoirs. The higher retention of the black oil is due to the lower initial adsorbed components content, leaving more storage capacity for the injected gas. Accordingly, the retention capacity increases with lower initial content of the light and intermediate components, especially the initial content of the injected solvent. Note that retention capacity is calculated as the percentage of injected fluid retained in the reservoir at the end of the production period of 15 years using Retention Capacity

Diffusion and Adsorption Effects
We investigated the impact of adsorption and diffusion on retention capacities, and the first observation was the correlation between the oil recovery factor and the retention capacity. As shown in Figure 10, since solvent injection improves oil recovery through molecular exchange, the higher the retention, the better the oil recovery performance. In volatile oil and condensate reservoirs, diffusion improved solvent retention as it increased the retention by 5% and 5.4% beyond the adsorption only. These results show that both diffusion and adsorption should be considered when retention capacity is considered. Figure 11 shows the consequence of using constant diffusion coefficients for each component computed using the initial composition. Hence, assuming constant diffusion coefficient values may result in an incorrect prediction. This observation also highlights the sensitivity of the results for diffusion coefficient values.

Multicomponent Diffusion Effect
The diffusion coefficient decreases with heavier hydrocarbon components and the diffusion of all the components should be considered when compositional changes occur in the reservoir. However, the inclusion of multicomponent diffusion increases computational time. Hence, an evaluation of the intermediate and heavy component diffusion impact is needed to define the associated error if the diffusion of the injected gas only is considered. Figure 12 presents the oil recovery factor for different reservoir fluids when CO 2 is injected, considering only CO 2 diffusion versus multicomponent diffusion. The reservoir modeled with the diffusion of injected gas-only caused an overestimation of 0.48%, 0.47%, and 0.01% for the black oil, volatile oil, and condensate. The impact of neglecting the multicomponent diffusion increases with decreasing GOR. With this low associated error, modeling huff-npuff injection with injected gas diffusion is accepted to reduce computational cost.    Figure 13 shows the impact of the adsorption of each component when CO 2 is injected into different reservoir fluids. The adsorption of the intermediate components has a lower impact on the recovery factor than the light components (CO 2 and Methane). Hence, the adsorption of the light components has to be included to model cyclic injection in shales adequately. Figure 14 reveals the hysteresis effect on oil recovery when CO 2 is injected into the different reservoirs. The increase in gas saturation close to the hydraulic fracture was not enough to impact production. This may occur due to the low average gas saturation through the process, as the average reservoir gas saturation was always below critical gas saturation. The hysteresis impact was expected to be more pronounced in the condensation because of the formed gas in the reservoir. However, the formed gas was trapped, but its effect was not enough to show differences.

Conclusions
Using a dual-porosity compositional model, we investigated the effects of diffusion, sorption, and hysteresis on huff-n-puff (10% of the pore volume) performance in shale reservoirs (e.g., Eagle Ford) by considering different fluid types (i.e., retrograde condensate, volatile oil, and black oil) and injection gases (i.e., CO 2 , C 2 H 6 , and C 3 H 8 ). The key conclusions are summarized as follows: • Huff-n-puff gas injection can potentially increase the oil recovery factor by 4.07% to 6.2% over 15 years. The increase in oil recovery is a function of injection gas and reservoir fluid types. Injection efficiency is higher for the reservoirs with lower GOR due to the lower primary production. Furthermore, CO 2 provides the highest recovery from condensate, and C 2 H 6 provides the highest recovery from black oil and volatile oil reservoirs. Hence, there is considerable potential for cyclic gas injection to overcome the rapid production decline in ultratight oil and condensate reservoirs.

•
Diffusion and adsorption mechanisms should be included in modelling EOR processes in shale reservoirs as they have substantial impacts on oil production. However, the relative permeability hysteresis impact was insignificant for the miscible injection. • Diffusion coefficient is a spatial-temporal function that has to be updated based on the changes in pressure, composition, and density. Assuming that the initial diffusion coefficient values are constant throughout the process may result in overestimating production. While diffusion decreases oil density during the soaking period, both diffusion and adsorption show higher density reduction. The captured density reduction trend due to adsorption and diffusion suggests that more extended soaking periods would increase oil mobility and production. • Desorption improves primary depletion oil recovery by 31%, 13%, and 8% in black oil, volatile oil, and condensate, respectively. The gas retention capacity ranges from 35 to 45%. The black oil reservoirs showed higher retention capacity, while condensate and volatile oil reservoirs showed the lowest. Retention capacity decreases with increasing GOR, neglecting adsorption and diffusion results in lower retention capacity.