Copyright © the author(s). This work is licensed under a Creative Commons Attribution 4.0 International License. 

DOI: 10.14800/IOGR.1184 

Received February 3, 2021; revised March 23, 2021; accepted April 2, 2021.  

*Corresponding author: neogi@mst.edu  
1 

Heavy Oil and VAPEX Process: A Brief 
Review 

 
Vijitha Mohan, University of Calgary, Calgary, Canada; Parthasakha Neogi* and Baojun Bai, 

Missouri University of Science and Technology, Rolla, USA 
 

 

Abstract 

Heavy oil recovery requires either heat or solvent assistance to reduce its high viscosity first, and then the less 

viscous oil can be recovered. Processes that use vapor assisted petroleum extraction (VAPEX) utilize vaporized 

solvent. Solvent vapors are introduced through a top injection well to induce viscosity reduction and the less 

viscous oil is collected at the bottom production well. It is noteworthy, that the solvent itself when recovered, 

can be seen as an integral part of oil. The exact molecular weight and chemistry of heavy oil are uncertain 

making analysis of thermodynamic properties difficult; however, this is overcome by making use of volume 

fractions of heavy oil and solvent used. It is possible to correlate all properties of heavy oil with or without 

solvent using the volume-based approach. To understand transport properties of heavy oil, studies have focused 

on free volume theory, which is also a volume-based approach. The diffusion of solvent in heavy oil is seen to 

be highly concentration-dependent. Both thermodynamic and transport properties using volume fractions, are 

discussed first. The efforts at introducing these properties to simulate production are then discussed. Also 

reviewed are work on sandpacks, which show good agreement with the theory when properly correlated. On 

the whole, it is found that the present day effectiveness of VAPEX has some important shortcomings. We 

mention some suggestions for improvement from literature. 

Introduction 

Highly viscous heavy oil (above 100 cp) cannot be recovered without the aid of an external resource like heat 

or solvents. The high viscosity provides impediment to easy recovery, making the reduction of the viscosity of 

oil a necessary precursor to oil recovery. Steam assisted gravity drainage (SAGD) process is a thermal process 

where steam is used as a heat source to reduce the viscosity of heavy oil to ~1cp. The less viscous oil flows 

under gravity to the drainage well. The water requirement, fuel usage to produce steam and discharge of a large 

amount of contaminated water, make the process, while still in use, expensive with future liabilities. Vapor 

assisted petroleum extraction (VAPEX) process uses gases (above critical temperature) or vapors (below 

critical temperature). These are usually called solvents even though they play the role of solutes in crude oil. 

They dissolve in the oil at the interface and diffuse into the bulk. In the process, the viscosity of the solution is 

brought down ~5cp, and it drains under gravity. The drainage is collected at the production well (Banerjee 

2012; Speight 2009).  

Quantitative understanding of the process is divided into two parts. In the first case, we look at the 

thermodynamic and transport properties of the oil-solvent mixture. In particular, the solubility has less variation, 

but the viscosity and diffusivity have very large variations with solvent concentrations. The solvent itself varies 

in its properties. At the reservoir temperatures, methane would be well above the critical temperature and would 

not condense. CO2 is below its critical temperature and will condense but appears to be immiscible with heavy 

oil even at very high pressures (Chung et al. 1988). Higher alkanes are below the critical temperature and will 

condense, and the condensate is miscible with oil. Prausnitz et al. (1999) have provided a list of specific 

mailto:neogi@mst.edu


 2 

volumes of mainly noncondensable gases in a liquid. They appear to be the same as those of organic liquids. 

Thus, the specific volume of CO2 at 25ºC is shown as 1.25 cm3/g.   

The difficulty in characterizing heavy oil is significant. It is a mixture of high molecular weight compounds 

that has a variety of structures. When we look at the molecular weight distribution, such fractionation does not 

imply that in a cut, the chemistry of oil molecules is the same. Similarly, if we fractionate the oil by chemical 

types, it does not mean that in a fraction the constituents have close molecular weights. In addition, different 

oil fields yield different oil types (Banerjee 2012; Speight 2009; Subramanian et al. 1996). All thermodynamic 

and transport properties are reported after they have been suitably weighted by the above distribution, where 

these weights change according to the property under consideration. 

Further along these lines, the higher molecular weight fraction contains a group of compounds called 

asphaltenes. These are primarily aromatic but quite complex in structures, often containing heavy metals, and 

vary not only from one molecule type to another but also from one oil field to another. Asphaltenes precipitate 

in presence of solvents of small molecular weights. The precipitates appear as small gelatinous masses, and 

the largest amount is seen with propane and the lowest with CO2. There has been considerable amount of work 

in the thermodynamics of asphaltene precipitation (Burke et al. 1990; Nghiem and Coombe 1997; Buckley et 

al. 1998). Asphaltene precipitation is not important in this work but will appear from time to time. The reason 

is that when the dynamics are studied for either measurement of transport properties or displacement process 

in the laboratory, the asphaltene precipitation is seen well after the experiments are over, that is, it appears to 

be a very slow process. 

When we look at the process itself, we exclude those with the same outlet well as inlet (Upreti et al. 2007), 

although many remarks of what we say for systems with separate injection and production wells also hold there. 

The schematic is shown in Figure 1. The original experiments by Butler and Mokrys (1989) in a vertical Hele-

Shaw cell showed small recovery. However, later work on two-dimensional sand packs showed much higher 

yields and since then has received the attention of many investigators. The yield is also low when CO2 is used. 

CO2 is still a suitable gas for sequestering in heavy oil reservoirs (Bachu and Shaw 2003; Shaw and Bachu 

2002). Small oil production can be enough to pay for the project. 

In sections below, we have looked at suitable correlations of both thermodynamic and transport properties, 

followed by models on the recovery process where these properties come into play. 

 

 

Figure 1—Schematic view of a VAPEX process. The height of the pay zone h is shown and a sharp solvent-

heavy oil interface is shown. 

Physical Properties of Oil with Solvent 

Physical properties in general can be divided into thermodynamic and transport properties. There have been 

good amount of measurements reported, but their correlations are scattered. As mentioned earlier, there is one 

Solvent vapor zone 

Solvent injection 

Interface 

