Pore-Scale Modelling of Solvent-Based Recovery Processes

: A need for a reduction in energy intensity and greenhouse gas emissions of bitumen and heavy oil recovery processes has led to the invention of several methods where mass-transfer-based recovery processes in terms of cold or heated solvent injection are used to reduce bitumen viscosity rather than steam injection. Despite the extensive numerical and experimental investigations, the ﬁeld results are not always aligned to what is predicted unless several history matches are done. These discrepancies can be explained by investigating the mechanisms involved in mass transfer and corresponding viscosity reduction at the pore level. A two-phase multicomponent pore-scale simulator is developed to be used for realistic porous media simulation. The simulator developed predicts the chamber front velocity and chamber propagation in agreement with 2D experimental data in the literature. The simulator is speciﬁcally used for vapor extraction (VAPEX) modelling in a 2D porous medium. It was found that the solvent cannot reach its equilibrium value everywhere in the oleic phase conﬁrming the non-equilibrium phase behavior in VAPEX. The equilibrium assumption is found to be invalid for VAPEX processes even at a small scale. The model developed can be used for further investigation of mass transfer-based processes in porous media.


Introduction
Vapor extraction (VAPEX) is the solvent variant of the steam-assisted gravity drainage (SAGD) process (steam is replaced by solvent) which is driven by the solvent diffusion/dispersion mass transfer. In this process, the solvent is injected into a horizontal well at a pressure close to its dew point, and bitumen viscosity is reduced by solvent diffusion/dispersion. The bitumen with lower viscosity drains toward the producer which is located below the injector. The first version of this process was co-injection of solvent with hot water in the horizontal well located at the top of the reservoir [1]. The main mechanism was to distribute the injected heat in regions far from the injector. The VAPEX process is analogous to SAGD but requires much less energy and is aimed at reducing energy input and greenhouse gas emissions, assuming it is as effective for oil production as steam. In addition, it reduces or eliminates the water treatment needed for the SAGD process; this would result in a significant reduction in capital investment and operational cost. Mokrys and Butler [2] stated that 0.2-0.5 kg of propane was needed per kg of the produced oil, whereas 3 kg of steam is needed for 1 kg of oil. Energy consumption in the VAPEX process is 3% of the energy required for SAGD for the same production rate in terms of latent heat of vaporization of water and the solvents [3]. Asphaltene precipitation can be beneficial in the VAPEX process through in situ upgrading of bitumen [4]. On the other hand, asphaltene precipitation could result in low recovery of bitumen due to pore plugging. The asphaltenes can flow downwards in the chamber and damage the formation [5]. In addition, VAPEX is beneficial in reservoirs with a gas cap, bottom water, low thermal conductivity, high water saturation, clay swelling, and formation damage [6]. Steam condensation in SAGD may damage the reservoir in sensitive clay formations, which would not occur in the dry VAPEX process [2].
Oil viscosity, interfacial tension in the vapor chamber, reservoir heterogeneity, chamber dimensions, and diffusivities are the parameters that determine the drainage rate in the VAPEX process [3]. Das (1998) [7] stated that the addition of solvent could reduce vaporliquid surface tension, which is the key to a reduction in the residual oil saturation.
Considering the benefits of VAPEX such as low energy intensity, reduced water treatment, and in situ upgrading and its application for thin reservoirs makes this process a good alternative to steam-based processes. However, the bitumen production rate is very low compared to steam-based processes such as SAGD. The previous investigations showed that the oil rate of the VAPEX was only 1/20 of the oil rate of SAGD. The oil mobilization in VAPEX is controlled by diffusion mass transfer, which is much slower compared to the heat transfer in SAGD [8]. The injection rate of the solvent in VAPEX can play an important role in the recovery of heavy oil. If the injection rate is higher than the flow rate by gravity drainage, the solvent will fill some part of the already swept region instead of reaching the chamber edge. On the other hand, if the injection rate is lower than the gravity drainage rate, there will be constant liquid hold up [9].
As noted above, VAPEX is controlled by molecular diffusion, which is a much slower process compared to thermal diffusion and therefore it is expected that the oil rate would be much lower than that for SAGD. However, the extraction rate in the porous media laboratory tests was about 10 times higher compared to the results predicted by the tests in Hele-Shaw models. This difference is due to the higher interfacial area of contact in porous media, transient mass transfer at the interface, and film drainage [10]. Singhal et al. (1996) [3] suggested using longer horizontal wells to increase reservoir contact and the oil production rate.
A variety of solvent-based processes and injection pattern have been proposed in the literature. Butler and Mokrys (1998) [11] suggested solvent injection into the aquifer noting that it can benefit the process in two different ways. First, the solvent will not dissolve in water, and there is no heat loss. As a result, the aquifer is used for initial injectivity. Second, the water layer will underride the mobilized oil toward a production well. Zeng et al. (2008) [12] suggested the VAPEX process in a Tee-well configuration. In this model, a second injector is drilled horizontally above the injector and producer, and perpendicular to them. They stated that this model could improve the performance of VAPEX through two different mechanisms: first, by having a higher injection pressure in the injector above, the gravity drainage force would be increased, and second, a new diluent profile will form perpendicular to the former diluent profile. Singhal et al. (1996) [3] note that in a lower-viscosity reservoir such as those in the Lloydminster area, a series of vertical injectors together with one horizontal producer can increase the size of the vapor chamber. VAPEX/Wormhole process was suggested by Metwally (1996) [13], wherein the wormholes induced by the sand production could be used to distribute the solvent through the reservoir. Similar to fractures, widespread wormholes would be harmful, as the solvent would channel to undesirable regions. Metwally (1996) [13] suggested that the injection of natural gas with vaporized propane in cyclic form. Butler and Mokrys (1991) [1] proposed wet VAPEX. In this process, hot water and solvent are injected from the top of the reservoir, and the producer is located at the bottom of the reservoir.
Cyclic production with continuous solvent injection (CPCSI) was proposed by Jiang et al. (2013) [14]. They suggested cyclic production from the producer with continuous solvent injection. Das (2008) [15] investigated the impact of non-condensable gas (NCG) addition to propane in VAPEX numerically. He showed that methane would accumulate near the top and expanding edge of the chamber, which could help to confine solvent in thin reservoirs. He stated that injection of hot solvent would reduce the solubility of the solvent in bitumen. It should be noted that none of the VAPEX versions has been commercially successful anywhere close to the experimental results.
There are very limited experimental and modelling studies of VAPEX process at the pore-scale level. James and Chatzis (2004) [16] conducted a pore-scale experiment for the VAPEX process using butane as a solvent. Kim (2017) [17] conducted a series of experiments using a micromodel to study the SAGD and SA-SAGD processes in porous media. She considered propane, butane, hexane, and condensate as solvents for the solvent-assisted steam-assisted gravity drainage (SA-SAGD) process. Experimental results indicated that condensate and hexane have the highest recovery factor while the ultimate recovery factor for the propane-SAGD process is even lower than for the SAGD process. The butane-SAGD process had a lower recovery factor than SAGD initially, but the ultimate recovery factor was higher than for pure SAGD. Microscopic analysis of experimental data provided insight into interface region and fingering, wettability alteration and phase (saturated, aromatic, resin, and asphaltene) partitioning due to the presence of a solvent. Mohammadzadeh et al. (2010) [18] conducted a series of experiments using glass micromodels to investigate the SA-SAGD process at the micro-scale level. They used pentane and hexane as a solvent for their experiment and co-injected with steam. They considered four different micro-models for their experiments. Their experiments showed two simultaneous drainage processes (capillary drainage displacement and film-flow drainage) for oil recovery in the SA-SAGD process. Both molecular diffusion and mass convection occur in porous media. However, due to the convective nature of the process in the mobilized region; convection is the dominant mass transfer mechanism. In addition, they observed asphaltene precipitation in the model which results in permeability impairment. Pore-scale experimental studies are summarized in Table 1. A pore-scale simulator is developed in this study for modelling interfacial mass transfer in porous media. The interfacial mass transfer phenomenon is a scale-dependent process in which the diffusion path and interfacial area play an important role. Solvent and gas extraction, sequestration of carbon dioxide into the subsurface, and wastewater treatment are some of the examples for the application of interface mass transfer. This mass transfer can change the composition of the fluids or it can affect fluid properties, e.g., viscosity and density. Alteration of bitumen viscosity through its dilution with solvent is the main mechanism in solvent-based recovery processes. To model this process at the porous media, a general volume of fluid (VOF) model is developed for interfacial mass transfer in Open Source Field Operation and Manipulation (OpenFoam) toolbox. This tool has been used for average properties calculation and describing the interface is beyond the scope of this work. In the development process, where needed, appropriate validation using analytical models and experimental data are provided.
VAPEX using butane as the solvent is modelled using the pore-scale simulator developed in this study. It is assumed that oil viscosity is 23,000 cp initially, which provides enough mobility for solvent movement in the model. Interface movement of vapor phase and chamber development in porous media are investigated at different times for VAPEX. Capillary effect and oil trapping phenomena are captured in the simulator. Nonequilibrium constant (contrary to equilibrium constant (K-value)) is calculated for cold VAPEX, confirming the need for non-equilibrium phase behavior in the reservoir.

