Numerical modeling of molecular diffusion and convection effects during gas injection into naturally fractured oil reservoirs

. Gas injection into a naturally fractured oil reservoir keeps the reservoir pressure and increments the initial recovery from the reservoir. The main aim of this work was to develop a numerical model to calculate the mass transfer (molecular diffusion and convection) between a gas injected in the fracture and residual ﬂ uid (gas and oil) in a matrix block. The dual continuum model is applied to describe ﬂ ow behaviour and ﬂ uid recovery in porous media. Finally, the model is validated by comparing the outcomes with the results of two experimental works available in the literature. The mathematical model results are in agreement with the laboratory data including recovery of each component, saturation pro ﬁ le, and the pressure gradient between matrix and fracture. Modeling results show that after 25 days of N 2 injection, the lighter and heavier components (C1 and C5) are recovered about 51% and 39%, respectively. These amounts for CO 2 injection are 49% and 27%. It is found that the convection mechanism has a great effect on preventing the pressure drop of the reservoir during injection operations. In the nitrogen injection, without considering the convection, after 30 days, the matrix pressure reaches 1320 Psi from 1479 Psi but after 30 days, considering the convection, the pressure reaches 1473 Psi from 1479 Psi.

Overall composition of c

Introduction
The major parameter for production from heavy oil reservoirs is high oil viscosity [1][2][3]. Gas injection into a naturally fractured oil reservoir not only maintains the reservoir pressure but also increases the reservoir initial recovery [4][5][6][7][8].
When gas is injected into an unsaturated fractured reservoir, the gas enters both fracture and matrix media [9,10]. It causes oil to swell and move from the matrix [11][12][13][14][15]. Furthermore, injected gas in the matrix causes some components of oil to evaporate and move to the gas phase. Therefore, oil is transported by both convection and diffusion mechanisms. Moreover, increasing the mobility of heavy oils by reducing their viscosity can be improved by injecting diluents such as gases. A lot of studies about injecting diluent gas have been carried out [11,[16][17][18].
The studies reveal that the diffusion of CO 2 into the core accounts oil to swell and subsequently reduce the viscosity which results in a higher gravity drainage rate [19,20]. Also, dissolving gas in oil causes the oil to expand and reduce its density, and thereby the oil exits from the rock. Moreover, when liquid and gas phases are not in equilibrium, one or more components move from one phase to another one due to molecular diffusion [4,11,[21][22][23][24]. The presence of dual porosity including porous media (low permeability/ high porosity) and fracture networks (high permeability/ low porosity) with different physical properties causes the structure of fractured reservoirs much more complex than conventional non-fractured reservoirs so that conventional methods which have been used for studying the performance of the reservoirs cannot present appropriate results to examine the behavior of reservoirs [10,14,25]. The molecular diffusion has been presented by Ertekin et al. [26]. Multiple mechanisms of fractured gas reservoirs have been considered by Ayala et al. [27], and they concluded that, unlike conventional reservoirs, molecular diffusion in fractured reservoirs has an important role in gas injection into oil reservoirs and gas production from gas condensate reservoirs. Moreover, secondary recovery and tertiary recovery have been studied by Karimaie [28]. According to the study, diffusion has a strong effect on secondary and tertiary oil recoveries. In addition, a compositional model has been presented by Kazemi and Jamialahmadi [29]. In their work, it has been concluded that at an early stage, swelling of oil has occurred, and gravity drainage has been followed by a slow extraction mechanism that can recover the heavy and intermediate components.
Diffusion coefficients of gas and oil are very important for any field simulation of a fractured reservoir. Moreover, CO 2 injection tests in systems of matrix/fracture by using a finely-gridded compositional simulator have been studied by Alavian [30], and a single matrix block has been represented to understand the impact of molecular diffusion on oil recovery by CO 2 injection in a fracture-matrix system. The outcome of his work shows that diffusion would have a crucial role in oil recovery, except at very low CO 2 injection rates. Hoteit [31] demonstrated analytically and numerically that classical Fick's law does not consider the total flux balance and is unable to take the correct direction of the diffusion flux in some situations. Moreover, a numerical model based on Fick's law has been suggested by him to predict the gas-oil transfer mechanism. In the suggested model, the gradient of chemical potential is considered as the driving force. Also, the implementation of the simulations that consider diffusion using various solution methods has been examined by Mohebbinia and Wong [32]. Also, optimal performance is achieved by the partially implicit approach. The simulation outcomes have revealed that constant diffusion coefficients may not predict an acceptable result during gas flooding and oil recovery might be overestimated or underestimated. It has been observed by Darvish et al. [33] for a higher porosity case that, the viscosity reduces later and the effect of the molecular diffusion on oil recovery will increase if the porosity of the matrix is big. Jin et al. [34] examined the CO 2 EOR process in tight reservoirs and concluded that it is preferred to inject CO 2 into the fractures. When CO 2 dissolves in oil, oil viscosity reduces and oil drains into the fractures. Experimental and field studies have shown that CO 2 has an important role in displacing oil in porous media via interfacial tension change, viscosity reduction, and subsequent reaction of carbonic acid with minerals as well as oil swelling [35,36]. For many conventional and unconventional reservoirs including heavy oil extraction in VAPEX (VAPour EXtraction) process [37,38] or oil extraction by solution-gas-drive [39], diffusion is an important parameter [4,17,29,[40][41][42][43][44][45][46].
Diffusion could be an important recovery mechanism in naturally fractured reservoirs while there are very few attempts to model diffusion and convection between a flowing gas through the fracture and oil and gas in the matrix. The first objective of this work was to mathematically model mass transfer between a gas phase in the fracture and oil and gas phases in the matrix in a single porosity model. There are several commercial models for modeling naturally fractured reservoirs. These models are classified as follows: Dual porosity model, dual porosity/dual permeability model, and dual continuum model. The dual continuum model is the model used for gas injection into naturally fractured reservoirs in this paper. The dual continuum approach is either the single-block model or the single-porosity model. This model will be examined in the following and finally, the proposed model will be validated with experimental data available in the literature.

Dual continuum theory
This theory is applied to model fractured reservoirs by considering a porous media as a single matrix block with an adjacent fracture. The fracture is considered a boundary condition for the matrix. This model is a good expression of the natural fracture reservoir because it considers the fluid flow between the matrix block and the fracture. Figure 1 shows an overview of this model. In this model, the matrix is separated but the fracture is not separated because the fracture acts as a boundary condition for the matrix.
Compositional simulators were created to study the oil recovery of a reservoir during the injection of gas. In the compositional simulation, it is supposed that the hydrocarbon phase and water phase are insoluble. So, the mass conservation equations are applied separately for hydrocarbon and water compounds. For the multi-phase compositional flow, viscosity, gravity, and capillary forces must be calculated. Moreover, if gas is injected, a diffusion mechanism should be taken into account to determine the mass transfer between phases. The equations governing the multiphase flow are presented below. In principle, these equations use the continuum equations for each component

Convection mechanism
The mass transfer due to convection is a mechanism that components are displaced by the motion of the fluid mass. The potential gradient is the driving force for convection (bulk flow). Inlet molar rate (mole/time) minus the exit molar rate of component c in three directions by convection transfer in the gas and oil phase as below:

See the equation (2) bottom of the page
Using the mathematical form of Darcy's law:

Molecular diffusion mechanism
Molecular diffusion is a mechanism that components are transferred by the random motion of molecules. Bird et al. [47] proved that for ideal or near ideal mixtures, gradients of concentration can be used instead of chemical potential gradients as the driving force. In this model, concentration gradient is considered as the driving force of molecular diffusion. If N c , o and N c , g (molar flux of Þ z ÁxÁy þ q g y c y g;z À Á z ÁxÁy À q 0 r c v 0;z ð Þ j zþÁz ÁxÁy À q g y c y g;z À Á zþy z ÁrÁy:

Flow equations
The flow equations could be acquired by substituting equations (2)-(6) in equation (1). Then, the obtained equations are divided by the volume of the element Dx Dy Dz and applying limit, we will have: Or, in vector notation, where q D, fm and q C, fm are the diffusion and convection mass transfer between the matrix and the fracture at the fracture-matrix boundary based on mol/s, respectively. The sum of the fluid moles in the volume element at any time is: Dx Dy Dz (q o S o + q g S g )/.
And the mole of component c is equal to: Therefore, the mole accumulation rate of component c is: By substituting m o and m g from equation (3) and N c,o and N c,g from equations (5) and (6) into equation (8), it becomes: Equation (9) is the general state of the composite multiphase flow inside the porous medium for each component in the gas and oil phases in the porous medium. The first bracket shows the convection mechanism in the gas and oil phases. The mechanism of diffusion is shown in the second bracket. q D, fm and q C, fm are the mass transfer of diffusion and convection between the matrix and the fracture. A mass balance equation only expresses the movement of water by the convection mechanism because we assume that the hydrocarbon phase is insoluble in the water phase. Thus, for the water phase, we have: The equations governing multi-phase compositional flow are derived from three sources: 1. The equilibrium equations of each component (Eq. (9)) and for water phase (Eq. (10)). 2. The phase equilibrium between the hydrocarbon phases is indicated as the below: 3. The auxiliary equations are the sum of the phase's saturation and the sum of the molar components in each phase is unity. Moreover, the oil, gas and water pressure should be related to the capillary pressure: The equations that governing multiphase compositional flow in the porous media are expressed by equations (9)- (14). This system of equations contains a set of (2n c + 6) equations with the same number of unknowns as the below: . . . ; x nc ; y 1 ; y 2 ; . . . ; y nc À Á :

Initial and boundary conditions
The initial and boundary conditions are presented next.

Initial conditions
It is supposed that there is gravity equilibrium in the model at time equal to zero. Also, pressure and composition at a reference are known. Since there is gravity equilibrium at time equal to zero, convective flow vanishes. Therefore, from Darcy's law: For a horizontal plane rD = 0, so equation (15) simplifies to: Equations (16) and (17) show that pressure, composition, and saturation are constant in a horizontal plane at time equal to zero. For the vertical direction, equation (15) becomes: Equation (18) means that vertical pressure distribution is given by the column weight.
Integrating equation (18) results in: where " c p is the average specific weight of phase p between z and z ref heights. If pressure at a reference height is given, then pressure at any point in the model can be calculated from equation (19).

Matrix boundary conditions
There are two boundary condition types for the matrix: 1-matrix sealed boundary condition and 2-matrix-fracture boundary condition (Fig. 1). These boundary conditions are attended as follows:

Matrix sealed boundary conditions
The total mass flux for all components in all phases vanishes at these boundaries. That is: Convection flux at the boundary: Diffusion flux at the boundary:

Matrix-fracture boundary conditions
The continuity equation in the fracture includes mass transport by diffusion and convection mechanisms in a laminar flow regime. The diffusion mass transfer rate at matrix-fracture surface is found by: where q g is gas stream density in the fracture. Convection between a matrix grid cell and the adjacent fracture (q C,fm,c ) is defined in the model based on Darcy's law as: See the equation (23) bottom of the page

Numerical Solution of equations
The numerical model of the resulting equation has been discretized using the IMplicit Pressure Explicit Saturation (IMPES) approach. We obtain the pressure implicitly and the saturation explicitly. Some of the equations in the multiphase compositional flow in the porous media are nonlinear. Therefore, using numerical methods, we substitute finite difference approximations in a nonlinear equation system for all derivatives. The equations are then linearized and solved using MATLAB software.

Validation of model
The mathematical model has been validated using laboratory data of Morel et al. [48], Le Romancer and Fernandes [49]. N 2 was injected in the first experiment and in the second experiment, CO 2 was injected into the porous media The matrix residual oil is C 1 and C 5 mixture. In the N 2 injection experiment, the oil and vapour phases were present initially, but in the experiment of CO 2 injection, the matrix was saturated with the oil phase. Experiments were carried out in a one-dimensional system to simulate the mass transfer between a gas in a fracture (N 2 or CO 2 ) and a residual oil (C 1 -C 5 ). The injected gas (N 2 or CO 2 ) diffuses into the residual oil in the porous matrix and evaporates the oil by convection and diffusion. The diffusion process alters the composition of the core fluid, resulting in changes in the phase's properties such as density. Recovery of each component (C 1 -C 5 ) during the matrix was obtained in both experiments.

Simulation of N 2 injection test
First, explanation of the N 2 injection experiment is presented followed by the results of mathematical model. Tables 1 and 2 show the model parameters and the relative permeability and capillary pressure, respectively. N 2 injection experiments were performed in a one-dimensional horizontal core by Morel et al. [48]. Figure 2 demonstrates the setup of the experiment. All parts of the core were insulated except one side where N 2 gas enters from. The sample of the Paris Basin Chalk was used as the core. N 2 was injected at an initial pressure of 1479 Psi, initial temperature of 38.5°C and the composition (C 1 (52.4 mole%) -C 5 (47.6 mole%)) were assumed constant. At the start of the experiment, the porous matrix containing a combination of C 1 and C 5 is distributed between the equilibrium liquid and gas phases. A gas with an initial saturation of 25% is presented within the matrix. Injected N 2 diffuses into the gas phase and solves at the boundary between the fracture and the matrix in the liquid phase. The pressure at the fracture-core boundary is fixed throughout the simulation and compositions are uniform along the core. The core is simulated with 20 grids in the x direction. The capillary pressure of gas-oil is corrected to a reference surface tension.
Nitrogen injection experiment was performed for 16 days and the simulation of the experiment was intended for 30 days. To better understand the importance of molecular diffusion in recovery, testing and simulation have been done in two ways: A -Regardless of the convection at the boundary of matrix-fracture (or q C, fm,c = 0). B -Considering the convection at the boundary of matrix-fracture.
The magnitude of the mechanisms of diffusion and convection in transfer of each component into the core is depicted in Figures 3a-3d. The positive sign in these figures means the mass transfer from the boundary of fracture-core to the core and the negative sign means the transfer of mass from the matrix to the boundary of fracture-core. The mass transfer of N 2 from the fracture to the core causes a countercurrent flow to the core. Figures 3a-3d reveal that N 2 is transferred to the end of the core by the gas convection and molecular diffusion from the boundary of fracture-core (positive values). The counter-current flow of the oil, from the end of the core to the fracture-core boundary, makes the N 2 to move to the fracture boundary by convection (negative values). Figures 3a and 3b depict that at the beginning of simulation, no matter whether the matrix-fracture boundary conditions are without convection (state A) or with convection (state B), N 2 is mainly transported in the gas phase into the core by molecular diffusion. Figure 3c describes if the convection between the matrix and the fracture is ignored (state A), then the molecular diffusion of the gas will stay as the dominant mechanism of N 2 transfer into the core for 28 days. But if the convection between the matrix and the fracture is considered (state B), then gas convection will be the dominant mechanism for mass transfer in 28 days (Fig. 3d).
Figures 4a-4d measure the value of convection and diffusion mechanisms of oil and gas in transfer of C 1 into the core. C 1 reaches to the boundary of fracture-core through the gas and oil phases and transported to the oil phase by convection (negative values). The gas flow that is moving from the boundary of fracture-core to the end of the core, moves methane in the reverse direction. Figures 4a and  4b show that for both matrix-fracture boundary conditions, without convection (state A) and with convection (state B), at the beginning of test times (t = 8 days), C 1 is mainly transmitted by molecular diffusion into the core because convection in the oil and gas phases moves C 1 in the reverse direction and finally neutralizes each other. If the movement between the matrix and the fracture is ignored (state A) the molecular diffusion of the gas stays the dominant mechanism of transport (Fig. 4c). However, Figure 4d reveals that at the end of times, gas convection is the most  important factor for transferring of C 1 from boundary to the core with considering convection between the matrix and the fracture in state B. Figures 5a-5d confirm that oil convection, disregarding to convection (state A) or considering convection (state B) between the matrix and the fracture, is the predominant factor of C 5 movement from the end of the core to the boundary of fracture-core. It is also the significant mechanism of pentane transport into the core for the whole simulation time.
It is obvious that C 1 recovery is greater than C 5 (because C 1 is the lighter component and also more mobile than the C 5 , which is a heavier component), as shown in Figure 6. Over time, the percentage of components recovery increases. For example, for component C 1 in state A, the recovery percentage has risen from 20% on the 5th day to about 43% on the 20th day. Recovery of both components C 1 and C 5 in state A (without convection) and in state B (with convection) shows that diffusion compared with convention has a crucial role in recovery of oil components. At the beginning of the simulation, the diffusion mechanism is the main recovery mechanism, but over time, the role of the convection mechanism in C1 recovery increases and this is the reason that why after t = 25 days, the boundary conditions of A and B no longer match together. Figures 7a-7c show saturated changes in 4, 8, and 16 days after starting the process. The differences between the modeling and experimental results are acceptable. At the beginning of the injection operation, at t = 4 days, the saturation of the gas inside the core is high and nitrogen has not yet been able to completely replace the gas inside the core. As the injection time increases, more nitrogen diffuses into the pores of the core and replaces the gas inside the core, reducing the percentage of gas saturation. At t = 4 days, the maximum saturation percentage for state A is 65%, while at t = 8 days, this value reaches to 55%. Also, at the end of the core, the percentage of gas saturation is high and less nitrogen has diffused to the end of the core. At t = 16 days,  the rate of gas saturation is 45% which is lower than in previous times. As we move to the end of the core, the percentage of gas saturation is higher as nitrogen has not been able to move to the end of the core. As nitrogen travels through the core due to pores inside the core, gases move to the end of the core and accumulate there as a result the percentage of gas saturation at the end of the core increases. The main    0.89 Oil initial total mole 1.30E-03 C 1 initial total mole 3.65E-04 C 5 initial total mole 9.37E-04 reason for the reduction in reservoir pressure is the release of oil-soluble gas. The pressure created by the gravity drainage mechanism keeps soluble the gas in the oil and maintains the pressure of reservoir. Figures 8a and 8b show the pressure of oil changes at the matrix-fracture boundary in the core for two states without and with convection, (A) and (B) respectively. Figure 8a demonstrates a significant drop in oil pressure from 1479 Psi at t = 0 day to 1320 Psi at t = 30 days when convection is ignored at boundary between matrix and fracture, While Figure 8b demonstrates that the convection of N 2 from the fracture into the matrix results the pressure of the core to stay constant around 1473 Psi for t = 30 days. Therefore, the presence of a convection mechanism prevents oil pressure to drop. Comparing Figures 8a and 8b, in the presence  Calculated and experimental C 1 and C 5 recoveries in CO 2 injection experiment. Fig. 11. Differential pressure between matrix and fracture in CO 2 injection experiment.
of convection (state B), the pressure drop at t = 30 days has finally reached from 1479 Psi to 1473 Psi, but without convection (state A) the pressure drop after t = 30 days, from 1479 Psi to 1320 Psi, it shows the great importance of the convection mechanism in preventing reservoir pressure drop.

Simulation of CO 2 injection experiments
Next, the mathematical model is investigated for CO 2 injection into the fractured reservoir and results are compared with experimental data carried out by Le Romancer and Fernandes [49]. In this case also, all the core surfaces except one is insulated. The Paris Basin Chalk was used as the matrix. CO 2 was injected into the fracture at a temperature of 38.5°C and pressure of 913.74 Psi. The matrix was saturated with a liquid phase of 28 mole% of C 1 and 72 mole% of C 5 in this test. Initially, the core has no gas saturation and there is an aqueous phase residing at 11% saturation in the core and the initial saturation of the gas is 0%, which is assumed to be constant throughout the core. The core is simulated in the x direction by 10 grids and at the fracturecore boundary. Model inputs for simulating of CO 2 injection experiments are given in Table 3. The specifications of the core in this experiment are presented in Table 4. Figures 9a-9c show the mass transfer of CO 2 , C 1 , and C 5 at the matrix-fracture boundary by the convection and gas diffusion mechanisms. Figure 9a proves that approximately 40 days after the start of the experiment, the CO 2 is transferred to the matrix-fracture boundary predominantly by diffusion. After 40 days, CO 2 leaves the core by the convection of gas phase and entering the fracture-core boundary by diffusion mechanism. The difference between the rate of gas convection and the CO 2 diffusion at the boundary of fracture-core increases with time. Figure 9b indicates that C 1 at the boundary of fracture-core is recovered mostly by diffusion about 40 days after the beginning of the experiment, while the gas convection begins to transfer C 1 from the boundary of fracture-core to the fracture after 40 days. In C 1 recovery after 70 days, mass transfer by gas convection and diffusion at the fracture-core boundary has the same importance. Figure 9c depicts that diffusion at the surface of fracture-core is the most substantial mechanism for C 5 recovery for up to 40 days. While the convection of gas from fracture-core to the fracture starts to participate in C 5 recovery after 40 days. Although recovery of C 5 with gas convection increases over time, C 5 transport by diffusion will remain as a significant mechanism. Like N 2 injection experiment, C 1 recovery is higher than the C 5 recovery (Fig. 10) and over time, the percentage of components recovery increases. Figure 11 proves the difference pressure between the fracture-matrix boundary. Due to the difference in the permeability of the fracture and the core, the pressure difference increases sharply and then remains at a constant value. The agreement between the experimental data and the obtained results from the model is acceptable. The pressure in the core drops sharply in the lack of convection between the fracture and the matrix. Figures 12a-12d compare saturated changes in 53, 67, 88, and 95 days. The more time passes from the start of injection, the oil saturates in the farther distance from the beginning of the core. For example, at t = 53 days, the oil reaches maximum saturation in x = 0.22 m from the beginning of the core. At t = 88 days, this distance is equal to x = 0.34 m. Also, this distance increases with the duration of the injection operation. Because saturation of water in the core is constant (11%), saturation of oil will not be higher than 89%. Figures 12a-12d show that the gas in the simulation results are more than experimental data.

Conclusion
Diffusion is the most important mechanism of mass transfer during N 2 injection between the matrix and the fracture. In the early simulation times, N 2 and C 1 are predominantly transported into the core by molecular diffusion of gas while C 5 is transferred into the core principally by oil convection. Moreover, during CO 2 injection test, diffusion and convection are both considerable. At the onset of the simulation, the diffusion mechanism is the dominant recovery mechanism, but over time, the role of the convection mechanism in C 1 and C 5 recovery increases. In addition, light component recovery is more than heavy component recovery. After 28 days of N 2 injection test, the lighter component recovery is about 50%, and the heavier component is recovered about 40%. Over time, the percentage of components recovery increases.
During N 2 injection experiment, at the primary simulation times, N 2 and C 1 are mainly transported into the core by molecular gas diffusion, and the role of the convection mechanism increases as the simulation is proceeded. C 5 is transferred into the core mostly by oil convection. Thus, for each component the mechanism of transport depends on the boundary conditions of matrix-fracture (state A or state B). For state A, boundary conditions of the fracture-matrix, gas molecular diffusion remains the most important mechanism in the transfer of N 2 and C 1 into the core at the simulation end times. Moreover, for state B, gas convection is the most important mechanism for the transfer of N 2 and C 1 into the core at the end of the simulation. Also, oil transfer, regardless of the boundary conditions between the matrix and the fracture (state A or state B), is the most important mechanism of C 5 transport into the core for the entire simulation time.
Ultimately, it is found that during the initial simulation times, when CO 2 is injected, CO 2 , C 1 , and C 5 are often transported into the core by the diffusion of oil and gas. In addition, the convection of oil and gas in the transfer of CO 2 , C 1 , and C 5 into the core increases with time, but the role of the diffusion mechanism during the simulation cannot be ignored.
As the time of nitrogen injection into the core increases, the maximum percentage of gas saturation in the core decreases and more gas escapes from the pores of the core. The convection mechanism has a great effect on preventing the pressure drop of the reservoir during injection operations.

Supplementary materials
The