Drained oil collection 

Oil zone 

h 



 3 

key difficulty in characterizing chemistry and molecular weight distribution of the components in heavy oil; 

hence, how to arrive at properly weighted physical properties continues to be debated. One feature of heavy 

oil is that the molecular weights of all molecules are large, a situation found in polymers where individual 

molecules are all large even though a molecular weight distribution exists. There can also be chemical 

differences among polymer molecules, such as chain branching. Flory (1953) has shown that in polymer 

solutions, the thermodynamic properties should be expressed in terms of volume fractions rather than mole 

fractions. In crude oils, the specific chemistry of oil is bundled under the name of the oil field from which it 

was extracted. Flory had problems in getting his point of view accepted, and one editor turned down request 

of our manuscript for review with a comment that they only look at materials with fully specified chemistry. 

As discussed later, the same volume based approached can be taken to understand the transport properties. 

Thermodynamic Properties 

One key assumption made here is that, for an oil-solvent pair, their volumes are assumed to be additive. When 

the solvent is CO2, it becomes difficult to determine what its specific volume will be in oil. In this case, it is 

back calculated from the data. Chung et al. (1988) have reported extensive data on Bartlett heavy crude and 

CO2. Under additivity of volumes, 

𝜌 = 𝜌𝑠𝜙 + 𝜌𝑜(1 − 𝜙),……………………………….………..…………………………………..……….(1) 

and  

𝑓𝑇 = 𝑐𝜙 + 𝑓(1 − 𝜙),………………….…………………………..………………………….……..……..(2) 

where ρ is the density, the volume fraction of the solvent and the subscripts s and o stand for solvent and oil. 

Now, oil and solvent can be split into the hard molecular dimensions and free volume. Here, f and c are the 

free volume fractions of oil and CO2 and the subscript T stands for total. Chung et al. (1988) provide the 

densities of oil, with or without CO2, at three different temperatures, and pressures going up to 5000 psi. These 

provided Tran et al. (2012) with a way of calculating f and c and the pressure and temperature dependence of 

f (namely, through the isothermal compressibility and the thermal expansion) but the data were sensitive 

enough to obtain the temperature dependence of c only (and not the pressure dependence). The free volume c 

was about four times the value of f.  

In CO2 oil recovery, CO2 swells oil and this feature is thought to be important in squeezing oil out of narrow 

pores and crevices. Tran et al. (2012) found that dissolved CO2 has a specific volume similar to organic liquids, 

and this volume adds to the volume of oil. The resulting swelling gave then the Welker and Dunlop (1963) 

correlation for swelling in Day crude, Kansas: 

𝑆𝐹 = 1 + 3.5 × 10−4𝑆,……………………….………………....…………………………………..…….(3) 

In Eq. 3 the swelling factor SF=volume occupied by oil containing CO2/(volume occupied by the same mass 

of oil without CO2 at 1 atm and the same temperature) and S=standard cubic feet of CO2/bbl. The result of 

Welker and Dunlop (1963) has been checked out with a number of other oils by Chung et al. (1988). In addition, 

Eq. 3 led them to the specific volume of CO2 in oil of 1.06 cm3/g, where Prausnitz et al. (1999) suggested 1.25 

cm3/g.  

Another discovery of Tran et al. (2012) was that the CO2 solubility data in Chung et al. (1988) fitted Henry’s 

law well. Henry law coefficients themselves at different temperatures showed an activation energy. 

Furthermore, for the subcritical case (temperature below Tc) Henry’s law held even during the phase change 

of CO2 from vapor to liquid. Due to the large types of solvents available, as discussed earlier, Mohan et al. 

(2017a) performed a volume fraction-based analysis of solubility following Flory-Huggins theory (1953) and 

the data on Lloydminster crude by Yang and Gu (2006a). The additivity of volumes was assumed. Supercritical 

gases, methane and CO2 which condensed but the condensate was not miscible with oil, both showed constant 

Henry’s law type behavior. However, the condensable higher alkanes did not, particularly when the saturation 

vapor pressure was approached. There are two parameters in the system. The first is x, the ratio between the 

size of an average oil molecule to that of a solvent molecule, which was found to be very high for all the above 



 4 

solvents. It can be set to infinity. The second is χ, Flory-Huggins coefficient, which was found to be greater 

than ½ for ethane and propane. This crosses the limit of solution stability, and it was assumed that asphaltene 

will precipitate. 

Finally, the surface tension of CO2-oil solution exposed to pure CO2 has been modeled (Tran et al. 2012). 

There are two aspects to the chemical potential of CO2 at the interface. The first is the usual mixing entropy, 

which can be modeled using simplified Flory theory using volume fractions. The second is the surface work, 

which uses work done by surface tension over a reference state (Miller and Neogi 2008). The results fit the 

data of Rojas and Ali (1988) of Aberfeldy heavy oil, Alberta, Canada, quite well. The model comes with three 

parameters, all of which have physical interpretations.  

It should be emphasized here that, in the correlations discussed so far, there is no need to know the molecular 

weight or molecular weight distribution, or information on species distribution. Like in polymers, volume 

fractions are all we need. Further, the assumption of additivity of volumes does not lead to problems where we 

end up with results that cannot be fitted to the data, or on fitting yields unphysical values of parameters. 

Transport Properties 

The diffusivity of a solvent in heavy oil is very concentration dependent. When the solvent enters the oil, the 

leading toe which is at very low solvent concentration shows a very low diffusivity. It severely limits the rate 

of penetration. Further, the solvent also decreases the viscosity of oil. Thus, the viscosity of oil fell to 2% and 

to 10%, both for CO2 in heavy oil (Welker and Dunlop 1963; Chung et al. 1988). Consequently, it is not just 

the viscosity, but also its concentration dependence, that is of importance, and the lack of penetration. 

Considered below are some theories of diffusivity and the difficulties in measuring those. 

When the liquid under consideration has a high density in its pure form, its free volume (total volume less 

the hard volumes of the molecules) is low. Generally, when a solvent molecule or a segment of oil moves, it 

has to overcome an Arrhenius type of activation energy barrier. As a result, both diffusivity and viscosity show 

such temperature dependence. However, when the free volume is low, the available space becomes restricted, 