Numerical Model Description
In the present study, the fluid dynamics is solved using phase tracking method of VOF. In this method, the position and form of an interface is tracked using a phase indicator field (α p ) which is volume fraction defined as below: where V p is the volume of phase p and V is the total volume. α p = 1 means the grid is occupied completely with phase p while α p = 0 means there is no phase p in the grid. Any value between these two boundary attributes to the phase interface in the VOF. It should be noted for a two-phase system, we always have In this method, the global variables are averaged over a cell as below: where p and v are pressure and velocity in a grid. α 1 is the phase indicator of phase 1 (e.g., vapor phase). Subscripts refer to two phases that occupy a gridblock (e.g., vapor and liquid). Other variables such as density and viscosity are calculated with the same method. Mass conservation for two phases would be In this equation, . m is the mass flux between two phases which is calculated from the mass balance equation. ρ is the density. The overall momentum equation is where µ is viscosity and F c is the interfacial force. In this study, the Continuum Surface Force (CSF) approximation method has been used for interfacial force calculation [26].
where σ is surface tension. This method is used by several authors to model interfacial mass transfer. Haroun [28] developed a model based on conditional volume-averaging technique and provided a detailed procedure for arithmetic and harmonic averaged diffusion coefficient. They also addressed artificial mass transfer using the Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) interface capturing scheme [29].
The first model for mass transfer was developed by Haroun et al. [27]. Later, their model was improved by Marschall where C, D, and J are concentration, diffusion coefficient and flux, respectively. j represents the phase of interest. They expressed the diffusion term as Assuming the continuity of fluxes and using Henry's law for describing equilibrium at the interface, where H is Henry constant. The diffusion rate can be expressed as As can be seen in the last term in the RHS of the above equation is non-zero only at the interface (∇α = 0). Using jump condition at the boundary by Henry's law, the concentration difference at the boundary is Using this expression, the final diffusion-advection equation will be as follows: The mass transfer flux at the interface is calculated based on concentration equation using following equation: where n 12 is the vector normal to the interface and w is speed of interface. Haroun et al.
(2010) [27] stated that their model can be used with any diffusion coefficient averaging scheme while using harmonic averaging can improve their prediction. However, Deising et al. (2016) [28] showed their model is only valid when harmonic averaging is used. A pure solvent-rich vapor phase is considered in VAPEX processes in which solvent only diffuses in the oleic phase. An example of this process assuming the volume change due to dissolution is negligible is shown in Figure 1. The model is a 2D capillary with dimensions of 1 mm × 0.4 mm with known concentration at the left and right wall of the tube. Solvent concentration in the oleic phase at different times is shown in Figure 2. As this figure shows, the concentration of solvent in the vapor phase is constant while the concentration of solvent in the oleic phase is increasing until it reaches the equilibrium value. It should be noted in this model, it is assumed the density and viscosity of the two phases are equal. A validation with analytical model steady-state and transient cases is provided in Appendix A.

