Mechanistic modelling of non-equilibrium interphase mass transfer during solvent-aided thermal recovery processes of bitumen and heavy oil

Mechanistic modelling of non-equilibrium interphase mass transfer during solvent-aided thermal recovery processes of bitumen and heavy oil

Fuel 241 (2019) 813–825 Contents lists available at ScienceDirect Fuel journal homepage: www.elsevier.com/locate/fuel Full Length Article Mechanis...

5MB Sizes 0 Downloads 20 Views

Fuel 241 (2019) 813–825

Contents lists available at ScienceDirect

Fuel journal homepage: www.elsevier.com/locate/fuel

Full Length Article

Mechanistic modelling of non-equilibrium interphase mass transfer during solvent-aided thermal recovery processes of bitumen and heavy oil

T

Abdullah Al-Gawfia, Hossein Nouroziehb, Ehsan Ranjbarb, Hassan Hassanzadeha, , Jalal Abedia ⁎

a b

Dept. of Chemical & Petroleum Engineering, Schulich School of Engineering, University of Calgary, Calgary, Alberta, Canada Computer Modeling Group Ltd, Calgary, Alberta, Canada

GRAPHICAL ABSTRACT

ARTICLE INFO

ABSTRACT

Keywords: Thermal recovery Bitumen Solvent-aided recovery Non-equilibrium Interphase mass transfer

One of the important mechanisms involved in solvent-aided thermal recovery processes is interphase-mass transfer phenomenon. However, in current reservoir simulation models a local equilibrium is commonly assumed such that a simulation grid block is at instantaneous equilibrium. The assumption of local equilibrium often fails at larger scales or in situations where flow velocities are considerable compared to that of mass or heat transfer which results in significant discrepancies between reservoir simulation predictions and field scale observations. In this work, solvent-aided gravity drainage of bitumen was simulated with propane as a solvent. The effect of non-equilibrium mass transfer was included in the model to simulate the process using a kinetic approach. The results show that the assumption of the local instantaneous equilibrium results in lower oil recovery for the typical field scale simulation models. The analysis shows that this mismatch can be mitigated fairly through the inclusion of the non-equilibrium interphase mass transfer. Correlations for the non-equilibrium interphase mass transfer coefficients for propane/bitumen mixture were developed which can be used as guidelines for field scale modeling of solvent-aided thermal recovery processes.

1. Introduction Canada ranks in the top ten producing countries in the world and the largest exporter of oil to the United States. However, the challenge we face in Canada and around the world is not simply producing



enough energy to meet the world increasing demand, but to reduce the environmental impact and greenhouse gas (GHG) emissions when developing energy resources. As Canada has the third largest reserves of crude oil, nearly 170 billion barrels of recoverable oil, 165 billion of which are located in the oil sands [1,2], Canada’s energy future depends

Corresponding author. E-mail address: [email protected] (H. Hassanzadeh).

https://doi.org/10.1016/j.fuel.2018.12.018 Received 14 September 2018; Received in revised form 4 December 2018; Accepted 5 December 2018 0016-2361/ © 2018 Elsevier Ltd. All rights reserved.

Fuel 241 (2019) 813–825

A. Al-Gawfi et al.

on the oil sands industry. Currently around 55% of the oil sands is being extracted with in-situ technologies. However, nearly 80% of the oil sands reserves is only recoverable using in-situ thermal recovery processes [3,4]. Two thermal recovery methods stand out to be the most viable and commercially practical for exploiting extra-heavy and highly viscous oil reservoirs: Steam-Assisted Gravity Drainage (SAGD), and Cyclic Steam Stimulation (CSS) processes. Despite the commercial success of these thermal recovery processes, solvent-aided thermal recovery processes recently gained increased industrial interest for their potential to achieve higher energy efficiency, reduced environmental impact, and increased economic viability. Alberta oil sands producers aim is to develop new technologies that will improve the performance of in-situ thermal recovery projects and make them more competitive with other oil plays in terms of operational costs, environmental footprint, and GHG emotions. Alberta oil sands operators have been experimenting on new and innovative technologies to reduce energy intensity of in-situ operations by reducing natural gas use and water consumption. These innovations include co-injection of solvents: Solvent Assisted-SAGD (SA-SAGD), using alternative sources of thermal energy such as electromagnetic energy: Enhanced Solvent Extraction Incorporating Electromagnetic Heating (ESEIEH), and pure solvents injection (Nsolv). Several field trials of solvent-based recovery processes have been carried out and field results were either mixed or inconclusive due to the lack of knowledge of the underlying physics involved in the recovery processes. SAGD with solvent co-injection was proposed and developed by Nasr et al. [5]. The fundamental mechanism of the process is that solvent condenses with steam around cold formation interface of the steam chamber causing oil dilution and viscosity reduction. SAGD with solvent co-injection processes needs to be operated in an optimized thermodynamic window (pressure, temperature and solvent composition) to be effective and provides the maximum productivity. SAGD with solvent co-injection technologies are considered among the promising technologies and is currently being piloted by several operators in Alberta. It is expected that energy intensity reduction would be in the range of 10–30% on a per barrel basis and it could reduce GHG emissions by 15–35% [6]. In the case of pure solvent injection technology (such as Nsolv), there is an expected 75% reduction in energy intensity and 75–80% reduction in direct fuel-derived emissions. In the case of solvent-assisted processes such as Solvent-Aided Process (SAP) and Expanding Solvent-SAGD (ES-SAGD), there is a forecasted reduction in steam-oil-ratio (SOR) of 33–36%, 35% reduction in natural gas consumption, and 15–20% emissions reduction relative to the SAGD base case [7]. Despite being one of the most promising technologies, solventaided processes are quite complex due to the inherent complexities of the associated mechanisms of heat and mass transfer coupled with complex phase behaviour. As steam delivers its latent heat to the reservoir and condenses at the boundary of the steam chamber, both bitumen viscosity and density are decreased by the combined effect of temperature increase and solvent dissolution into the oil [8]. In the latter mechanism, mass transfer parameters become very important to understand the process. Diffusion coefficient and interphase mass transfer coefficient are the basic mass transfer parameters in understanding the dissolution of vapor solvents into heavy oils. One of the most basic assumptions when modeling multiphase flow in porous media is the local thermodynamic equilibrium within an averaging volume. It is assumed that concentration of species in contact are reached instantaneous equilibrium everywhere in the averaged volume. Such an assumption may be acceptable in the cases of fast mass transfer processes. However, local thermodynamic equilibrium assumption will give erroneous results in the cases of slow mass transfer processes where the characteristic time of mass transfer is small compared to flow velocities of the system. Also, since interphase mass transfer process is inherently a pore-scale process as it takes place across fluid–fluid interfaces [9–11], local thermodynamic equilibrium