and this restriction controls mobility. Fujita (1961) argued that as the free volume of the solvent was much 

larger (earlier we had said that c was four times larger than f for CO2 in heavy oil) the free volume grows as 

the solvent enters oil and viscosity of the solution drops.  

The free volume theory provides the expression for viscosity of the form (Cohen and Turnbull 1959), 

𝜇 = 𝑅𝑇𝐴𝜇𝑒
𝐵𝜇

𝑓𝑇 ,…….……………………………………..………..………………………………...……..(4) 

where Aμ and Bμ are constants and fT is given by Eq.2. Eq.4 can be rewritten for convenience as  

𝜇 = 𝜇𝑜𝑒
𝐵𝜇(𝑓−𝑐)𝜙

𝑓.𝑓𝑇 ,……….…………………………...……………...……………………………….………(5) 

where μo is the viscosity of the pure oil. Since c is larger than f, the argument in the exponential of Eq. 5 is 

negative. 

Tran et al. (2012) compared the extensive data by Chung et al. (1988) on Bartlett oil to Eq. 5. At the end, 

after all parameters were obtained, they could actually predict the viscosity at one temperature as a function of 

CO2 partial pressure which showed excellent comparison with the actual data. Two fitted results are of interest, 

Bμ = 1.000 and f* the reference free volume fraction in pure oil at 1 atm and 75ºF is 0.02088. In general, the 

free volume fraction of oil decreases with increasing pressure and increases with temperature. 

Under free volume theory diffusivity becomes (Vrentas and Duda 1979)  

𝐷 = 𝐷𝑜𝑒
𝐵𝑑(𝑐−𝑓)𝜙

𝑓.𝑓𝑇 ,………………………………..…..……………...………………..…………….….……(6) 

Note that this time the argument is positive and Bd is generally set to Bµ. There are a large number of 

diffusivity data of solvents in heavy oil. The difficulty with these is that to obtain diffusivity, one needs the 

solution to a boundary value problem where the diffusivity is a strong function of solvent concentration. One 

such solution exists (Neogi 1988). If a mathematical solution is obtained by assuming that the diffusivity is a 



 5 

constant, it leads to a concentration averaged diffusivity which is not possible to interpret and difficult to use. 

In addition, there are some items in solving a boundary value problem that need to be addressed. To write the 

diffusive flux as 

𝑗𝐴
𝑜 = −𝐷∇𝜌𝐴 = −

𝐷

𝑣𝐴
∇𝜙,………………..…………..……………..………...………….……….……….. (7) 

where 𝜌𝐴and 𝑣𝐴mass concentration and specific volume of the diffusing species A, it is necessary to use 

volume average velocity, and assume that the specific volumes are constants, that is, volumes are additive 

(Bird et al. 2002) . In general, there will be a convective flux (from the volume average velocity above) in 

addition to Eq. 7, whether or not convection is being forced. Consider also a special case where there is a drop 

of heavy oil exposed to a vapor or half-filled test tube with oil exposed to vapor. The vapor will enter the oil 

through dissolution and diffusion and swell it. If the volume of oil under consideration is increasing, then that 

alone will give rise to a convective term. When the problem under consideration is a moving boundary problem, 

the boundary condition at the liquid-vapor interface becomes a jump boundary condition. Slattery (1972) has 

provided proofs of the above. As mentioned earlier, there are a large number of experimental works, a review 

by Ghanavati et al. (2014), and papers by Afsahi and Kantas (2007), Fedai et al. (2013), Zhang et al. (2007), 

Yang and Gu (2006a, 2006b, 2008) all of which contain the above errors in some form or another. However, 

we do not know the magnitude of the error.  

Returning to free-volume theory and Eq. 6, if numbers are used there, then Eq. 6 approximates as  

𝐷 = 𝐷𝑜𝑒𝛼𝜙,…………….………………...….……………………………..………………………..……..(8) 

where 𝛼 ≈
𝐵𝑑

𝑓
, which is property of the oil. Tran et al. (2012) found �̄� = 𝛼𝜙𝑜~ 10 for CO2, and about the same 

value was found by Mohan et al. (2017b) for hexane, heptane and toluene in heavy oil. Here, 𝜙𝑜is the solubility 

of the solvent in oil. One peculiarity of Eq. 8 is that at moderate solvent concentrations not much variation in 

diffusivity with concentration is observed with diffusivities ~10-5 cm2/s. However, it plunges sharply to low 

values to ~10-9 cm2/s at less than 1% solvent concentrations. This is also seen in solid polymers, where free 

volume theory holds (Vrentas and Duda 1979). The consequence is that when the solvent enters oil, the toe is 

at this small concentration and the small diffusivity controls the rate of penetration into oil. It leads to self-

sharpening of the solvent profile which gives rise to a pseudo-interface, sometimes called a concentration shock.  

Instead of looking at the concentration profile in this pseudo-shock region, Mohan et al. (2017b) monitored 

the location of the shock front as a function of time. The changes in the location of the interface is  

Δ𝑥𝑜 = [
2𝐷𝑜(𝑒𝛼𝜙𝑜−1)𝑡

𝛼(𝜙𝑜−𝜙𝑖)
]

1/2

,………………………..…………………..……...…………………….….……..(9) 

where 𝜙𝑜 = 1 as pure solvent was introduced and the initial concentration in oil was 𝜙𝑖 = 0. They contacted 

heavy oil with solvent (hexane, heptane and toluene) for that purpose and followed the front shown in Figure 

2. They supplemented this set of experiments with desorption with the object that between the two experiments 

they would get the two parameters, 𝐷𝑜and �̄� = 𝛼𝜙𝑜. However, desorption was too slow to use and an estimate 

of 𝐷𝑜was obtained from Stokes-Einstein equation.  

Note that as in Eq. 8, Eq. 5 yields 

𝜇 = 𝜇𝑜𝑒−𝛼𝜙,……………………………………..……………………...………………………….….…..(10) 

where viscosity decreases with solvent concentration. Like in the Stokes-Einstein equation, Eq. 8 and Eq. 10 

share exact inverse relation.  

If only Arrhenius type behavior is expected, 𝐷 ∝ 𝜇−2/3 for very viscous liquids (Hiss and Cussler 1973). In 