Physical Model Description
The VAPEX performance in the porous media is investigated assuming pure gaseous solvent and compared with James's (2009) [31] 2D etched glass micromodel experiment. Reported parameters (Table 2) are used to build the porous media model in AutoCAD software. In this porous medium, a permeable channel was constructed at the left side of the model to facilitate oil and solvent injection. The butane vapor was injected through a sidechannel and the mobilized bitumen produced through the production outlet at the bottom of the model. Therefore, there is no forced convection mass transfer in the system, but the flow is mainly due to gravity similar to VAPEX process. A small model with 15 × 15 grains is considered to optimize simulation execution time. The snappyHexMesh command in OpenFOAM has been used in three steps: castellated mesh (removal of grains grids in cartesian coordinates), snapping the surface (removal of the jagged surface), and Mesh Layer (addition of grid layers in the vicinity of boundary surface). Using snappyHexMesh, the model was meshed, resulting in 102,900 grid blocks, where 5628 cells are prismatic, and the rest are hexahedral (Figure 3).  The boundary conditions were selected based on operational and mechanistic conditions and also those in the experimental study. A solvent channel is considered at the left of the model, providing the solvent influx into the model. It is assumed that the top and right sides of the model have no-flow boundary conditions and the oil drained only through the lower part of the model. It is assumed that the solvent is not injected into the model, rather it moves in the model by gravity drainage. The porous medium simulator developed in this study is used to model solvent propagation and bitumen drainage in this model. No-slip boundary condition is assumed for grains, meaning the velocity is zero adjacent to the solid surface. A least squares cell-point-cell stencil has been used in the present study for gradient scheme as depicted in Figure 4. A mesh sensitivity analysis has been done and the results are shown in Table 3 confirming sufficient resolution for this study.