assumption may be valid when modeling the process at small scale grid sizes. However, local thermodynamic equilibrium assumption may yield erroneous predictions when modelling solvent-aided recovery processes at the macro-scale level (reservoir simulation grid block scale). Given its importance, the concept of non-equilibrium interphase mass transfer is not new to literature. Several researchers investigated the non-equilibrium phenomena and tried to capture the physics of the non-equilibrium mass transfer when modelling multiphase flow in porous media. Mayer et al. [12] expressed the non-equilibrium mass transfer as a first-order kinetic as: Q K = k K a (C K, s C K ) , where Q K [M / L3t ] is the rate of interphase mass transfer of component K from phase to phase , k K [L / t ] is the interphase mass transfer coefficient, [1/ L] is the specific interfacial area separating phases and , C K, s [M / L3] is the solubility limit of component K in phase , and C K [M / L3] is the concentration of component K in phase . Hunt et al. [13] applied this model to estimate the non-aqueous phase liquids (NAPLs) dissolution rates. Due to the complexities exerted by the heterogeneity of the porous media, it is difficult to estimate the specific interfacial area parameter . Therefore, several authors lump interphase mass transfer coefficient k K and the specific interfacial area parameters into a single parameter k, effective interphase mass transfer coefficient, [14–18], which yields to a simplified form: Q = k (Cs C ) where Cs is the solubility limit, and C is the actual concentration and k is the effective interphase mass transfer coefficient. Nghiem and Sammon [19] modeled the non-equilibrium mass transfer using a compositional simulator in which the oil and gas phases approach thermodynamic equilibrium through a rate process. The results of their simulations of a gas injection displacement process in a core showed that non-equilibrium cases lie between equilibrium cases and nomixing cases. Indrupsky and Lobanova [20] and Lu et al. [21] modeled non-equilibrium mass transfer using component mass transfer rate, which is proportional to the difference of the component chemical potentials in the phases, i.e.: i, L V = (µi, L µi, V ), where, µi, L and µi, V are the chemical potentials of the component i in the liquid and vapor phases, respectively, and is a function of reservoir characteristics and dynamic flow properties. Through their analysis, they showed that the necessity for non-equilibrium phase behaviour depends on simulation scale. Wu et al. [22] assessed the effect of the non-equilibrium mass transfer on the productivity of a single well producing from a gas condensate field. Their model incorporated the non-equilibrium mass transfer term on a compositional simulator and they used the model of Wilkins et al. [23] described as ja0 = k 0 (Cs C ) , where Cs is the solubility limit, and C is the actual concentration and k 0 is the effective interphase mass transfer coefficient. Their simulation results reveal that the non-equilibrium mass transfer leads to a reduction in the condensate saturation in the near wellbore region and therefore slower reduction in well productivity. Ivory et al. [24] conducted reservoir simulation and experimental studies of solvent cyclic injection (CSI) and concluded that the oil production rate depends on rate of solvent dissolution and exsolution from the process. Chang and Ivory [25] developed numerical simulation models of the CSI process which incorporate non-equilibrium rate equations representing the delay in solvent reaching equilibrium concentration as it dissolves or exsolves in the oil in response to changes in the pressure and/or gas-phase composition. Their formulations are described by: Ni, L/ t = ki Non1 (x i, g x i, L) n1Ngn2 yi, g , where ki is rate constant for solvent dissolution, Ng is the moles of the gas phase, Ni, L is the moles of solvent in the oleic phase, No is the moles of the oil phase, n1 and n2 are constants. They concluded that the assumption of instant equilibrium results in a 23% reduction in oil production compared to the nonequilibrium cases. Soh et al. [26] conducted experimental and numerical studies of pure methane, methane/propane mixture, and CO2 as CSI in heavy oil. In their numerical simulation study, they modeled the nonequilibrium behaviour using two reaction terms. The first reaction illustrates the clustering of dissolved gas bubbles to be trapped in pore 814

Fuel 241 (2019) 813–825

A. Al-Gawfi et al.

space, and the second reaction explains the development of free gas phase. Both reactions are defined in the numerical simulation using Arrhenius equation: r = Ae E / RT C n, where r is the reaction rate, A is a preexponential factor, E is activation energy, R is the universal gas constant, T is temperature, and C is the component concentration. They concluded that methane-propane mixture reduces the foaminess effect of pure methane solvent which turns the process close to instantaneous equilibrium. Jia et al. [27] conducted numerical simulation to model solvent vapor extraction process (SVX) incorporating non-equilibrium solvent solubility phenomenon. They used a time-dependent mass transfer term of the solvent vapor dissolving into oil phase using a chemical reaction term. Their study showed that simulations conducted using non-equilibrium method resulted in more realistic relative permeability curves that closely resembles curves obtained from classical core floods. In addition, they concluded that numerical dispersion effect for larger grid block sizes can be compensated by tuning the reaction rate frequency factor. Nourozieh et al. [28] conducted numerical simulations to model noncondensable gas injection for a hybrid SAGD process. They incorporated rate-dependent dissolution and ex-solution of solution gas in their simulation using chemical reactions. They noted that by incorporating rate-dependent dissolution and ex-solution terms for the co-injection of non-condensable gases with steam, a slight improvement in oil production with reduction in steam-oil-ratio was observed. Also, they concluded that real field behaviour was captured, and steam chamber conformance was not affected by initial solution gas when using the non-equilibrium dissolution and ex-solution reaction terms. Saha and Peng [29] proposed a correlation to model non-equilibrium mass transfer through the use of a partition coefficient to modify the equilibrium k-values as given by: K i = Ki (1 Ei ), where Ei is a measure of the degree of non-equilibrium of the system with respect to component i. The results of their work show that the equilibrium model shows slightly more optimistic results compared to both experimental and non-equilibrium cases. Rangriz-Shokri and Babadagli [30] conducted lab experiments to assess the feasibility of heavy oil recovery by CO2 injection while considering non-equilibrium foamy oil behaviour. They further numerically simulated depletion lab experiments to determine the parameters of reaction kinetics using Arrhenius equation: r = Ae Ea /RT C n . They developed a power-law scaling relation for the reaction kinetics parameters for field scale grid size modelling of the process. They further used the scaling rule to upscale the kinetic parameters to history match Cold Heavy Oil Production with Sand (CHOPS) field case with 15 wells in Alberta. Niessner and Hassanizadeh [9–11] presented an interfacial-area-based model in which they included interfaces explicitly in their formulations to describe kinetic interphase mass and energy transfer in a thermodynamically-based approach. In their formulations, they included a micro-scale diffusive flux through which kinetic interphase mass transfer is captured. The aim of this work is to get a better understanding and proper modelling of solvent-assisted recovery processes at the macro-scale level (reservoir simulation grid block scale) while capturing the physics of the interphase mass transfer mechanism by avoiding the local thermodynamic equilibrium assumption in the simulation models. The mechanism of solvent dissolution into bitumen due to solvent diffusion are first addressed and the key parameter of diffusive dominant interphase-mass transfer coefficient for several solvent/bitumen binary mixtures is determined. Then the non-equilibrium of mass transfer is modelled through the use of rate dependent reaction terms substituting the typical equilibrium k-values method currently used in commercial simulators.

Fig. 1. Mass transfer by diffusion – model description.

system. When the dominant mass transfer mechanism is pure diffusion an analytical treatment is possible subject to few simplifying assumptions. It is assumed that the liquid phase in Fig. 1 is a non-volatile component (bitumen) and therefore only gas is diffusing into the liquid phase. Also, swelling of the bitumen is assumed negligible. Sheikha et al. [31] analysis of pressure decay tests show that ignoring swelling of bitumen when estimating diffusion coefficients for CH4, CO2 and N2 resulted in a range of error between 5.9% and 9.3%. Finally, it is assumed that the diffusion coefficient is constant, and the system is at isothermal condition. In the numerical simulations (see Section 4) the gas phase is contacted with the immobile bitumen from the bottom of the domain based on the application of our interest. However, in the analysis of pure diffusion the gas phase is placed in top of bitumen to avoid the influence of convective mixing. Through employing Fick’s law, and the initial and boundary conditions presented in Fig. 1, the diffusive interphase mass transfer coefficient can be expressed by:

k=

D H (1

CD (zD = 0) zD

CD )

(1)

where k is the diffusive interphase mass transfer coefficient [L/t], D is the diffusion coefficient [L2/t], H is the thickness of the liquid column [L], and finally CD is the dimensionless concentration of the diffusing component which can be found analytically using separation of variables method with the initial and boundary conditions presented in Fig. 1, [32]:

CD (zD , tD) = 1

4 n=0

2n + 1 2

1 sin 2n + 1 2

tD

2n + 1 2

zD exp

(2)