many cases, the entire concentration profile is measured. However, low solute concentrations are not usually 

reached. Fadaei et al. (2013) using light absorption and Ghanavati et al. (2014) using x-ray absorption go down 

to 10 wt % solvent but not less. They find, with some exceptions, diffusivities ~10-5 cm2/s as mentioned earlier. 

At lower solvent concentrations than above, Afsahi and Kantas (2007) reach D~10-7 cm2/s using NMR, that is, 



 6 

in the right direction but still a long way to go. Where needed, the Arrhenius equation is used for viscosity 

(Chun et al. 1988). 

 

 

Figure 2—A case where liquid hexane is contacted with heavy oil ( A Hauser, Kansas) is shown by Mohan et al. 

(2017b) with permission. 

Process Design 

The first to be considered are the recovery and the mechanism of oil production. Before offering results on oil 

recovery in processes in a porous medium as in Figure 1, we first look at displacement inside a single pore 

modeled as a cylinder. Tran et al. (2015) observed that the penetration of CO2 into oil was extremely low under 

dynamic conditions. Consequently, the line of thought that dissolved CO2 penetrated the oil and decreased the 

oil viscosity does not quite work. However, the displacement velocity seemed relatively unaffected when the 

viscosity of the liquid increased. A thin film of oil was left clinging to the wall during displacement, which 

probably offered much less resistance to the displacing CO2 than the bare solid wall. So, the conclusion was 

that heavy oil to very heavy oil does not change the displacement process and remain virtually unaffected by 

carbonation, actually by the lack of it. Note that the velocity of displacement front in an oil reservoir medium 

is same as that in a single pore.  

One other feature in mass transfer in porous media is that we should be using dispersivity (Gardner et al. 

1981). However, when we look at interphase mass transfer, it shows that dispersivity which is an averaged 

value cannot be used and molecular diffusivity will have to be used (Tran et al. 2016). Further, Lake (1989) 

has shown that the velocity associated with oil displacement (~1ft/day) is so low that the dispersivity is equal 

to the diffusivity. Consequently, molecular diffusivity is all we need for simulation. 

The convective-diffusive transport in an oil reservoir of the configuration in Figure 1 was first given by 

Butler and Mokrys (1989). It provides the recovery rate of 

𝑄𝑏 = √2𝑘𝑔ΦΔ𝑆𝑜ℎ𝑁𝑠,……………………………….……………….……………….…..………………(11) 

where 𝛥𝑆𝑜is the difference between the fractional pore volume containing oil before and after displacement. 

For perfect displacement 𝛥𝑆𝑜= 1. Others are k permeability, g acceleration due to gravity, Φ porosity and h is 

the pay zone height. Further, 𝑁𝑠 = ∫
𝛥𝜌𝐷(1−𝜙)

𝜇

1∫ 𝑙𝑛 𝜙

𝜙𝑚𝑖𝑛
 where Δ𝜌 is the density difference between pure oil less 

the solvated oil.  ϕmin is the solvent concentration at the end of the shock front of the solvent that has penetrated 

the oil. Finally, h is the height of the pay zone. 

Mohan et al. (2019a) found the concentration profile assumed by Butler and Mokrys (1989) to be 

inconsistent with the conservation of solvent equation.  However, their final results are not that different even 

after the concentration dependence of the transport properties as described above, are used. These are 



 7 

𝑄𝑏 = 2√
𝑘ΦΔ𝜌𝑔𝐷𝑜ℎ

𝜇𝑜
∫ 𝑒𝛼𝜙𝑜𝜓𝑢

0
(1 − 𝜙𝑜𝜓)𝑑𝜓,………...……………………………....………..………….. (12) 

where Δ𝜌is the density difference between solvated oil and solvent (vapor). Here,𝑢 =
𝛿

√4𝜇𝑜𝐷𝑜𝜂𝛷/𝑘𝛥𝜌𝑔 𝑠𝑖𝑛 𝜃
 

where δ is the thickness over which solvent concentration in oil falls from ϕo to ϕmin in the direction normal to 

the interface and 𝜓 = 𝜙/𝜙𝑜. This definition of the boundary layer and the profile itself are used to obtain the 

value of u. If we use a vapor which has a condensate that is miscible in oil, then the ϕo=1. The resulting drainage 

is calculated in this case to be mainly the solvent. For most solvents, known to us �̄� = 𝛼𝜙𝑜~10. This case 

corresponds to low solvent penetration and low oil recovery in general. The region occupied by can be divided 

roughly into a dry part and a thin sliver of oil at the interface with very high amount of solvent. The end result 

is that if free volume theory and 𝐷 ∝ 1/𝜇are used then both production and fraction oil at the effluent end will 

be very low, perhaps unacceptably low. If Arrhenius equation and 𝐷 ∝ 1/𝜇2/3are used, then there will be a 

small amount of improvement. Finally, if D~10-5 cm2/s and Arrhenius equation are used then both quantities 

will improve to possibly acceptable results. 

The second feature is that there are many sandpack results but they are not fully organized. Butler and 

Mokrys (1989) conducted experiments in Hele-Shaw cells and verified k and h dependence. Many other 

experiments using two-dimensional sand packs followed and showed often a result that was proportional to h 

rather than ℎ
1/2

as in above. Nenninger and Dunn (2008) correlated results of many experiments and showed 

that the recovery was given by mass flux (or mass velocity) 

𝑊𝑏 = 43550(𝑘Φ/𝜇𝑜)0.51,. …………………………..…………………….……………..……………(13) 

where 𝑊𝑏 = 𝑄𝑏𝜌𝑜/ℎ. Note that the variables on the right appear in the same manner as in Eq. 12. Eq. 13 is in 

a form of results of dimensional analysis. For a more formal dimensional analysis, if we take 𝑊𝑏to be a function 

of other variables 𝜇𝑜 , 𝐷𝑜 , 𝛼, 𝛷, 𝛥𝜌, 𝑔, 𝑘, 𝜙𝑜 , ℎ then there are 7 variables and 3 dimensions (mass, length and 

time) leading to 4 dimensionless groups (leaving aside the quantities that are already dimensionless) according 