Fluid Properties
The density of the oleic phase is taken to be constant and its viscosity is considered to be a function of composition via the following equation: where x and µ are mole fraction and viscosity, respectively; bs, b, and s subscripts refer to the oleic phase (solvent + bitumen), bitumen, and solvent. Oil density and viscosity are set to 23,000 cp and 987.0 kg/m 3 , respectively, in the simulation model. Bowman (1967) [32] presented the surface tension of bitumen as a function of temperature in the presence of air. Baptista (1983) presented these data in a tabular form which indicated a value of 34.7 mN/m at a temperature of 25 • C [33]. Das and Butler's (1996) diffusion model is used for butane diffusivity in oil leading to a value of 9.76 × 10 −11 m 2 /s. A K-value of 2 is considered for butane in Athabasca bitumen at a temperature of 25 • C and atmospheric pressure. In this simulation, the porous model is open to oil drainage from the bottom side. The rock is oil-wet with a contact angle of zero and vapor is the non-wetting phase. The contact angle is assumed to be static and ascending and descending wetting effects are ignored.

Results
The model developed is used for investigating VAPEX process in a 2D porous media as described in the methodology section. The main purpose of the simulation is to track solvent path in the porous media and its dissolution in the oleic phase. The efficiency of solvent injection is investigated in terms of two phenomena: first, solvent interfacial transfer from vapor to oleic phase, second, its diffusion in the oleic phase and reduction in viscosity. The VAPEX interface advancement can be seen in Figure 5, initially and after 2000 s, 4000 s, and 6000 s, respectively. Butane started diffusing into the oleic phase from the interface of oleic and gaseous phases. The simulated solvent fraction distribution of butane in the porous medium is shown in Figure 6 after 1.6 h. As the figure shows, the solvent started diffusing into the oleic (oil-rich) phase from the two-phase interface. Considering the K-value of 2, the equilibrium mole fraction of butane in the oleic phase at the equilibrium condition is 0.5. As Figure 6 shows at the interface and in the trapped oil within the swept regions, the mole fraction of butane reached the equilibrium value. However, after 1.3 h, even for a model with the size of 3 cm × 3 cm, the butane mole fraction in the oleic phase did not reach the equilibrium value and is close to zero away from the interface. Solvent diffusing into the oleic phase reduced the viscosity at the interface of vapor and oleic phases. Viscosity distribution in the pore-scale simulator at 1.6 h can be seen in Figure 7. Based on Figures 6 and 7, it can be observed that the solvent can reach the equilibrium value only in the vicinity of the interface. Figure 6 indicates not only that solvent cannot dissolve in the oleic phase, but its distribution is far from homogeneous, which is the main assumption in numerical methods used for solvent recovery process modelling. Based on the mole fraction and viscosity distribution, three regions can be observed in the VAPEX process; vapor chamber, chamber edge (mobilized oil), and immobile oil. The vapor chamber comprises pure butane which propagated upward and sideways and includes residual (trapped) oil. At the chamber edge, the mobilized oil drains downward due to gravity. Oil was produced mainly from this region and as a result, the rate of chamber advancement is key to the performance of the VAPEX process. As Figure 7 shows that the mobilized oil region is narrow and solvent diffusion occurs over a very thin interval.
The capillary effect can be observed in the pore-scale simulation in the curved interface at the throat junction with pore body. Figure 8 shows the effect of capillarity and wettability on the film of gaseous phase interface position. Propagation of the solvent phase in the model is in agreement with the experimental results. The oil drainage at the chamber edge can also be observed using the simulator developed in the present study ( Figure 9). The figure shows the vapor phase interface in the second row from the top of the porous media and sixth to eight columns from the left. As the figure indicates, the non-wetting phase (butane vapor) drained the wetting phase of mobilized oil. The invasion will start by filling the first pore by butane vapor and continue by invading more pores. As seen, the propagation of the vapor phase is not uniform in the porous medium as a result of capillary force and there is residual oil trapped inside the butane invaded zone.  Figure 11 shows the residual oil saturation in the second row from the top and seventh column from the left. As the figure shows, the low viscous oil starts flowing in the throat, and the residual oil saturation is decreasing.   Figure 12 shows the butane mole fraction distribution after 4.2 h. The mole fraction of butane in the oleic phase is equal to its equilibrium value at the interface and its value is decreasing with the distance from the interface. Butane mole fraction in the oleic phase has not reached its equilibrium mole fraction after 4.2 h. In this case, viscosity distribution in the oleic phase is solely dominated by the diffusion coefficient. In conventional VAPEX, where it is assumed a cold solvent is injected into the reservoir, the rate of pore invasion of the solvent is controlled by the rate of diffusion where mobilized bitumen can flow to the well and let the solvent proceed in the porous medium. Using the model developed in the present study, non-equilibrium constants are calculated after 4.2 h by averaging mole fractions in gaseous and oleic phases. It was found that after this time, the mole fraction ratio is 4.65 which is much higher than the equilibrium value of 2 ( Figure 13). Using the results of the pore-scale model, the averaged viscosity in the oleic phase is also calculated ( Figure 14). As the figure shows the solvent was not able to reduce oil viscosity significantly after 4.2 h for cold solvent. Time-dependent variables of non-equilibrium constant and viscosity can be implemented in upscaled models through appropriate upscaling approach, e.g., volume averaging.  The results in this section show low interface advancement velocity in conventional VAPEX. In this process, chamber growth mainly dominated by diffusion mass transfer, which leads to low oil rates.
where C is concentration. Applying the boundary condition at (−l) for C Left and +l for C Right , the equations will transfer to The extra equations for the analytical solution come from the continuity of concentration and its derivative at the interface.
where H is Henry constant and D L and D G are diffusion coefficients in gas and liquid phase. Solving these equations together, the resulting solution is For the comparison, a 2D capillary is considered with dimensions of 1 mm × 0.4 mm with known concentration at the left and right wall of the tube is considered ( Figure A1). The grid blocks are assumed to be 1000 × 40 with a total number of 40,000. The left and right boundary conditions are assumed to be 1 and 0; respectively, with a diffusion coefficient of 2 m 2 /s and 1 m 2 /s. The simulation continued until it reached a steady-state condition. The numerical results are compared with the analytical solution in Figure A2. As the figure shows, the numerical results for the steady-state case are in good agreement with analytical results. For transient case, the tube is considered to be infinite which is only valid when the numerical results are compared at early times when the concentration gradient is only located near the interface which makes this assumption valid. This equation needs initial conditions as well as boundary conditions. The initial conditions are as follow: C(x < 0, t = 0) = C Next, a 1D model with length 2e is considered. It has been assumed that D L = 0.1D G and the Henry constant is 25. The model discretized by 500 grid blocks. Similar to the previous test case, the concentration distribution is modelled for early times and compared with analytical results [27]. The local mass transfer coefficient for each phase (K p,local ) is calculated as follows: In this equation, p is the phase under study, and n p is normal to the interface pointing to the phase p. ∆C is concentration difference at the interface and a reference concentration far from the interface. Φ is interface flux for concentration jump condition. Considering the analytical model for concentration distribution in the 1D model, the local mass transfer coefficient (K p ) and Sherwood number (Sh p ) at early times are The analytical model compared with numerical data at two dimensionless number of 0.005 and 0.01, where dimensionless time is defined as t * = tD r e 2 (A15) In this equation, D r is the diffusion coefficient in the reference phase, which is considered to be liquid in this study.