The analytical solution for the diffusive interphase mass transfer process presented in Eq. (1) can be used for solvent/bitumen mixtures to estimate the interphase mass transfer coefficients at different pressure and temperature settings. The solvents studied are C1 to C5 and the bitumen samples used are JACOS, MacKay River and Surmont [33–35]. As noted from Eq. (1), all parameters are in dimensionless form except for H [m] and D [m2/s]. H is the column height of the liquid phase and is a property of the physical model which can be set to any arbitrary value for the analytical solution. To unify the solution of the analytical and numerical models that will be addressed herein, the value for H was

2. Pure diffusive interphase mass transfer coefficient When a gaseous solvent is brought into contact with liquid phase, as depicted in Fig. 1, it diffuses into the liquid phase as described by Fick’s law. Equilibrium concentrations are established instantaneously at the gas–liquid interface. As gas is further dissolving into the liquid, an equilibrium concentration, C , eq is eventually achieved throughout the 815

Fuel 241 (2019) 813–825

A. Al-Gawfi et al.

set to 1 m. D is the molecular diffusion coefficient of the diffusing solvent into liquid. Analytical determination of the diffusive interphase mass transfer coefficients requires prior knowledge of the corresponding molecular diffusion coefficients. However, experimental molecular diffusion coefficients are scarce for solvent/bitumen systems at different temperature and pressure conditions. Therefore, a well-known correlation developed for prediction of molecular diffusion coefficient for binary systems will be used to determine the molecular diffusion coefficients. The empirical correlation used in this study to estimate the molecular diffusion coefficient for liquids at infinite dilution is the one proposed by Hayduk and Minhas in 1982: DAB = 13.3 × 10 8T1.47µB10.2/ VA 0.791 / V A0.71, where DAB is the infinite dilution diffusion coefficient [cm2/s], MB [g/mol] is the molecular weight for bitumen in this case, T is for temperature in Kelvin, µB is the viscosity of bitumen in cP, and VA is the molar volume of the diffusing solvent at its normal boiling point [cm3/mol]. Several other correlations were proposed for estimation of molecular diffusion coefficients. However, it is worth noting that such empirical correlations are mainly accurate for the range of experimental data points, temperature and pressure settings used in their developments. For example, Wilke and Chang developed a correlation to estimate the molecular diffusion coefficient at infinite dilution. However, since their correlation is mainly designed for liquids, Hayduk and Minhas correlation was used. Furthermore, Hayduk and Minhas correlation was developed using solvents ranging from C5 to C32 diffusing into normal paraffins, which is similar to the system of our study, and only produces 3.4% error [36]. Fig. 2 presents the results for the pure diffusive interphase mass transfer coefficient for different solvent/bitumen binary mixtures using Eqs. (1) and (2). The dotted lines are for gaseous solvents and solid lines are for liquid solvents. As can be observed from the figure, lighter solvents have relatively higher mass transfer coefficient compared to the heavier solvents. This indicates that more deviations from equilibrium is expected for gaseous solvents. Hence, for modelling of solventaided recovery processes such as ES-SAGD and SAP, inclusion of the interphase mass transfer phenomena is relatively more important for gaseous solvents compared to heavier/liquid solvents. The results show that for a light solvent such as methane, mass transfer coefficient can be expressed as k = 2 × 10 8µB 0.593 , and for a heavier solvent such as pentane it can be expressed as k = 1 × 10 8µB 0.771 . Also, Fig. 2 shows that when heat is involved (low viscosity), interphase mass transfer is higher and thus more important to be properly modelled when designing solvent-aided thermal recovery processes regardless of the

solvent type used in the process. Nuske et al. [37] observed the same and mentioned in their paper that the assumption of local thermodynamic equilibrium is questionable when a heat source is present in the porous medium. The analysis described above forms the basis for the prediction of the mass transfer coefficient when more complex processes are involved as discussed in Section 4 of this paper. In the following section, fluid model required for the numerical simulation of solvent dissolution into bitumen is presented. In the numerical simulation that follows propane is used as solvent. 3. Thermophysical properties of propane/bitumen mixture The experimental data used to generate the fluid model for Surmont bitumen/propane mixture is obtained from the work of Nourozieh [33]. In the following, modelling of fluid density, viscosity and k-values is presented. 3.1. Density of bitumen/propane mixture To find density of bitumen/propane mixture as a function of pressure and temperature applicable in reservoir simulator CMG STARS [38], the following equation was used: i

=

o i exp

c pi (P

ct1i (T

Tref ) + ct 2i

Pref )

cpti (P

[(T + 273.15) 2

Pref )(T

(Tref + 273.15) 2] 2

Tref )

(3)

= Mwi / is the partial molar volume of component i, Mwi is where the molar mass of component i, io is the reference density obtained from equations given in Appendix (A1) and (A2) for both bitumen and propane, respectively, at reference pressure (101.325 kPa) and temperature (15.6 °C). It is worth noting that i (or io ) for propane is the pseudo partial molar volume. The mixture molar volume is calculated based on the assumption of no volume change upon mixing using the 2 following mixing rule: vo = i = 1 x i vi , where xi is the mole fraction of the components and vi is the component partial molar volume. The oil phase molar density in kg-mole/m3 can be obtained by 1/v0. Regressing experimental data of density for propane/bitumen system to find the constants: ct1i , ct2i , c pi , and cpti resulted in the values presented in the Table 1, along with the corresponding AARD%. vio

o i

Fig. 2. Estimated mass transfer coefficients for several solvent/bitumen mixtures, the dotted lines are for gaseous solvents and solid lines are for liquid solvents. 816

Fuel 241 (2019) 813–825

A. Al-Gawfi et al.

Table 1 Fluid density constants of Eq. (3) for propane/bitumen fluid model.

Ct1 (1/K) Ct2 (1/K2) Cp (1/kPa) Cpt (1/kPa K)

Table 3 Fitting parameters of viscosity correlation.

Bitumen density constants

Propane density constants

AARD%

0.0010524 −1.09 × 10−6 3.75 × 10−7 4.40 × 10−9

0.0097915 −1.58 × 10−5 1.0 × 10−6 2.02 × 10−7

0.285

It is worth noting that Eq. (3), A1, and A2 are valid for predicting bitumen densities and pseudo liquid densities of propane for temperature and pressure ranges of 50 T 190 C and up to 13.8 MPa, respectively.

Bitumen viscosity correlation parameters

c1

c2

Solvent viscosity correlation parameters

d1

d2

−16.6828

0.148144

24.10688

c3

AARD (%)

0.005782

13.94

d3

d4

d5

−0.00025

−2.09824

0.004773

−3.77505

(a)

3.2. Equilibrium K-values

Bitumen

The K-values for the simulation model constructed for this study are obtained through regressing the experimental K-value data for the saturated bitumen/propane system using CMG STARS equation: K = (k1/ P + k2 P + k3) exp(k 4 /(T k5)), where P and T are the pressure and temperature of the system in kPa and oC, respectively. k1 to k5 are the constants found through regression of the experimental data. Table 2 summarizes the K-values parameters for the Surmont bitumen/ propane system, which will be used in the simulation models. The temperature and pressure range of applicability for the saturated bitumen/propane K-values are: 100–190 °C and 1–8 MPa, respectively.

Solvent

Steam and solvent condensing Oil and condensate draining

3.3. Viscosity of bitumen/propane mixture

Fig. 3. (a) Schematic of typical solvent-aided process and (b) numerical simulation geometry of mechanistic model.