to Buckingham-pi theorem (Buckingham 1914). These groups are found to be Reynolds’ number 𝑅𝑒 =

𝑊𝑏ℎ/𝜇𝑜the ratio between the inertial and the viscous forces, Froude number 𝐹𝑟 =
𝛷𝑘𝛥𝜌

𝜇𝑜
√

𝑔

ℎ
 the square-root of 

the ratio between kinetic energy and potential energy due to gravity, a square of aspect ratio 𝐴𝑟 = 𝛷𝑘/ℎ2
, and 

Schmidt’s number 𝑆𝑐 =
𝜇𝑜

𝛥𝜌𝐷𝑜
. Because of the way Φ is associated with k in Eq. 13, Φ has been considered only 

in a product Φk. It is reasonable to assume that in most cases 𝛼𝜙𝑜~ 10 and it is not considered to be a variable 

(Mohan et al. 2017b). Do that occurs is difficult to find, and Stokes-Einstein’s equation is used to calculate this 

value (Mohan et al. 2017b). Note that because of the inverse relation between viscosity and diffusivity in 

Stokes-Einstein equation, 𝐷𝑜 = 𝐵/𝜇𝑜where B is a constant. 

Fr/Re has been plotted against √(Ar√.Sc) using the data compiled by Nenninger and Dunn (2008) in Figure 

3a for sandpacks and fit to 

𝐹𝑟/ 𝑅𝑒 = 4 × 10−6𝐴𝑟1/2𝑆𝑐1/4 ,………..……..…....……………..….……………………..……..….. (14a) 

There are 8 outliers in their group of 43 on sandpacks that have been omitted. The fit has been stretched in 

Figure 3b to show that the fit is excellent at small values near the origin. This is good, because field results 

will have values of Ar close to zero on this scale. The data on Hele-Shaw had scatter that was too high and 

were dropped. Eq. 14a can also be re-expressed as 

𝑊𝑏 = 25 × 104√(𝑘𝛷/𝜇𝑜)
(Δ𝜌)5/4𝑔1/2𝐵1/4

ℎ
1/2 ,………..…………….……………..……………..……..…..(14b) 

or, in terms of Qb, 

𝑄𝑏 = 25 × 104√(𝑘𝛷/𝜇𝑜)(Δ𝜌)1/4𝑔1/2𝐵1/4ℎ1/2,………………..………………………..……..……...(14c) 

which supports the result that 𝑄𝑏 ∝ ℎ
1/2

 (Mohan et al. 2019b). 

The above theoretical results assume that flow is given by Darcy’s law, 



 8 

𝑣 = −
𝑘

𝜇
𝛻𝑃,……………………………………………….………………………………..………..…….(15) 

where P is the dynamic pressure, a sum of static pressure p and the potential energy due to gravity. Brinkman 

equation is (Brinkman 1947) 

∇[𝜇∇𝑣] −
𝜇𝑣

𝑘
− ∇𝑃 = 0,……………………………………………………………………..………..…..(16) 

Quantities in bold in the above equations are vectors. Note that if the first term in Eq. 16 is omitted, Eq. 15 

results. The velocity here is the superficial velocity. However, in conservation of species equation, the 

interstitial velocity 𝒗/𝛷 is used. 

When flow profile is introduced, Mohan et al. (2019a) used an approximation to show that 𝑄𝑏could reach a 

height dependence of h but only at larger values of k. Cuthiell and Edmunds (2012) using Tetrad, have provided 

the results of their simulation that shows 𝑄𝑏is proportional to h. It is concluded that as sandpacks offer large 

values of k only, those result may well support h instead of h1/2. The vapor oil-interface is not sharp, and 

Cuthiell and Edmunds (2012) put together a model expression for fractional flow as a function of fractional 

volume occupied by the oil, to obtain a diffuse interface. Mohan et al. (2019a) obtained an approximate shape 

of the vapor-oil interface as linear except at very large times and near the top end.  

Finally, there is the issue of stability of the interface. Tran et al. (2016) have shown that when displacement 

is accompanied by mass transfer, the interface is stable to disturbances of large wavelengths but unstable to 

disturbances of small wavelengths. Mass transfer and accompanied reduction of viscosity of the oil play 

important roles. If there is no mass transfer, the situation is reversed where disturbances of large wavelengths 

are unstable and those at small wavelengths are stable. There are some differences between the above stability 

results and the present case in that the directions of flow here are mainly parallel to the interface, and the 

displacement velocities are conventionally about 1 ft/day which appears to be a lot larger than in the present 

case.  

 

       
(a)                                                                        (b) 

Figure 3—(a) Plot of Fr/Re against √𝑨𝒓√𝑺𝒄 in sandpacks from the compilation by Nenninger and Dunn (2008). 

(b) The same data as in Figure 3(a) replotted differently to show the region near the origin. From Mohan et al. 

(2019). 

On Solvent Species 

In the immiscible displacement study of heavy oil in a sand pack with CO2, Rojas and Ali (1988) found that 

the reduction of surface tension was the most important reason for improved oil recovery. The data compiled 

by Nenninger and Dunn (2008) contain measured values of surface tension where the vapor phase is almost 

saturated. The values are often low. Since the vapor phase is posed to condense, the question arises as to what 

the dynamics of an oil drop floating on surface the solvent liquid will look like. Heavy oil when contacted with 

0.0E+00

5.0E-06

1.0E-05

1.5E-05

2.0E-05

2.5E-05

3.0E-05

3.5E-05

0 2 4 6

F
r/

R
e

√(𝐴𝑟√𝑆𝑐) 

0.00E+00

5.00E-06

1.00E-05

1.50E-05

2.00E-05

2.50E-05

3.00E-05

0.01 0.1 1 10

F
r/

R
e

√(Ar√Sc)



 9 

a solvent, shows a very clear interface as shown in Figure 4a and 4b. The interface looks clean but 

subsequently deformation show a wavelike character that suggests Rayleigh-Taylor instability due to gravity 

imbalance, Figure 4a, observed when the liquid on the top is heavier than the liquid at the bottom (Miller and 

Neogi 2008). The wavelengths associated with this instability should be of the order of basic dimension of the 

drop ~0.5 mm. However, one wave very quickly develops which is pointed and sharp, Figure 4b. This “wave” 