Bitumen viscosity was modeled using the correlation reported by Mehrotra and Svrcek [39]: ln(ln(µ )) = [c1 + c2 ln(T )] + c3 P , where µ is the viscosity of the dead bitumen in cP, P is the gauge pressure in MPa, T is the absolute temperature in K. The raw bitumen sample used is Surmont bitumen and its measured molecular weight is 540 g/g-mol. The pseudo liquid viscosity of solvent dissolving in bitumen was determined using the correlation reported by Zirrahi et al. [34]. It is assumed that the effective solvent viscosity in the oleic phase is a function of temperature and pressure of the system: µs = d1 + d2 T + d3 T 2 + (d 4 + d5 T ) P , where T and P are the absolute temperature and pressure in K and MPa, respectively. It is worth mentioning that this correlation is valid for prediction of the pseudo liquid viscosity of solvent for temperature and pressure ranges of 50 T 175 C and up to 4 MPa, respectively, and cannot be used to predict pure solvent viscosity. The pseudo-liquid viscosity correlation provides hypothetical viscosities for solvent (propane) to only get a good match for the binary mixture experimental data. Nonetheless, the trend of the estimated propane pseudo-liquid viscosities are physically sound and in agreement with actual pure propane viscosities. The mixture viscosity is found through employing a log mixing rule based on mole fractions: ln µm = xs ln µs + xb ln µb , where the subscripts m, s, and b denote mixture, solvent, and bitumen, respectively. Constants c1 to c3 and d1 to d5 are the fitting parameters which can be found through regressing the experimental viscosity data of raw bitumen and bitumen/propane mixture. The results obtained for the fitted parameters of the bitumen and effective solvent viscosity equations are summarized in Table 3.

4. Design of numerical simulation models The numerical simulations represent mechanistic study of solventaided thermal recovery processes. In these simulations, bitumen is exposed to heat and solvent from bottom as shown in Fig. 3a which results in reducing its viscosity and consequently increasing its mobility. Diluted bitumen will then flow downward by adverse density (gravity) effects. In order to understand the impact of non-equilibrium interphase mass transfer on solvent-aided thermal recovery processes, two-dimensional simulation models were constructed for bitumen/propane system using CMG STARS as shown in Fig. 3b. A 2D 1 m by 2 m mechanistic numerical model was constructed where the upper half of the model is the bitumen zone and the lower half of the model is the solvent zone. As solvent dissolves into bitumen, the mobilized bitumen will flow downward due to higher mobility and the adverse density (gravity) effects. There are no injection/production wells designed in the simulation model of this study. The numerical model is designed such a way that the process is governed only by gravity drainage. This was conducted intentionally to avoid interference of injection and production wells control on the performance of the process. The numerical model is assumed homogenous with constant properties. Table 4 summarizes the typical reservoir and rock properties used in the numerical model. Mechanistic studies have been conducted on the above model where pressure, temperature, and grid size are varied, and the interphase mass transfer is calculated from the simulation results. In order to study the impact of the non-equilibrium interphase mass transfer, two main cases are studied herein.

Table 2 Coefficients of K-values correlation obtained for saturated bitumen/propane system. k1 (MPa) 3.10 × 10

4

(b)

k2 (1/MPa)

k3

k4 (K)

k5 (K)

AARD (%)

0.000483

−2.267099

−279.601

−49.91377

5.5

4.1. Equilibrium case In the equilibrium case, phase distribution of components will reach 817

Fuel 241 (2019) 813–825

A. Al-Gawfi et al.

component i. The transfer term i can be expressed as first-order Fickian expansion term given by [9–11]:

Table 4 Summary of reservoir and thermal properties used in the numerical model. Parameters Reservoir properties

Thermal Properties

Value

Unit

Porosity Permeability Initial water saturation Mole fraction of Solvent

0.3 3000 0 0

Thermal Conductivity-Rock

2.47 × 105

ik

– mD – –

4

Thermal Conductivity-Water

5.35 × 10

Thermal Conductivity-Oil

1.15 × 104

Thermal Conductivity-Gas

4500

Volumetric Heat Capacity Formation Compressibility Over/Under burden capacity Over/Under burden compressibility

1.2 × 106 7 × 10−6 2.74 × 106 2.47 × 105

J/m/day/ C J/m/day/ C J/m/day/ C J/m/day/ C J/m3 1/kPa J/m3/day J/m/day/ C

Ki (P , T ) =

Solventliq

R2

(4)

Solventgas,

(5)

o

i

+ qio

V (Nion + 1 t

Nion) = 0

(6)

Tig yig

g

i

+ qig

V (Nign + 1 t

Nign ) = 0

(7)

yik ),

(8)

k1i k 4i + k 2i P + k3i exp . P (T k5i )

(9)

1 , xi

(10)

where x is the equilibrium mole fraction of the component i in each grid block. Therefore, the reaction term specified by keyword: RXEQFOR is: ik

= ±K (x ik

x ik ),

(11)

where K is the frequency factor for dissolution or exsolution of the component [1/day], x ik is the mole fraction of the component i in phase k in each grid block, and x ik is the solubility limit (equilibrium) of component i in phase k in each grid block. 4.3. Methodology of the numerical analysis and mechanistic studies The impact of non-equilibrium interphase mass transfer term is most prominent for slow processes. However, from a modelling perspective, for very fine grid size blocks, the non-equilibrium interphase phenomenon diminishes as the term (xik x ik ) in Eq. (11) becomes insignificant. Also, the dissolution/exsolution rate frequency factor approaches zero for a very fine grid size. Hence, theoretically, the assumption of instant equilibrium is valid when modelling using very fine grid size blocks. Further to note, it is computationally intensive and time consuming to model field scale cases with fine grid size blocks. Therefore, by including the dissolution/exsolution reaction terms when modeling field scale processes, it is possible to capture the non-equilibrium interphase mass transfer phenomenon and study its relative impact on oil recovery for proper design and implementation of solventbased processes. To begin with, several base cases for propane/bitumen system with fine grid size block (1 cm) are generated at different pressure and temperature conditions using the equilibrium K-values method (instant equilibrium assumption). Then, other cases with the equilibrium Kvalues are generated at coarse grid block sizes and compared to their respective base cases to quantify the impact of non-equilibrium mass transfer phenomenon. Afterward, non-equilibrium cases are generated using the rate-based reaction terms for dissolution and exsolution of solvent from bitumen. Through tuning the exsolution/dissolution frequency factors, a match with their respective base-cases at fine grid size blocks will be achieved. The frequency factors obtained intends to capture the necessary information pertaining to the non-equilibrium

where R1 and R2 represent the reaction rates of the dissolution and exsolution, respectively. The non-equilibrium mass transfer is modelled through the use of rate dependent reaction terms substituting the commonly used equilibrium K-values approach. Typically, the reservoir simulator solves a set of mass balance equations pertaining to a component in all grid blocks. For the sake of illustration, when mass transfer at grid scale is not instantaneous, mass balance equation for component i in each phase can be written as [19]:

Tio yio

alg (yik

Ki (P , T ) = k3i =

In the second set of simulations, the fluid model used is the same as the equilibrium simulations with the difference of adding another component to the fluid model to represent the solvent (propane) in the form of free gas and solution gas. Due to the high viscosity of bitumen, the dissolution/exsolution of solvent from bitumen is not instantaneous, and therefore, to simulate the kinetics of the dissolution/ exsolution of solvent from bitumen, CMG STARS allows for the definition of different forms of the gaseous solvent. Thus, for the nonequilibrium cases, propane is defined as solution gas (dissolvable in oil) and as a free gas. With the addition of the new component, it is now possible to simulate the transformation kinetics of the solvent through the use of reactions described as:

Solventliq,

ik

Since bitumen is a non-volatile pseudo-component, the solubility of a component is defined as the inverse of Ki (P , T ) . Therefore, by setting k1i , k 2i , k 4i , and k5i to be zero, Eq. (9) becomes:

4.2. Non-equilibrium case

R1

kD d ik

where Dik is the micro-scale Fickian diffusion coefficient for component i in phase k, k is the molar density of phase k, dk is the diffusion length of component i, yik is the solubility limit of component i in phase k, yik is the micro-scale mole fraction of component i in phase k, and alg is the specific interfacial area separating gas and liquid phases. In CMG STARS, it is possible to model the non-equilibrium interphase mass transfer term ik through the use of a rate-based reaction terms, where the term Dik alg / dik in Eq. (8) is lumped into one variable, K [1/s], which is called the dissolution/exsolution rate frequency factor. Also, it is worth noting here that the term Dik / dik in Eq. (8) is equivalent to the non-equilibrium interphase mass transfer coefficient, k [m/s]. The keyword “RXEQFOR” defines the chemical reaction that deviates from equilibrium described by:

their respective equilibrium limits in each grid instantly. Phase composition distributions are controlled by equilibrium K-values and calculated using the vapor–liquid equilibrium flash calculations. In CMG STARS, K-value are entered either in tabulated format or as a correlation. K-values correlation used in CMG STARS is presented in Section 3.2, and the parameters (k1 to k5) of the K-values correlations are determined through fitting the experimental data (see Table 2).

Solventgas



where Tio and Tig the transmissibility of component i in oil and gas phases, respectively, yio and yig are the mole fraction of component i in oil and gas phases, respectively, o and g are the flow potentials for both phases, qio and qig are the molar injection/production rates, V is grid block volume, Nio and Nig are moles of component i in oil and gas phases, respectively, and finally i is the interphase transfer term of 818

Fuel 241 (2019) 813–825

A. Al-Gawfi et al.

Fig. 4. Flowchart of mechanistic studies methodology.

mass transfer coefficients. The obtained mass transfer coefficients are used to develop correlations needed to mitigate the impact of nonequilibrium interphase mass transfer when modeling field scale solventbased processes. Fig. 4 shows the flowchart of the applied methodology. The base case is a 2D vertical cross section of 1 m by 2 m domain where the upper half of the model is the bitumen zone and the lower half of the model is the solvent zone. The grid block size in the bitumen zone is 1 cm in both x- and z-direction and the total number of grid blocks in the bitumen zone is 10,000. The solvent zone was gridded non-uniformly in the z-direction where at the interface the grid block size is 1 cm, then it was gradually increased in the solvent zone to reduce the run time of the simulations. The total grid block size of the solvent zone is 4100. Several base cases were conducted at different temperature and pressure settings. The temperature setting ranges between 50 and 100 °C, and pressure ranges between 1600 and 2500 kPa. 5. Simulation results

Fig. 6. Cumulative oil drained for propane/bitumen system, black line shows fine grid (1 cm) grid block size and red line shows 10 cm grid block size both using equilibrium K-values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 shows mole fraction of propane in oil phase for one of the base cases of propane/bitumen system at 1600 kPa and 75 °C at different times during the gravity drainage process. The results illustrate the dissolution of propane into bitumen as time progresses. As the base cases conducted using fine grid block size (1 cm), the physics of diffusion and convective mixing (density-driven fingering) are captured and observed as shown in Fig. 5. For illustration purposes, Fig. 6 shows the cumulative drainage of oil for the base case (1 cm) of propane/bitumen system at 1600 kPa and 75 °C using equilibrium K-values. The results show that dissolution of

t=0

solvent into bitumen results in mobilization and subsequent drainage of the bitumen. As mentioned before, all base cases are fine grids (1 cm grid sizes) and modelled using equilibrium K-values where a grid block is assumed to tend equilibrium condition instantaneously. The results for the 10 cm grid block size for the same temperature and pressure are also shown for comparison. As can be noted, the base case (1 cm) has a

t=5 days

t=3 days

Fig. 5. Propane/bitumen system base case at 1600 kPa and 75 °C – mole fraction of propane in oil phase at initial condition (t = 0), t = 3, and t = 5 days. 819

Fuel 241 (2019) 813–825

A. Al-Gawfi et al.

Fig. 7. Mole fraction of propane in oil phase at different times (initial condition, t = 3, and t = 5 days) for 1 cm (top row) and 10 cm (bottom row) grid block sizes at T = 75 °C and P = 1600 kPa.

higher cumulative drainage compared to the 10 cm grid size case. This difference can be attributed to several reasons among which is the nonequilibrium interphase mass transfer phenomenon. In other words, the non-equilibrium phenomenon can be utilized in coarse grid simulation to achieve the fine grid simulation results. In this case the impact in oil drainage is about 8%. This shows the necessity of including non-equilibrium interphase mass transfer phenomenon when modelling field scale (coarse grid) solvent-based processes. Fig. 7 shows a comparison between the mole fraction of propane in oil phase at different times for 1 cm grid block size as compared with 10 cm grid block size at T = 75 °C and P = 1600 kPa. The results show density-driven fingering is evident in the fine grid simulations. The results reveal that the coarse grid simulation is not able to capture the real physics of the problem. It is noted that a common choice of grid block in field scale numerical simulation of thermal recovery processes is 1 m. However, it is impractical to model field scale cases with such a high-resolution grid block size and that is why inclusion of non-equilibrium interphase mass transfer aim at capturing the details of fine grid simulations.

Fig. 8. Interphase mass transfer coefficient for propane/bitumen mixture obtained using ji = keff A i , from direct numerical simulations.

transfer coefficient, which will be discussed in the next section. It can be noted that the mass transfer coefficients for propane/bitumen system at different temperature and pressure settings is in the range of 1 × 10−8–4 × 10−6 m/s. The dashed trend line provides an average for the mass transfer coefficient obtained for propane/bitumen system. The upper and lower limits show the mass transfer coefficient for larger and smaller grid blocks, respectively.

5.1. Numerical determination of the equilibrium mass transfer coefficients for propane/bitumen system The results of interphase mass transfer coefficients for propane/bitumen mixture are presented in this section. 2D numerical simulations were conducted at different pressure and temperatures by considering three grid block sizes of 1, 5, and 10 cm and the mass transfer coefficients of the propane/bitumen system were determined using ji = keff A i , where ji is the mass flux of the component i [kg/s] obtained from direct numerical simulations, keff is the effective mass transfer coefficient [m/s], A is the area normal to flow [m2], and is the density of the mixture [kg/m3], and is the driving force mass fraction difference for component i. The results are shown in Fig. 8. The effect of temperature and pressure is captured through viscosity. The scatter observed is as a result of grid block dependency of the mass

5.2. Numerical simulation results of modelling non-equilibrium interphase mass transfer In this section, non-equilibrium interphase mass transfer is included in the numerical modelling by including the ex-solution and dissolution reaction term described in Section 4.2. CMG CMOST software [40] was used in which the coarse grid non-equilibrium cases were conducted to match the respective fine grid equilibrium base cases through tuning 820

Fuel 241 (2019) 813–825

A. Al-Gawfi et al.

cases. The error is calculated based on the cumulative drained oil at the end of the simulation run. In almost all cases, it is possible to tune the dissolution/exsolution frequency factors to obtain a good match to the base cases and reduce the error. In general, the assumption of instant equilibrium modelling of field scale solvent-based processes would result in error ranging between 4% and 10% in predicting oil recovery as shown for the case of propane/ bitumen system. This is for the situation when grid blocks are upscaled 5–10 times. However, further upscaling of the model may result in higher discrepancy in predicting oil recovery. Also, using other solvents with bitumen would result into different range of errors in oil recovery. It is also worth pointing that the last case shown in Table 5 above did not yield a better result when including non-equilibrium terms and that is likely due to numerical instability. Through analyzing the results of numerical simulations conducted for the propane/bitumen system at different pressure, temperature, and grid size settings, the scaling relations of dissolution rate frequency factors and interphase mass transfer coefficients were developed. The temperature setting ranges between 50 and 100 °C, pressure ranges between 1600 and 2500 kPa, and grid block sizes used in the development of the scaling relations are 5, 10, and 20 cm. Fig. 12 shows the frequency factor for dissolution (K ) versus grid block size. These scaling relations can be used in designing and modeling of propane/ bitumen processes while capturing the impact of non-equilibrium interphase mass transfer phenomenon. As observed from Fig. 12, the dissolution rate of propane into bitumen increases with grid block size regardless of the temperature and pressure of the system. Typically, for very fine grid block sizes, the dissolution rate values are insignificant and in such cases the equilibrium and non-equilibrium modelling of mass transfer cases considered to be similar. For large grid block sizes, however, the dissolution rate increases confirming the importance of the non-equilibrium interphase mass transfer when modelling field scale solvent-aided processes. As the dissolution rate by definition is the product of the interphase mass transfer coefficient with specific interfacial area (i.e. grid block dimension), non-equilibrium interphase mass transfer coefficient is directly proportional to dissolution rate frequency factor. Also, it can be observed from Fig. 12 that for large gird block sizes (1 m), the extrapolated dissolution rate tends to approach to a range of values between 8 and 10 (1/day) for the propane/bitumen system depending on the temperature and pressure condition of the system. It is important to note that the exsolution frequency factors did not exhibit a consistent trend with varying grid block sizes and it takes large values ranging 5.5 × 108–2.5 × 1018/day. It is worth noting that the process studied here includes both dissolution and exsolution kinetics. However, simulation cases seem to be more sensitive to dissolution kinetics compared to exsolution. Therefore, it can be inferred that

two parameters only: dissolution frequency factor and exsolution frequency factor. The objective function is defined to minimize the difference in oil rates of the coarse grid non-equilibrium and fine grid base cases. The design of experiment implemented has three temperature settings: 50 °C, 75 °C, and 100 °C, and three pressure settings: 1600 kPa, 2000 kPa, and 2500 kPa. The levels of grid block size variables selected for the experiment are: 5 cm, 10 cm, and 20 cm. In total, there were 50 experiments to be conducted for this design of experiment, 20 of which were conducted using CMOST software, the remaining were done manually as it was easy to achieve the match without the need to use CMOST software. As the fine grid simulations are computationally intensive, CMOST software was used to reduce the time required for optimization process. For the coarse grid cases, manual tuning of parameters was done to achieve the match following a systematic approach to avoid any bias that may arise in the obtained results. In the manually tuned cases, first step the dissolution frequency factor is adjusted (since it is the dominant factor) until a good match is achieved, then, as a second step, the exsolution factor is tuned to observe if it can provide further optimization to the case. If the exsolution frequency factor provided a better match in oil drainage, a second optimization trial is conducted in which steps one and two are repeated. The process is continued till no further optimization is achievable. The results shown in Figs. 9–11 provides a comparison between the volume of the oil drained for equilibrium and non-equilibrium cases at different pressures and temperatures. The base case, which is modelled using fine grid block size, represents the ideal case that captures the physics of dissolution of solvent into bitumen, density-driven fingering and convective mixing. When using larger grid blocks, mechanisms that are inherently fine scale processes such as density-driven fingering effects are lost. Therefore, field scale modelling using the assumption of instant equilibrium will lead to erroneous results. The equilibrium cases presented in Figs. 9–11 show the difference in oil recovery compared to the base cases. The non-equilibrium cases presented in these figures show the optimized results obtained using CMOST optimization through tuning of dissolution/exsolution rate frequency parameters. It is noted that when adjusting the dissolution/exsolution rate frequency factors for the non-equilibrium cases, it is possible to mitigate the impact of non-equilibrium interphase mass transfer. It is worth noting that the step-wise profile seen in the equilibrium cases of Figs. 9–11 is due to the larger grid size impact on oil drainage. As fluid properties change per simulation grid block, for coarse grid block size cases, it results in a clear step-wise oil drainage profile. All the cases discussed here are for propane/bitumen mixtures at different pressures and temperatures. Two grid block sizes of 5 cm and 10 cm are presented in the examples shown in Figs. 9–11. Table 5 compares the percent error in the cumulative drained oil for the equilibrium and non-equilibrium cases against their respective base

0.14 0.12

Cumulative Oil Drainage (m3)

Base Case (1 cm) Non-Equilibrium (5 cm) Equilibrium (5 cm)

0.10 0.08 0.06 0.04 0.02 0.00

0

10

20

30

40

0.14

Base Case (1 cm) Non-Equilibrium (10 cm) Equilibrium (10 cm)

0.12 0.10 0.08 0.06 0.04 0.02 0.00

0

10

20

30

Fig. 9. Cumulative oil drainage for propane/bitumen system at 1600 kPa and 50 °C, (a) comparison of fine grid (1 cm) with 5 cm grid block, (b) comparison of fine grid (1 cm) with 10 cm grid block. 821

Fuel 241 (2019) 813–825

A. Al-Gawfi et al.

0.25 0.20

Cumulative Oil Drainage (m3)

Base Case (1 cm) Non-Equilibrium (5 cm) Equilibrium (5 cm)

0.15 0.10 0.05 0.00

0

10

20

30

0.25 0.20 0.15 0.10 0.05 0.00

40

Base Case (1 cm) Non-Equilibrium (10 cm) Equilibrium (10 cm)

0

10

20

30

Fig. 10. Cumulative oil drainage for propane/bitumen system at 2000 kPa and 75 °C, (a) comparison of fine grid (1 cm) with 5 cm grid block, (b) comparison of fine grid (1 cm) with 10 cm grid block.

dissolution is the controlling kinetics and exsolution is instantaneous. This can be attributed to the fact that the process studied here is dissolution dominant. However, in the presence of pressure variation and other mechanisms the exsolution kinetics may play a more important role. This will be explored further in another study through applying the non-equilibrium approach to ES-SAGD process. While the maximum grid block size used in our simulations is 20 cm, the trend of dissolution frequency factor suggests that a constant frequency factor range of 8–10 (1/day) will be able to capture the fine scale processes as the grid block approaches the typically used block size in field scale simulation of solvent-aided recovery processes. The results presented in Fig. 12 serves as a guideline for history matching of field data by incorporation of dissolution kinetics into thermal simulations. A correlation for the dissolution rate frequency factor of the propane/bitumen system was developed based on the data obtained from the non-equilibrium simulation models. First, the temperature, pressure, grid size variables were transformed and normalized to −1 and + 1 range, and then a regression analysis was done using Essential Regression and Experimental Design for Chemist and Engineers tool [41] to find a correlation for the scaling relation of the dissolution rate frequency factor as a function of the temperature, pressure, and grid block size of the system. The correlation can be expressed as:

K

= 4.956 + 5.372 + 0.239

3

3.339

+ 2.588

3

0.459 0.286

+ 0.626

2

Table 5 Cumulative oil drainage % error – equilibrium versus non-equilibrium compared with the fine grid (base case). Case 1600 kPa 1600 kPa 2000 kPa 2000 kPa 1600 kPa 1600 kPa

(12)

where K is the dissolution frequency factor 1/day, is the scaled grid block size, is the scaled temperature, and is the scaled pressure. The 0.25 0.20

Cumulative Oil Drainage (m3)

Base Case (1 cm) Non-Equilibrium (5 cm) Equilibrium (5 cm)

0.15 0.10 0.05 0.00

0

10

20

30

50 °C 50 °C 75 °C 75 °C 75 °C 75 °C

% Error- Nonequilibrium

% Error- Equilibrium

5 10 5 10 5 10

4.6 2.0 6.2 3.9 0.5 9.8

10.1 10.2 7.4 8.6 4.8 8.1

predicted R-squared is 0.961, which shows how well the model can predict responses for new observations. The correlation is valid for propane/bitumen system for the following range of temperature and pressure: 50 °C to 100 °C and 1600 kPa to 2500 kPa, respectively. The temperature, pressure, and grid size parameters were normalized before conducting the regression. Therefore, in order to use the scaling relation presented above with dimensional values for temperature, pressure, and grid size, the following normalization relations must be used: = 0.1333( x ) 1.6667 , = 0.04T 3, = 0.0022P 4.5556 , and where , , and are for dimensionless scaled temperature, pressure, and grid size, respectively. T is for temperature in oC, P is for pressure in kPa, and x is for grid size in cm. Normalization of experimental data was used in order to standardize the range of input variables, accurately present the statistical significance of input variables, and mitigate any misleading influence of higher magnitude input variables.

2.750

2,

, , , , , ,

Grid size (cm)

0.25 0.20 0.15 0.10 0.05 0.00

40

Base Case (1 cm) Non-Equilibrium (10 cm) Equilibrium (10 cm)

0

10

20

30

Fig. 11. Cumulative oil drainage for propane/bitumen system at 1600 kPa and 75 °C, (a) comparison of fine grid (1 cm) with 5 cm grid block, (b) comparison of fine grid (1 cm) with 10 cm grid block. 822

Fuel 241 (2019) 813–825

A. Al-Gawfi et al.

Fig. 14. Interphase mass transfer coefficient for propane/bitumen mixture – pure diffusion versus full numerical solution. Fig. 12. Propane/bitumen system dissolution rate versus grid size for different temperatures and pressures.

grid block simulation models. Therefore, pure diffusion cases are simulated initially to establish the base cases at fine grid block sizes and then the non-equilibrium interphase mass transfer was introduced at larger grid block simulation cases to capture the true behaviour of the process and mitigate the errors. Fig. 14 illustrates that the analytical solution of the pure diffusive interphase mass transfer coefficients is in the range of 1 × 10−11–1 × 10−8 m/s whereas the numerical results shown to be in the range of 1 × 10−8–1 × 10−6 m/s. It presents a comparison between the interphase mass transfer obtained for the pure diffusive processes and the numerical solution, which considers other aforementioned mechanisms. These results are for propane diffusing into bitumen as a function of bitumen viscosity. The scaling relations shown in Fig. 14 finds applications in numerical simulations of field scale solvent-aided recovery processes when propane is utilized as a solvent. 6. Summary and conclusions This study presented the analysis of non-equilibrium interphase mass transfer phenomenon with applications to solvent-aided enhance oil recovery processes. The contributions can be categorized into the following: analytical determination of pure diffusion interphase mass transfer coefficient for several solvent/bitumen binary mixtures, quantifying the impact of the non-equilibrium interphase mass transfer on the performance of solvent-aided processes using propane as the working solvent, and finally presenting a correlation that can be used as a guideline for modeling the non-equilibrium interphase mass transfer for field scale simulations of solvent-aided processes. The non-equilibrium interphase mass transfer was incorporated into multiphase thermal simulation model using a kinetic approach in order to quantify the impact of non-equilibrium mass transfer on the performance of solvent-aided thermal recovery processes through the use of reaction terms. Mechanistic studies were conducted using commercial thermal reservoir simulator (STARS) [38] for the propane/bitumen mixture to determine scaling relations for the non-equilibrium interphase mass transfer coefficients. In the case of propane/bitumen systems when upscaling the models 5–10 times the assumption of instantaneous equilibrium introduces 4–10% error in oil recovery. The analytical solution of the pure interphase mass transfer coefficient for propane/bitumen system was found to be in the range of 1 × 10−11–1 × 10−8 m/s whereas the numerical results shown to be in the range of 1 × 10−8–1 × 10−6 m/s. The developed scaling relations find application in field scale modelling of the non-equilibrium mass transfer for properly designing and predicting the process performance.

Fig. 13. Mass transfer coefficient versus grid size for different temperatures and pressures for propane/bitumen system.

The interphase mass transfer coefficient can be calculated from the dissolution frequency factor obtained above. The interphase mass transfer is the product of the dissolution rate with the grid size block dimension. The results for the interphase mass transfer coefficient for the propane/bitumen mixture at different temperatures, pressures, and grid block sizes are shown in Fig. 13. The mass transfer coefficient versus the grid block size results in a linear relationship since the grid size is the dominant parameter. Similar to the dissolution rate behaviour, the importance of the interphase mass transfer is quite prominent at field scale models, which indicates the necessity of including non-equilibrium mass transfer when modeling field scale solvent-based processes. The interphase mass transfer coefficient obtained analytically in Section 2 of this study is for pure diffusion. Therefore, those analytically estimated results are considered the lower limiting values for diffusive mass transfer in solvent-aided processes. When comparing the analytical results of the pure diffusive interphase mass transfer coefficients with the numerical simulation counterparts, Fig. 14, the latter is much higher as it considers other important mechanisms influencing mass transfer such as convective mixing, density-driven fingering, and the inherent numerical dispersion. It demonstrates that pure diffusive mass transfer process is relatively insignificant and accordingly will not be sufficient to capture the non-equilibrium process particularly in large 823

Fuel 241 (2019) 813–825

A. Al-Gawfi et al.

Acknowledgments

Energy, Computer Modelling Group Ltd., ConocoPhillips Canada, Devon Canada Co, Foundation CMG, Husky Energy, Imperial Oil Limited, Japan Canada Oil Sands Limited, Nexen Energy ULC, N-Solv Co., PennWest Energy, Statoil Canada Ltd., Suncor Energy, and Total E &P Canada. The support of the Department of Chemical and Petroleum Engineering and the Schulich School of Engineering at the University of Calgary is also acknowledged.

The authors would like to thank four referees for their constructive comments. The authors also wish to express their appreciation for the financial support from the Natural Sciences and Engineering Research Council of Canada and all member companies of the SHARP Research Consortium: BP Canada Energy Group ULC, Brion Energy, Cenovus Appendix

To model the density of raw bitumen, the correlation presented in Zirrahi et al. (2017) can be used.

=

o exp(

(A1)

P ),

where o = 1 + 2 T + 3 and = 4 exp( 5 T ) . T and P represent temperature and pressure in K and MPa, respectively, and α1-α5 are constants. To find the pseudo liquid density of propane in the oleic phase, the following equation was used:

T 2,

=

0

(A2)

exp( P ), 3

where, o = 1 + 2 T + 3 , and = 4 + 5 T , where ρ is density in kg/m , T is temperature in K and P is pressure in MPa, respectively, and 1 5 n are constants. The mixture density is found by using the mole-fraction-based mixing rule: m = i = 1 xi ei , where m is the mixture density, and x i is the mole fraction, n is for number of components, and ei is the pseudo density of the component in the mixture which can be found using Eqs. (A1) and (A2) for bitumen and propane, respectively. The parameters 1 to 5 and 1 to 5 were found through regression of the experimental data of Surmont raw bitumen and Surmont bitumen/propane mixture densities. The regressed results of the parameters 1 to 5 and 1 to 5 for Surmont bitumen and the associated AARD % is presented in the Tables A1 and A2, respectively.

T2

Table A1 Regressed constants of density correlation Eq. (A1) for Surmont bitumen/propane system. 3

1(kg/m

)

1299.651

2 (kg/K

m3)

−1.17221

2

m3)