would correspond to that of zero wavelength and zero surface tension. In general, there can be a dynamic 

surface tension, that is small but positive, but quickly disappears to the equilibrium value of zero which is the 

accepted value for miscible systems. All photographs were taken almost instantaneously after layering the drop 

of oil on the surface of the solvent. 

 

 
 

  
(a)                                                                     (b) 

Figure 4—Inverted drop of heavy oil on glass substrate at 30oC flooded with (a) hexane and heptane, where 

heptane shows a pointed end corresponding to zero surface tension for Rayleigh-Taylor instability and (b) 

shows the Rayleigh-Taylor waves 

 

Eq. 12 shows that the choice of the solvent species is felt through 𝐷𝑜, 𝜙𝑜 and 𝛼𝜙𝑜. There are practically no 

data on diffusivity at infinite dilution in species of interest, however, the estimates show that there is not much 

difference among species (Mohan et al. 2017b). Similar estimates show that usually 𝛼𝜙𝑜~10 (Tran et al. 2012; 

Mohan et al. 2017b). Further, as the temperatures of operation are below the critical temperatures of solvents 

used, the vapor will condense and be miscible with oil. Thus, in that case 𝜙𝑜 = 1. Thus, one does not expect 

any differences in performance among species in the first approximation. This is observed by Nenninger and 

Dunn (2008) who have surveyed results for CO2, ethane, propane and butane. However, one point needs to be 

cleared up. Activity is defined as the pressure of the vapor p/saturated vapor pressure pv. In most cases the data 

tabulated by Nenninger and Dunn (2008) show activity of ~1. However, there are some which are below to 

well below that activity and are expected to show lower values of 𝜙𝑜 and 𝛼𝜙𝑜. To the first approximation, 

activity ~𝜙𝑜, yet even at below 0.5, this CO2 data show no significant deviation from Nenninger and Dunn 

(2008) correlation, Eq.12. On the other hand, the outliers to Nenninger and Dunn (2008) correlation have 

activity below 1 but not so much below as 0.5. It is necessary to go through experiments which attempt to 

correlate Qb to the activity of the vapor. Such data do not exist at present. 

None of these show asphaltene precipitation which according to Mohan et al. (2017b) is very slow. They 

also discuss the bounds of operation for an extremely fast rate of precipitation and for extremely slow rate of 

precipitation. 

If steam is used (in conjunction with vapor) the end results could be oil-in-water (o/w) or water-in-oil (w/o) 

emulsion. This issue has been addressed by Abdulmohsein et al. (2021) using sodium dodecyl sulfate or SDS, 

an anionic surfactant, cetyl trialkyl ammonium bromide, CTAB a cationic surfactant, and IGEPAL CO-530, a 

nonionic surfactant C9E6 close to the hydrophile-lipophile balance. The emulsion samples, viewed under 

microscope of o/w emulsions are shown in Figure 5. Experiments also had SiO2 or Al2O3 nanoparticles. The 

systems containing the nonionic surfactant were difficult to flocculate in presence of nanoparticles, even on 

raising temperatures when the surface active property of nonionics, disappear. In general, the phase separation 



 10 

was very quick, but some water was left in the oil and did not separate in a day’s time, because of high viscosity 

of oil. 

 

                         

Figure 5—Confocal microscopy images of o/w emulsion using SDS on the left and a bicontinuous structure using 

IGEPAL CO-530 at high oil: water ratio on the right. Oil is fluorescent and water appears black. 

Discussion 

We were not able to find examples of VAPEX process used in industry. It appears that the fall in price of oil 

has made this technology into one that is unsuitable for a long time, but trade journals suggest that view is 

changing (Rassenfoss 2017). One issue is where one would get the solvent from, and one interesting feature 

here is an unconventional oil source of shale oil, which on pyrolysis provides naphtha, typically 5-6 or 6-12 

carbon aliphatics (Beckwith 2012).  

If we look at the process itself, the correlation by Nenninger and Dunn (2008), Eq. 13 speaks a volume that 

higher the oil viscosity, the output will be lower. Further, the formation cannot be “tight.” One key feature that 

has been obtained is the fact that solvents like propane and butane would have a very high value of 𝛼𝜙𝑜 ~10. 

If we use CO2 at a very high pressure of ~3000 psi, there too, we will see high value of 𝛼𝜙𝑜. Solvents with 

high values of 𝛼𝜙𝑜do not penetrate the oil much hence give rise to low recovery. Consider now CO2 where the 

pressure is lowered. The solubility 𝜙𝑜 is proportionately lowered. Now, a lower value of 𝛼𝜙𝑜implies more 

penetration of CO2 and better recovery but lower value of 𝜙𝑜 also implies lower carbonation. It would seem to 

suggest that there could be an optimum pressure.  

Other than changing pressures, other features have been suggested (and we merely note these without review) 

are that the gas used be heated. Heat transfer by itself would be good in lowering the viscosity of the crude oil, 

and then mass transfer would lower the viscosity still further. Gas mixtures have also been proposed for use. 

However, if a gas mixture is made of one condensable and other non-condensable (such as methane) gases, 

then the non-condensable gas would form a “vapor blanket” and the condensable gas has to diffuse through it. 

It sets up a significant mass transfer resistance in the gas phase where for a single component there is none. 

Further, some more complex molecules such as dimethyl ether, are receiving attention at present, but nothing 

is known about its impact here. 

Finally, it is worth reminding that no report of field study has appeared yet although there are some 

indications that some start has been made. The field test by itself will show another set of problems and 

improvements.  



 11 

Conclusions 

VAPEX process uses no water to extract heavy oil, and that is its main advantage. More recent results than the 

seminal ones by Butler and Mokrys (1989) have shown yields to be significant in sandpacks and hence have 

been pursued vigorously. We have reviewed here the components needed for a simulation and find the perhaps 

the present call for success may be premature but show that the process could be improved.  

Conflicts of Interest 

The author(s) declare that they have no conflicting interests. 

Nomenclature 

𝐴𝜇  = constant 

𝐵𝜇  = constant 

𝐵𝑑 = solvent segment needed for a jump during diffusion 

𝑐 = free volume fraction of CO2  