3 (kg/K

0.000734

4 (1/MPa)

5 (1/K)

8.99 × 10−5

0.005595

AARD (%) 0.0682

Table A2 Regressed constants of density correlation Eq. (A2) for Surmont bitumen/propane system. 3

1 (kg/m

)

−2000.57

2 (kg/K

m3)

15.20045

2

3 (kg/K

m3)

4 (1/MPa)

−0.01988

−0.22803

References [15]

[1] Canadian Association of Petroleum Producers. Canada’s Oil Sands 2016:57. doi:10. 1002/14356007.a26. [2] Canadian Association of Petroleum Producers. Crude Oil Forecast, Market and Transportation 2017:54. [3] International Energy Agency. Energy policies of IEA countries. Canada 2015. [4] IEA. World Energy Outlook 2017. Int Energy Agency 2017. doi:10.1016/03014215(73)90024-4. [5] Nasr TN, Isaacs EE. Process for enhancing hydrocarbon mobility using a steam additive. US 6,230,814 B1, 2002. doi:10.1145/634067.634234. [6] Acedemies C of C. Technological Prospects for Reducing the Environmental Footprint of Canadian Oil Sands. 2015. [7] Israel B. Using solvents in the oilsands, Pembina Institute 2017. http://www. pembina.org/blog/using-solvents-oilsands. [Accessed 12 March 2018]. [8] Faradonbeh MR. Mathematical Modeling of Steam-Solvent Gravity Drainage of Heavy Oil and Bitumen in Porous Media;2013. [9] Niessner J, Hassanizadeh SM. Non-equilibrium interphase heat and mass transfer during two-phase flow in porous media-Theoretical considerations and modeling. Adv Water Resour 2009;32:1756–66. https://doi.org/10.1016/j.advwatres.2009. 09.007. [10] Niessner J, Hassanizadeh SM. Modeling kinetic interphase mass transfer for twophase flow in porous media including fluid-fluid interfacial area. Transp Porous Media 2009;80:329–44. https://doi.org/10.1007/s11242-009-9358-5. [11] Niessner J, Hassanizadeh SM. Mass and heat transfer during two-phase flow in porous media – Theory and modeling. Mass Transf Multiph Syst Appl 2011. [12] Mayer AS, Hassanizadeh SM, Illangasekare TH, Javandel I, Jensen KH, Oostrom M, et al. Soil and groundwater contamination: nonaqueous phase liquids – principles and observations. Washington, DC: American Geophysical Union; 2005. [13] Hunt JR, Sitar N, Udell KS. Nonaqueous phase liquid transport and cleanup: 1. Analysis of mechanisms. Water Resour Res 1988;24:1247–58. https://doi.org/10. 1029/WR024i008p01247. [14] Nambi IM, Powers SE. Mass transfer correlations for nonaqueous phase liquid

[16] [17]

[18] [19] [20] [21] [22] [23]

[24] [25] [26]

824

5 (1/K)

0.000477

AARD (%) 0.4477

dissolution from regions with high initial saturations 2003;39. doi:10.1029/ 2001WR000667. Zhang H, Schwartz FW. Simulating the in situ oxidative treatment of chlorinated ethylenes by potassium permanganate. Water Resour Res 2000;36:3031–42. Imhoff PT, Jaffe PR, Pinder GF. An experimental study of complete dissolution of a nonaqueous phase liquid in saturated porous media. Water Resour Res 1993;30:307–20. https://doi.org/10.1029/93WR02675. Powers SE, Abriola LM, Weber WJ. An experimental investigation of nonaqueous phase liquid dissolution in saturated subsurface systems – Steady-state masstransfer rates. Water Resour Res 1992;28:2691–705. https://doi.org/10.1029/ 92WR00984. Miller CT, Poirier-McNeil MM, Mayer AS. Dissolution of trapped nonaqueous phase liquids: mass transfer characteristics. Water Resour Res 1990;26:2783–96. https:// doi.org/10.1029/WR026i011p02783. Nghiem LX, Sammon PH. A non-equilibrium equation-of-state compositional simulator. SPE Reserv Simul Symp 1997. Indrupsky IM, Lobanova OA. Modelling non-equilibrium phase behavior of hydrocarbon mixtures. Soc Pet Eng 2015;463:695–8. https://doi.org/10.1134/ S1028334X15070016. Lu XH, Ji YH, Liu HL. Non-equilibrium thermodynamics analysis and its application in interfacial mass transfer. Sci China Chem 2011;54:1659–66. https://doi.org/10. 1007/s11426-011-4308-9. Wu W, Wang P, Delshad M, Wang C, Pope GA, Sharma MM. Modeling non-equilibrium mass transfer effects for a gas condensate. Field Spe 1998. Wilkins MD, Abriola LM, Pennell KD. An experimental investigation of rate-limited nonaqueous phase liquid volatilization in unsaturated porous media: steady state mass transfer. Water Resour Res 1995;31:2159–72. https://doi.org/10.1029/ 95WR01677. Ivory J, Chang J, Coates R, Forshner K. Investigation of cyclic solvent injection process for heavy oil recovery. J Can Pet Technol 2010;49:22–33. https://doi.org/ 10.2118/140662-PA. Chang J, Ivory J. Field-scale simulation of cyclic solvent injection (CSI). J Can Pet Technol 2013;52:251–65. https://doi.org/10.2118/157804-PA. Soh Y, Rangriz-Shokri A, Babadagli T. Optimization of methane use in cyclic solvent

Fuel 241 (2019) 813–825

A. Al-Gawfi et al.

[27] [28] [29] [30]

[31] [32] [33]

injection for heavy-oil recovery after primary production through experimental and numerical studies. Fuel 2018;214:457–70. https://doi.org/10.1016/j.fuel.2017.11. 064. Jia P, Knorr KD, Imran M. Advanced numerical simulation of solvent vapour extraction (SVX) processes. SPE Heavy Oil Conf – Canada 2013;2013:1–20. https:// doi.org/10.2118/165522-MS. Nourozieh H, Ranjbar E, Kumar A. Modelling of non-condensable gas injection in SAGD process –Important mechanisms and their impact on field Scale simulation models. SPE Heavy Oil Conf – Canada 2015;2015:1–15. Saha S, Peng D. A non-equilibrium phase behaviour model for compositional reservoir simulation. Pet Soc Canada 1992. Rangriz Shokri A, Babadagli T. Feasibility assessment of heavy-oil recovery by CO2injection after cold production with sands: lab-to-field scale modeling considering non-equilibrium foamy oil behavior. Appl Energy 2017;205:615–25. https://doi.org/10.1016/j.apenergy.2017.08.029. Sheikha H, Pooladi-Darvish M, Mehrotra AK. Development of graphical methods for estimating the diffusivity coefficient of gases in bitumen from pressure-decay data. Energy Fuels 2005;19:2041–9. https://doi.org/10.1021/ef050057c. Hassanzadeh H, Pooladi-Darvish M, Keith DW. Modelling of convective mixing in CO2 storage. J Can Pet Technol 2005;44:43–50. https://doi.org/10.2118/05-10-04. Nourozieh H. Phase partitioning and thermo-physical properties of athabasca

bitumen/solvent mixtures; 2013. [34] Zirrahi M, Hassanzadeh H, Abedi J. Experimental and modelling studies of MacKay River bitumen and light n-alkane binaries. Can J Chem Eng 2017;95:1417–27. https://doi.org/10.1002/cjce.22775. [35] Azinfar B, Haddadnia A, Zirrahi M, Hassanzadeh H, Abedi J. A thermodynamic model to predict propane solubility in bitumen and heavy oil based on experimental fractionation and characterization. J Pet Sci Eng 2018;168:156–77. https://doi.org/ 10.1016/j.petrol.2018.04.065. [36] Poling BE, Prausnitz JM, O’Connell JP. The Properties of Gases and Liquids. New York: McGraw-Hill; 2001. 10.1300/J111v23n03_01. [37] Nuske P, Joekar-Niasar V, Helmig R. Non-equilibrium in multiphase multicomponent flow in porous media: an evaporation example. Int J Heat Mass Transf 2014;74:128–42. https://doi.org/10.1016/j.ijheatmasstransfer.2014.03.011. [38] Computer Modelling Group. STARS User Guide. Comput Model Gr Ltd; 2016. p. 1511. [39] Mehrotra AK, Svrcek WY. Viscosity of compressed athabasca bitumen. Can J Chem Eng 1986;64:844–7. https://doi.org/10.1002/cjce.5450640520. [40] Computer Modelling Group. CMOST User Guide. Comput Model Gr Ltd; 2016. [41] David D. Steppan, Joachim Werner RPY. Essential Regression and Experimental Design for Chemists and Engineers 1998. http://www.jowerner.homepage.t-online. de/ERPref.html#preface. [Accessed 20 August 2018].

825