𝐷 = diffusivity 

𝐷𝑜  = diffusivity at infinite dilution 

𝑓 = free volume fraction of oil 

𝐹𝑟  = Froude number  

𝑓𝑇   = total free volume fraction of oil 

𝑔  = acceleration due to gravity  

ℎ  = total pay zone height of the system  

𝑗𝐴
°   = diffusive flux  

𝑘  = permeability 

𝑃  = dynamic pressure  

𝑄𝑏  = recovery rate  

𝑅  = gas constant  

𝑅𝑒  = Reynold’s number  

𝑡 = time  

𝑇  = temperature  

𝑆𝐹 = swelling factor 

𝑆𝑐  = Schmidt’s number  

𝑣𝐴 = specific volume of the diffusing species A 

𝑊𝑏  = mass flux  

 

Greek Letters 

𝛼  = concentration dependence term  

Δ𝜌  = density difference 

Δ𝑆𝑜  = difference in fractional pore volume  

Δ𝑥𝑜 = changes in location of interface  

𝜇  = viscosity  

𝜇𝑜  = viscosity of pure oil 

𝜌 = total density  

𝜌𝐴 = mass concentration of the diffusing species A 

𝜌𝑜  = density of pure oil  

𝜌𝑠  = density of pure solvent dissolved in oil 



 12 

𝜙  = volume fraction of solvent 

𝜙𝑜  = solubility of solvent in oil 

𝜙𝑖  = initial concentration of solvent in oil     

𝛷 = porosity  

References 

Afsahi, B. and Kantzas, A. 2007. Advances in Diffusivity Measurement of Solvents in Oil Sands. J. Can. Pet. Technol. 

46(11): 56–61. PETSOC-07-11-05. 

Bachu, S. and Shaw, J. 2003. Evaluation of the CO2 Sequestration Capacity in Alberta’s Oil and Gas Reservoirs at 

Depletion and Effect of Underlying Aquifers. J. Canad. Pet. Tech. 42(9):51-61. PETSOC-03-09-02. 

Banerjee, D.K. 2012. Oil Sands, Heavy Oil and Bitumen, From Recovery to Refinery. PennWell Corp, Tulsa.  

Beckwith, R. 2012. The Tantalizing Promise of Shale Oil. J Pet. Tech. 64(1):30-36. SPE-0112-0030-JPT. 

Bird, R.B, Stewart, W.E, and Lightfoot, E.N. 2002. Transport Phenomena, 2nd ed. New York: John Wiley.  

Brinkman, H.C. 1947. On the Permeability of Media Consisting of Closely Packed Porous Particles. Appl. Sci. Res. 

A1:81-86. 

Buckingham, E. 1914. On Physically Similar Systems; Illustrations of the Use of Dimensional Equations. Phys. Rev 

4:345-376. 

Buckley, J. S., Hirasaki,G. J., Liu, Y., et al. 1998. Asphaltene Precipitation and Solvent Properties of Crude Oil. Pet. Sci. 

Tech.16: 251-285. 

Burke, N.E., Hobbs, R.E., and Kshou, S.F. 1990. Measuring and Modeling of Asphaltene Precipitation. J. Pet. Tech. 

42(11):1440-1520. SPE-18273-PA. 

Butler, R. M. and Mokrys, I. J. 1989. Solvent Analog Model of Steam-Assisted Gravity Drainage. AOSTRA J. Res. 

5(1):17-32. 

Chung, F. T. H., Jones, R. A., and Nguyen, H. T. 1988. Measurements and Correlations of the Physical Properties of 

CO2-Heavy Crude Oil Mixtures. SPE Reservoir Eng. 3(3): 822-828.SPE-15080-PA. 

Cohen, M. H. and Turnbull, D. 1959. Molecular Transport in Liquids and Glasses. J. Chem. Phys. 31(5):1164-1169. 

Cuthiell, D. and Edmunds, N. 2012.Thoughts on Simulating the VAPEX Process. J. Can. Pet. Technol. 52(3):192-203. 

SPE-158499-PA. 

Fadaei, H., Shaw, J. M., and Sinton, D. 2013. Bitumen-Toluene Mutual Diffusion Coefficients Using Microfluidics. 

Energy & Fuels. 27(4): 2042–2048. 

Flory, P.J. 1953. Principles of Polymer Chemistry. Ithaca: Cornell University Press. 

Fujita, H. 1961. Diffusion in Polymer-Diluent Systems. Adv. Polym. Sci. 3(1): 1-47. 

Ghanavati, J., Hassanzadeh, M., and Abedi, H. 2014. Critical Review of Mutual Diffusion Coefficient Measurements for 

Liquid Solvent+bitumen/Heavy Oil Mixtures. Can. J. Chem. Eng. 92(8):1455–1466. 

Gardner, J. W., Orr, F. M., and Patel, P. D. 1981. The Effect of Phase Behavior on CO2-Flood Displacement Efficiency. 

J. Pet. Tech. 33(11): 2067-2428. SPE-8367-PA. 

Hiss, T. G. and Cussler, E. L. 1973. Diffusion in High Viscosity Oils. AIChE J. 19(4): 698-703. 

Lake, L. W. 1989. Enhanced Oil Recovery. Englewood Cliffs, New Jersey: Prentice-Hall. 

Mohan, V., Neogi, P., and Bai, B. 2017a. Flory-Huggins Solution Theory for Heavy Oils. Canad. J. Chem. Eng. 

95(4):796-798. 

Mohan, V., Neogi, P., and Bai, B. 2017b. Concentration Dependent Diffusivities of Model Solvents in Heavy Oil. Diff. 

Fundam. 27(3):1-27. 

Mohan, V., Neogi, P., and Bai, B. 2019a. Revisiting Butler-Mokrys Model. Soc.Pet.Eng. J. 20(2): 511-521. SPE-194212-

PA. 

Mohan, V. and Neogi, P. 2019b. Extraction Rate in Vapor Assisted Extraction of Heavy Oil (VAPEX). Improved Oil 

and Gas Recovery 3:1-6. 

Miller, C. A. and Neogi, P. 2008. Interfacial Phenomena: Equilibrium and Dynamic Aspects. Cleveland: CRC Press, Inc. 

Neogi, P. 1988. Nonlinear Elastodiffusion at Small Deformations in Polymer Membranes. Chem. Eng. Comm. 68(1):185-

195. 

Nenninger, J.E. and Dunn, S.G. 2008. How Fast ID Solvent Based Gravity Drainage? Paper presented at Canadian 

International Petroleum Conference, Calgary, Alberta, June 2008. PETSOC-2008-139. 

Nghiem, L.X. and Coombe, D.A. 1997. Modeling Asphaltene Precipitation During Primary Depletion. Soc. Pet. Eng. J., 

2(2):170-176. SPE-36106-PA. 



 13 

Prausnitz, J.M., Lichtenthaler, R.N., and De Azevedo, E.G. 1999. Molecular Thermodynamics of Fluid Phase 

Equilibrium, 3rd edition. New Jersey: Prentice-Hall Inc. 

Rassenfoss, S. 2017. Shale Oil EOR Works, But Will It Make a Difference? J.Pet.Tech. 69(10):34-40.SPE-1017-0034-

JPT. 

Rojas, G.A. and Ali, S.M.F. 1988. Dynamics of Subcritical CO2/Brine Floods for Heavy-Oil Recovery. Soc.Pet.Eng. J. 

3(1):35-44. SPE-13598-PA. 

Shaw, J. and Bachu, S. 2002. Screening, Evaluation, and Ranking of Oil Reservoirs Suitable for CO2-Flood EOR and 

Carbon Dioxide Sequestration. J. Canad. Pet. Tech. 41(9):51-61. PETSOC-02-09-05. 

Slattery, J.C. 1972. Momentum, Energy and Mass Transfer in Continua,2nd edition. New York: McGraw-Hill Book Co. 

Speight, J.D. 2009. Enhanced Recovery Methods for Heavy Oil and Tar Sands. Houston: Gulf Publishing House. 

Subramanian, M., Deo, M.D., and Hanson, F.V. 1996. Compositional Analysis of Bitumen and Bitumen Derived 

Products. J. Chromatographic Sci. 34(1):20-26. 

Tran, T.Q.M.D., Neogi, P., and Bai, B. 2012. Free Volume Estimates of Thermodynamic and Transport Properties of 

Heavy Oils With CO2. Chem. Eng. Sci. 80(1):100-108. 

Tran, T.Q.M.D., Neogi, P., and Bai, B. 2015. Single Pore Model for Displacement of An Immiscible Heavy Oil in a 

Reservoir. Soc. Pet. Eng. J. 21(3):864-872. SPE-178425-PA. 

Tran, T.Q.M.D., Neogi, P., and Bai, B. 2017. Stability of CO2 Displacement of An Immiscible Heavy Oil Reservoir. Soc. 

Pet. Eng. J. 22(2):539-547. SPE-184407-PA. 

Upreti, S.R., Lohi, A., Kapadia, R.A., et al. 2007. Vapor Extraction of Heavy Oil and Bitumen: a Review. Energy & 

Fuels 21(3):1562-1574.   

Vrentas, J.S. and Duda, J.L.1979. Molecular Diffusion in Polymer Solutions. AIChE J. 25(1):1-24. 

Welker, J.R. and Dunlop, D.O. 1963. Physical Properties of Carbonated Oils. J. Pet. Technol. 15(8):873-873. SPE-567-

PA. 

Yang, C. and Gu, Y. 2006a. Diffusion Coefficients and Oil Swelling Factors of Carbon Dioxide, Methane, Ethane, 

Propane, and Their Mixtures in Heavy Oil. Fluid Phase Equilibria. 243(1):64-73. 

Yang, C. and Gu, Y. 2006b. A New Method for Measuring Solvent Diffusivity in Heavy Drop by Dynamic Pendant 

Drop Shape Analysis (DPDSA). Soc.Pet.Eng. J. 11(1): 48-57. SPE-84202-PA. 

Yang, D. and Gu, Y. 2008. Determination of Diffusion Coefficients and Interface Mass-Transfer Coefficients of the 

Crude Oil-CO2 Systems By Analysis of The Dynamic and Equilibrium Interfacial Tensions. Ind. Eng. Chem Res. 

47(15): 5447-5455. 

Zhang, X., Fulem, M., and Shaw, J.M. 2007. Liquid-Phase Mutual Diffusion Coefficients for Athabasca 

Bitumen+Pentane Mixtures. J. Chem. Eng. Data. 52(3):691–694. 

 
Vijitha Mohan received her PhD in Chemical Engineering from Missouri University of Science and 
Technology and served as a Lecturer after graduation. Her research interests lie in heavy oil recovery. She 
acquired her B.Tech in Chemical Engineering from University of Madras and M.S in Chemical Engineering 
from Mississippi State University. She is at present a Research Associate at University of Calgary, Canada. 
She has three prior publications in this area. 
 
Parthasakha Neogi, SPE, is a Professor of Chemical Engineering at Missouri University of Science and 
Technology, where he has worked as faculty for the last 36 years. His research interests are in wetting, 
surfactants and polymers, and in interfacial transport phenomena. He holds B.Tech. (Hons.) from the Indian 
Institute of Technology Kharagpur, M. Tech. from the Indian Institute of Technology Kanpur, and Ph.D. from 
Carnegie-Mellon University, all in chemical engineering. 
 
Baojun Bai, SPE, is the Lester R. Birbeck endowed Chair professor of petroleum engineering at Missouri 
University of Science and Technology. Previously, Bai was a reservoir engineer and head of the conformance 
control team at the Research Institute of Petroleum Exploration and Development, PetroChina. He also was 
a post-doctoral scholar at the California Institute of Technology and a graduate research assistant at the New 
Mexico Petroleum Recovery Research Center for EOR projects. Bai has over 20 years of experience in the 
area of EOR. He holds PhD degrees in petroleum engineering from New Mexico Institute of Mining and 
Technology and in petroleum geology from China University of Geoscience-Beijing. Bai has published more 
than 130 papers in peer-reviewed journals and international conferences. He served on the JPT Editorial 
Committee for the feature “EOR Performance and Modeling” during 2007-2013. He is a technical editor for 
SPE Journal and SPE Reservoir Evaluation and Engineering.