Experimental and numerical investigations into the effect of heterogeneities on the recovery of heavy oil by VAPour EXtraction (VAPEX)

We used a combination of experimental, analytical and numerical approaches to examine the oil drainage rate obtained from the VAPour EXtraction (VAPEX) process to recover heavy oil. In particular we investigated the influence of macroscopic scale heterogeneities through a series of experiments. These heterogeneities comprised layers, a quadrant model and two cases with discontinuous shale barriers above the injection well. All experiments were performed in well-characterised glass bead packs using glycerol and ethanol as analogues of heavy oil and solvent respectively. All the porous medium and fluid properties (including permeability, porosity, viscosity, density, diffusion and dispersion) were measured independently. The experimentally measured rates were compared to the estimates derived from the Butler– Mokrys analytical model. Numerical simulations were validated using the experimental observations from the homogenous and heterogeneous systems. The findings confirmed the square root dependency of the oil rate on permeability (at least for the range of permeabilities used here) which is consistent with findings from previous studies. Despite that, the results indicated that the Butler–Mokrys derived expression under-predicts the physical oil rate, even when the effects of convective dispersion and end-point density difference (as suggested by other works) were factored in. Results from the heterogeneous models suggest that layering may not reduce VAPEX oil drainage rates significantly. They also showed that oil was not recovered from the lower permeability layers. For models with discontinuous shale barriers, the simulations tended to over-predict the oil rate compared with the experiments. Overall the rates estimated from the simulations were more comparable to the physical rates than those estimated from the original Butler–Mokrys analytical derivation. Moreover the simulator was able to capture the pattern of solvent–oil distribution for both homogeneous models and heterogeneous models. This would suggest that the Butler–Mokrys model does not properly describe all the physical processes that are controlling the oil drainage rate even in homogeneous models. 2014 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http:// creativecommons.org/licenses/by/3.0/). 414 M.M. Al-Hadhrami et al. / Fuel 135 (2014) 413–426


a b s t r a c t
We used a combination of experimental, analytical and numerical approaches to examine the oil drainage rate obtained from the VAPour EXtraction (VAPEX) process to recover heavy oil. In particular we investigated the influence of macroscopic scale heterogeneities through a series of experiments. These heterogeneities comprised layers, a quadrant model and two cases with discontinuous shale barriers above the injection well. All experiments were performed in well-characterised glass bead packs using glycerol and ethanol as analogues of heavy oil and solvent respectively. All the porous medium and fluid properties (including permeability, porosity, viscosity, density, diffusion and dispersion) were measured independently. The experimentally measured rates were compared to the estimates derived from the Butler-Mokrys analytical model. Numerical simulations were validated using the experimental observations from the homogenous and heterogeneous systems.
The findings confirmed the square root dependency of the oil rate on permeability (at least for the range of permeabilities used here) which is consistent with findings from previous studies. Despite that, the results indicated that the Butler-Mokrys derived expression under-predicts the physical oil rate, even when the effects of convective dispersion and end-point density difference (as suggested by other works) were factored in. Results from the heterogeneous models suggest that layering may not reduce VAPEX oil drainage rates significantly. They also showed that oil was not recovered from the lower permeability layers. For models with discontinuous shale barriers, the simulations tended to over-predict the oil rate compared with the experiments.
Overall the rates estimated from the simulations were more comparable to the physical rates than those estimated from the original Butler-Mokrys analytical derivation. Moreover the simulator was able to capture the pattern of solvent-oil distribution for both homogeneous models and heterogeneous models. This would suggest that the Butler-Mokrys model does not properly describe all the physical processes that are controlling the oil drainage rate even in homogeneous models.

Introduction
The VAPour EXtraction (VAPEX) of heavy oil by solvents is considered to be one of the most energy efficient, economically attractive and pollution-free alternatives to thermal extraction methods . VAPEX is suitable for thin, shallow, low permeability reservoirs or those underlain by aquifers, since these can result in excessive heat loss and thus poor recovery efficiency for thermal Enhanced Oil Recovery (EOR), methods [14,16,29,30]. CO 2 based solvents have been shown experimentally to be more productive than methane and ethane based solvents [17,31]. VAPEX may therefore be a way of improving oil recovery whilst simultaneously storing excess CO 2 in the subsurface and thus reducing greenhouse gas emissions.
The mechanism was first proposed by Butler and Mokrys [9] as an analogue of Steam Assisted Gravity Drainage (SAGD) in that diffusion driven mass transfer between the solvent-heavy oil reduces oil viscosity in a similar way to the heat diffusion between the steam and oil. Two long horizontal wells are drilled parallel to each other (as in SAGD) in order to maximise the exposure to the reservoir. Solvent is then injected into the upper well while the diluted oil from the solvent-oil diffusion layer drains under gravity to the lower well (see Fig. 1).
The main uncertainty limiting the field application of this process is the oil drainage rate and the effect of reservoir heterogeneity on the rate. Oil drainage rates are much lower than in SAGD because the oil viscosity is reduced via mass diffusion and dispersion rather thermal diffusion and it is unclear whether these low rates will be economic. It is very difficult to perform accurate numerical simulations of the process due to the levels of grid refinement needed to resolve the diffusion and drainage processes occurring at the solvent-oil interface. To circumvent this issue Butler and Mokrys [9] derived a semi-analytical equation to estimate oil drainage rate without the need for simulation. Their analysis assumed that the mixing between the solvent and the heavy oil driven by molecular diffusion acts to reduce the viscosity of the oil in a similar way to the way that thermal diffusion reduces the viscosity of oil during SAGD. The oil drainage rate during VAPEX can then be estimated from: where k is permeability, g is the acceleration due to gravity, h is the reservoir thickness, / is the porosity and DS o is the initial mobile oil saturation. The constant N s in Eq. (1) is determined from the fluid properties as: where Dq is the density difference between the diluted oil-solvent mixture and vapour, D m is the molecular diffusion coefficient between the vapour and the oil, C s is the concentration of the solvent, lðc s Þ is the viscosity of the solvent-oil mixture, c max is the maximum solvent concentration at the outermost boundary of the diffusion layer, and c min is the minimum solvent concentration at the innermost boundary of the diffusion layer. This analysis assumes first contact miscibility between the oil and the vapour. This is clearly much simpler than performing a fine grid simulation to predict VAPEX performance but unfortunately seems to systematically underpredict the drainage rate seen in laboratory experiments [1][2][3][4][5][6][7][8]13,14,[16][17][18][20][21][22][23][24][25][26][27][28] In contrast, the only published field pilot (in Soda Lake, Saskatchewan) obtained oil production rates a factor of 10 lower than predicted. A further uncertainty in the field scale application of VAPEX is the impact of geological heterogeneity. Small (mm to cm) scale heterogeneities may be expected to increase dispersive mixing and thus the oil drainage rate. Larger scale heterogeneities (greater than 10 m) may alter the shape of the solvent cone and thus either improve or, more likely, reduce the oil drainage rate, depending on the nature of the heterogeneity.
Worldwide, most extra-heavy oil reservoirs exhibit considerable geological heterogeneity at different scales (bearing in mind the additional fluid heterogeneities which are beyond the scope of this study). In commercial SAGD projects, reservoir geological features and associated petro-physical properties are the main influences on successful heavy oil extraction [32]. For instance, the McMurray Formation, within the Athabasca oil sand deposit in Canada, is considered to be one of the most complicated extraheavy oil (more than 10 6 m Pa s) and low gravity (4-8°API) geological deposits. It is made up of a continental sequence of unconsolidated sands and shales overlaying an unconformity surface of Devonian limestone. The Lower McMurray zone is a fluvial dominated succession consisting of highly braided channels and sand bars of clean blocky sands associated with intraclast mudstone breccias whereas the Upper McMurray zone is a highly channelized mixed fluvial-tidal estuarine deposit. The scale and style of these channels varies wildly so that lateral continuity between wells is almost impossible (even for short distances of <50 m). The channel sands in this unit grade upwards into inclined heterolithic stratified units of rippled sandstones and are heavily burrowed with silty mudstone and clays. This zone is very fine grained and argillaceous with, often, very low vertical permeability as a result of the widely dispersed shale lenses formed from clays and low permeability siltstones [33,34].
There is a significant literature in SAGD investigating the impact of discontinuous shales on oil recovery [35][36][37][38], using either physical or numerical heterogonous models. The shales disrupt the ideal V shape of the steam chamber but do not always change the oil rates and/or recovery. Chen et al. showed that the growth of the steam chamber and the resulting oil drainage were highly sensitive to shales between the injection and production wells [35]. The growth of the steam chamber was adversely affected by discontinuous shales above the injection well but this did not seem to affect the oil drainage rate. The expansion of the steam chamber to the top of the reservoir, however, was reported to be adversely affected in cases where the shale barriers in the above well region were either continuous or small but with high density. Similarly Yang and Butler's SAGD experiments suggested that models with a high permeability layer at the top of the reservoir provided better vertical communication and hence more production [38]. They also observed that heat was still able to penetrate by conduction through the shale barriers, although the heated oil was not drained to the producer and was trapped due to the pressure build-up of undiluted oil underneath the shale.
It is likely that VAPEX performance will respond differently to heterogeneity from SAGD. This is because the rate of mixing between the solvent and the oil will be reduced by the tortuosity of the pore space of the rock (which is affected by pore connectivity, the porosity and connate water saturation) whereas heat diffusivity depends only upon the conductivity and specific heat capacity of the rock (which is controlled by the mineralogy of the rock as well as the porosity and connate water saturation). This is exemplified by shale barriers. There will be minimal molecular diffusion through a shale barrier because of the low porosity, low permeability and the fact that the pore space of many shales is filled with water. In contrast thin shales may allow heat to be con-ducted through them because water has a relatively high conductivity but thicker shales will result in a higher heat loss because of the thermal heat capacity of the shale and the water contained within the shale. The influence of permeability variations on the VAPEX solvent chamber, however, is likely to be very similar to their influence on SAGD.
A limited amount of research has been performed to investigate experimentally the effect of vertical and lateral variations in permeability and of barriers on VAPEX mechanisms. An experimental study was carried by Jiang and Butler [39] in a 2D (35.6 Â 21.6 Â 3.2) cm pack, where physical systems with thin low permeability continuous layers, low permeability discontinuous layers and high permeability sand lenses, were investigated. In this study the low permeability sand was nearly 43.5D while the high permeability sand was 217D. As would be expected the homogeneous low permeability packs produced at a lower oil drainage rate than systems with high permeabilities. It was also observed that continuous shale/low permeability barriers prevented the formation of the solvent-vapour chamber above these barriers (as might be expected) and therefore reduced the VAPEX performance significantly. Different patterns of vertical and horizontal fractures were also investigated in Jiang and Butler study. Similar to SAGD, it was observed that vertical fractures tended to improve the recovery, especially in layered systems.
In a different study performed using numerical simulations to evaluate the impact of reservoir permeability distributions on VAPEX, Zeng et al. concluded that a random permeability distribution had a very minor effect on overall VAPEX performance compared to homogeneous cases [40]. It was also noted that the highest oil drainage rates were obtained for models with a high permeability layer close to the producer.
Frauenfeld et al. meanwhile, tested the impact of sand lenses and layering, as well as bottom aquifers [29]. Their physical model was constructed from a field scale reservoir model applying Pujol and Boberg's gravity and diffusion scaling criteria [41]. The model simulated a 30 m thick reservoir underlain by a 10 m thick aquifer with a geometric ratio of 100:1 (i.e. 30 cm reservoir thickness, 10 cm bottom water zone and 25 cm horizontal well offset). In these experiments, Kerrobert oil (50,000 mPa s) was used with butane as the solvent. It was observed that oil drainage rates were lower in layered systems compared to homogeneous models with uniform permeabilities, and that the layering resulted in the formation of mini-vapour chambers above the injectors. Interestingly, it was noticed that low permeability lenses did not severely inhibit the oil rate as the solvent was diverted sideways and around these features.
In this paper we investigate the influence of macroscopic or reservoir scale heterogeneities on VAPEX processes in order to understand the likely influence of heterogeneity on VAPEX performance  and to test the ability of commercial simulations to predict these outcomes. We achieve this by using a combination of well characterised laboratory experiments and numerical simulation. We examine the impact of idealised permeability heterogeneities representative of those encountered in real fields on both oil drainage rate and cumulative oil recovery and compare these results with those obtained from homogeneous models. The fluid and porous media properties were independently characterised in all the laboratory experiments so that the results could be compared directly with the predictions of a commercial reservoir simulator, without the need for history matching.

Bead pack design
The experiments were all performed in a 2D bead pack constructed from clear glass Perspex, with internal dimensions 40 cm Â 20 cm Â 0.6 cm. The length (40 cm) and height (20 cm) of the pack were chosen to be large enough to capture the influence of macroscopic heterogeneity (keeping the aspect ratio at 2:1) but small enough to ensure uniform packing. The thickness of the pack (0.6 cm) was chosen to be thick enough to allow multiple beads across the dimension but thin enough to ensure that the fluids were uniformly distributed across this thickness by diffusion.
The material of the pack was selected to permit visualisation of the fluid flow and to be strong enough to withstand experimental conditions (20°C, 101.32 kPa). The pack represented a symmetry element of a vertical cross-section through a reservoir undergoing VAPEX, orientated perpendicular to the horizontal injection and production wells. The producer was located at the lower corner of the model while the injector was directly above it, halfway between the producer and the top of the model reservoir (see Fig. 2).
The top cover of the model was designed to be removable so as to allow the beads to be packed and removed. The internal surfaces and edges of the pack were coated with a thin layer of soft silicone sealant material to mitigate as far as possible edge effect problems and any lack of consolidation of the packing near the boundaries.
A number of different permeability patterns were created within this pack using different sized beads (Grades 6 and 11). These are shown in (Fig. 3). Three of these patterns were formed of layers of contrasting permeabilities. The aim was to investigate the impact of layer ordering on oil rate and recovery. The fourth permeability pattern investigated was that of a quadrant. This was chosen as it is both geologically relevant (e.g. when a fault causes an offset between layers) and particularly challenging to model flow in this pattern using conventional reservoir simulators [42]. The low permeability layer was formed from the very small beads (Grade 11) while the high permeability layer was formed from coarser glass beads (Grade 6). This resulted in a 1:15 permeability ratio between the layers.
Two packs were also constructed to investigate the impact of the length of a discontinuous shale just above the injection well on rate and recovery (Fig. 4). The shale lengths and location were  selected based on some preliminary numerical simulations. According to the simulations the chosen shale location would interrupt the segregation of the V-chamber to the top of the reservoir and would have most effect on oil rate and recovery. The first model consisted of a homogenous high permeability pack with a single discontinuous shale barrier that was 25% of the total reservoir length ($10 cm), while in the second model the length of the barrier was 50% of the reservoir length ($20 cm).
In order to generate these various patterns, very thin sections of glass Perspex with a thickness equal to the internal thickness of the model were positioned gently into the pack (avoiding scratching the surface). Beads were then gradually loaded into the vertically clamped pack using the method described by Caruana and Dawe [43]. We ensured the beads were closely packed by tapping the pack gently with a plastic hammer at intervals during packing. Each layer was packed individually while sealing the other layers from the open side so as to avoid any mixing between the beads/ layers. The Perspex sections were then slowly removed again while tapping the pack to re-distribute the beads and fill all the empty spaces. A compressible rubber strip was used to model the shale. In each case the strip was slightly wider than the model's internal thickness to ensure a complete sealing at the edges and to avoid movement either of beads or of the barrier during the experiment.

Analogue fluids
Following Alkindi et al. [5,6] we used glycerol to represent the oil and ethanol to represent the vapour. These analogue fluids were chosen so that: (a) Their properties were representative of the real fluids used in previous experimental studies [11]. In particular we needed to ensure that the viscosity ratio was such that the model oil was essentially immobile compared with the solvent and that the density difference between the model oil and model solvent were typical of real fluids. (b) The fluids were first contact miscible as assumed by Butler and Mokrys [9] in deriving their analytical expression for oil drainage rate (Eq. (1)). (c) We were able to fully characterise the fluid properties for input into the simulator, to avoid the need for history matching. The fluid properties are given in Table 1. The effect of ethanol concentration on glycerol viscosity and density has been determined previously by Alkindi et al. [3] and can be determined using the empirical correlations: As shown in Fig. 5 below, the experimentally measured oil-solvent mixture viscosities as a function of solvent concentration were well represented by the ideal mixing rule used in the simulator.
The concentration of ethanol and glycerol in the produced fluids from the experiments was determined from the refractive index of these fluids and a calibration curve of concentration of ethanol as a function of refractive index. The calibration curve was determined by measuring the refractive indices of 32 different known ethanol-glycerol concentrations using a digital refractometer (Stanley Abbe Refractometer, model 60/ED). This has an accuracy of ±0.0002 accuracy. The prism of the refractometer was carefully rinsed with either ethanol or acetone. The refractometer was first calibrated using pure samples of water, ethanol and glycerol, and the measured values at room conditions compared well with the values in literature [46]. During the sample preparation, the volumes of the ethanol and glycerol were gently controlled using a one millilitre syringe. To eliminate experimental errors the measurements for each concentration were repeated at least three times.

Petro-physical properties
The values of the porosities, permeabilities and longitudinal dispersion associated with each grade/size of beads were measured using a smaller bead pack that was 17.3 Â 0.6 Â 9 cm. This mini pack had an inlet and an outlet along opposite sides (to represent a line drive).
The porosity, permeability and longitudinal dispersion associated with each bead size are given in Table 2. The porosity of each size of bead was measured by weighing the pack before and after fully saturating it with a liquid and dividing the pore volume by the total internal bulk volume (assuming that the inlet and outlet lost fluid volumes was negligible). The absolute permeability of the different beads, meanwhile, was simply calculated from Darcy' Law by measuring the pressure difference across the inlet and outlet of the pack using pressure transducers while varying the injection rates. The longitudinal dispersion coefficient for each bead size was determined using the method of [47] from effluent profiles resulting from vertical upwards, glycerol-ethanol displacements in the mini-pack [48]. These displacements were performed at the same injection rates as used in the VAPEX experiments. This flow configuration with glycerol displacing ethanol was chosen in order to eliminate viscous fingering (see discussion in [3,49]). Transverse dispersion was not measured as previous work by Alkindi et al. [5] has shown that its effect in these experiments is negligible.

Experimental investigations of VAPEX
Before starting an experiment, CO 2 was injected at a very low rate into one side of the pack whilst leaving the outlet on the opposite side of the pack open (see Fig. 4), in order to remove any air within the pack. The pack was then saturated with glycerol. This was achieved by mounting the pack vertically an injection port at the bottom and an outlet port at the top (i.e. in a line drive configuration) and injecting glycerol very slowly through the inlet whilst allowing the displaced CO 2 to escape from the outlet port. Once the pack was saturated with glycerol (note than any remaining CO 2 would have dissolved in glycerol) then this outlet port was sealed and the outlet port underneath the inlet port was opened (to represent the cross-section through the injection and production well). Finally the pack was orientated so the inlet and outlet port were both on one side.
During all the analogue VAPEX experiments ethanol was injected at a constant rate using a ''Series 12 Â 6 HPLC'' pump (flow rate precision of ±2%), while the producer was held at atmospheric pressure. Simulations have previously shown that this results in a constant pressure drop between the injector and producer as is required for a VAPEX process [6]. The ethanol was dyed with Lissamine green (C 27 H 25 N 2 NaO 7 S 2 ) in very small concentrations (0.03 g/400 cc) in order to be able to observe the solvent chamber segregation. This allowed images to be taken of the fluid distributions and then analysed. This dye was chosen because it does not affect the measured refractive index and, unlike methylene blue, does not bleach over time once in contact with glycerol [6]. The effluent produced was collected using 1 ml measuring cylinders for accuracy in measuring the production volumes. The solvent injection rates were chosen so that all flows within the packs were gravity dominated (one of the requirements of the VAPEX process).
This was achieved by calculating the viscous to gravity number for horizontal flow using Fayers and Muggeridge [50] expression: , where Q is the injection rate, H is the height of the pack, w is the internal thickness, and / is the effective porosity of the pack.
Based on these calculations, and also bearing in mind the minimum rates that could be achieved by our pump, the solvent injection rate was set to 0.11 cm 3 /min in the homogeneous experiments and 0.07 cm 3 /min in the heterogeneous experiments. These result in gravity numbers (see Table 3) that are less than 10 but still greater than 0.1, except for the very lowest permeability pack. Based on Eq. (5) these results suggest that flows in our experiments will be influenced both by gravity and viscous effects. Preliminary screening experiments with these rates showed that a V shaped solvent-oil interface (as is expected in VAPEX gravity dominated flows) formed even in the lowest permeability experiment. These experiments suggest all our experiments were gravity dominated despite the fact that our calculations give N g=v < 10. We note that this number was derived for line drives in horizontal or slightly tilted beds and thus may not be applicable to VAPEX flows.

Numerical simulation
The experiments were simulated using the semi-compositional numerical simulator CMG-STARS (Computer Modelling Group, Canada). A 2D (x-z) Cartesian grid was used with the same dimensions as the experimental physical model. We used a 200 Â 160 Â 1 grid in x, z and y directions, respectively. This was selected after a grid refinement study to ensure that physical dispersion dominated the numerical dispersion and that the gravity dominated flow, particularly at the leading edge of the vapour chamber.
The porous medium and fluid properties were used directly as input data so that there was no history matching involved. The only exception to this was that we increased the permeability of the top five rows of grid blocks by a factor of 1000 in order to properly capture the advance of the leading edge of the vapour chamber. The oil and solvent were modelled as two components and no water or gas phases were included in the simulation. The ideal mixing rule was used to calculate the viscosity of solvent-oil mixtures. The vapour-oil relative permeabilities were entered as  straight lines with zero capillary pressures as ethanol and glycerol are first contact miscible. The reference pressure was chosen to be atmospheric (101 kPa) and the reservoir temperature was set to be room temperature (20°C). STARS does not have velocity dependent dispersion at the time of writing so the physical dispersion of ethanol in glycerol was modelled by: (a) finding the maximum solvent-oil frontal advance rate, (b) calculating the corresponding longitudinal dispersion and then, (c) using that value for the diffusivity in the simulation. The injector and producer were operated under constant rate and constant bottom-hole-pressure (BHP) respectively as in the experiments. The injection well was located in the centre of one face of the model and the production well was located in the corner of the model directly underneath the injection well in order to replicate the locations of inlet and outlet in the experiment.

Results and discussion
We first discuss the results from our homogeneous experiments as we can only really evaluate the impact of heterogeneity on VAPEX if we compare the results with those obtained from the homogeneous models.

Homogenous systems
The physical cumulative oil production, oil drainage rates and solvent oil ratios (SORs) as a function of time for the three different homogeneous packs are shown in Fig. 6. Based on these data, the average stabilised oil drainage rates were estimated to be 1.20 ± 0.51, 2.90 ± 0.47 and 4.2 ± 0.34 cm 3 /h, for the models with grade 11, 9 and 6 beads, respectively. As expected the highest drainage rate, highest recovery and lowest SOR is obtained from the highest permeability pack.
The measured oil drainage rates were compared with the analytically predicted rates estimated from the Butler and Mokrys [9] equation (Eq. (1)). It is well known that this equation correctly predicts the oil rates obtained from experiments performed in Hele-Shaw cells but underpredicts the oil drainage rate obtained from experiments in porous media [1-8,12,13,15-17,20-23, 25-28,50,52]. Various modifications have been proposed to the Butler and Mokrys [9] to improve its predictions including replacing the   (2)). Table 4 compares the oil drainage rates calculated from Eq. (1) using the molecular diffusion, the longitudinal dispersion and the longitudinal dispersion + endpoint density difference in the calculation of N s . The table also gives the values of N s used (C min was 0.01 and C max was 1). These values of N s are comparable with those obtained in a range of other studies (e.g. [8,9]).
It can be seen that all the experimental oil drainage rates are systematically higher than those obtained from the Butler-Mokrys [9]equation even when using longitudinal dispersion and the endpoint density difference. This difference is similar to that seen by Das and Butler as can be seen from fig:6 where we have plotted the observed and predicted oil drainage rates as a function of the square root of permeability [16]. This discrepancy has also been observed by other workers [1][2][3][4][5][6][7][8]12,13,[15][16][17][20][21][22][23][25][26][27][28]50,52] and suggests that there is some problem with the Butler-Mokrys equation.
Various explanations for this underprediction of oil rate have been proposed including higher than expected mixing. This has lead many authors to suggest that the mixing between solvent and oil is not driven purely by molecular diffusion. Various processes by which mixing may be increased have been proposed including convective dispersion and increased surface area between the fluids due to the formation of oil films on the pore walls caused by capillary imbibition. In contrast other authors have suggested that the dependency of oil rate on reservoir thickness is greater than the square root proposed by Butler and Mokrys [9]. Alternatively Zainee et al. [53] suggested that the gravity drainage rate was higher than expected because it depends upon the end point density difference beetle the solvent and oil rather than the concentration dependent difference suggested by Butler and Mokrys. Our results show that including convective dispersion and the endpoint density difference improve our estimates of drainage rate but they are still a factor of $4 too small. Capillary pressure cannot be an explanation in our experiments as our fluids are first contact miscible. Fig. 7 does at least confirm that oil drainage rate does depend on the square root of permeability as predicted by the Butler-Mokrys equation.
The experimental results were also compared to the predictions obtained using the semi-compositional numerical simulator from Computer Modelling Group (CMG) STARS. Fig. 8 compares the solvent-oil distributions predicted by the simulator (second and fourth rows) with the observed saturation profiles (first and third rows) for the highest permeability and lowest permeability homogenous packs. It can be seen that overall the simulations captured the pattern of the solvent segregation and the physical diffu- Fig. 8. Comparison of solvent-oil distributions for both physical and numerical models of homogeneous packs (a) highest permeability model (formed from grade 6) and (b) lowest permeability (formed from grade 11). However the solvent injection rates were different in these experiments. Fig. 7. The asymptotic underestimation of oil rates observed in both Das and Butler [14] experiments and this work. sion boundaries observed in the experiments. Looking more carefully at these figures it can be seen that the gravity override of the solvent chamber in the experiments is more advanced than that observed in the simulations (see Fig. 9).  Table 5. It can be seen that the agreement between the numerical simulations and the experiments is much better than the agreement between the Butler-Mokrys estimates and the experiments but that the simulation still under-predicted the oil rates by at least $22 %. This contrasts with the findings of Alkindi et al. [5,6] who found an excellent match between the experimental and simulation results when using the same fluid system but in a smaller pack and applying lower injection rates.
It is possible that the more advanced gravity over-ride in the experiments (compared with the simulations) could be due to the edge effect, despite our use of silicone sealant. Alternatively this could be due to the way in which gravity is treated in all the simulator. STARS, in common with most commercial simulators, uses upstream weighting to determine fluxes into and out of grid cells. Standard upstream weighting can lead to errors when flow is gravity dominated as stagnant points or sonic points can occur [54]. Similar problems in the modelling of gravity dominated flows in vertical miscible displacements through a system with a single shale have been observed by [55]. Methods to improve the modelling of gravity dominated have been discussed by [56].

Impact of layering
Figs. 10 and 11 compare the measured oil drainage rates and cumulative recoveries as a function of time for the homogenous and layered systems. The highest oil rates and recovery were observed in the homogenous model formed from the highest permeability beads (Grade 6). The lowest oil rates and recovery were observed in the homogenous model formed of the lowest permeability beads (Grade 11). The oil rates for the homogeneous model with intermediate permeability (Grade 9) were between these two. The interesting observation here is that the oil rate and recoveries for the various layered models were also just between the low and high case and exceeded the recoveries obtained from the intermediate permeability pack for some systems. We also note that a higher oil drainage rate and oil recovery are obtained when the higher permeability layer is above the lower permeability layer.  Table 5 Average oil drainage rates for high and low permeability packs obtained experimentally and predicted by simulation. The simulator underpredicts oil rate but not as badly as the Butler-Mokrys model.

Model description
High permeability case (injection rate of 1.     12 compares the solvent-oil distributions seen in the experiments with those predicted at the same times by simulation. It can be seen that in both simulation and experiments the solvent chamber grows preferentially in the higher permeability layer, regard-less of the position of that high permeability layer. This is particularly obvious in the quadrant model where the juxtaposition of a lower permeability quadrant with the higher permeability quadrant prevents spreading of the vapour chamber that has been  growing in the higher permeability quadrant. This is despite the fact that the lower permeability layers has a permeability of 10D. Clearly the flow is controlled by the permeability contrast between the layers rather than the absolute permeability. The cumulative oil recovered from these layered systems is still higher than that obtained from the lowest permeability homogeneous model (Fig. 11) despite this bypassing of oil in the lower permeability layers.
The simulator is less successful at predicting the experimental flows in layered systems than in the homogeneous systems. Once again the simulated flows appear to be less gravity dominated than in the experiments, although the difference appears to be worse in these layered systems. Table 6 compares the oil drainage rates obtained experimentally with those predicted by the simulator. It is possible this is due to there being a higher permeability streak along the top of the bead pack due to inconsistent packing, however, as noted in the method section, the permeability of the top few rows of grid blocks in the simulation has been increased to replicate this. Again the simulator underpredicts oil drainage rate but the difference is much larger than in the homogeneous systems. Moreover this also occurs in the case of the higher permeability layer underneath the lower permeability layer which suggests that it is not caused by the edge effect.
One challenge in the modelling of VAPEX in heterogeneous reservoirs is how best to upscale the absolute permeability. Jiang and Butler [39] proposed that this should be the arithmetic mean permeability. This would be consistent with the recovery being driven by the horizontal expansion of the vapour chamber. Alternatively it could be argued that the most appropriate upscaled permeability should be something between the harmonic mean permeability (appropriate for vertical flow) and the arithmetic mean (horizontal flow) as the oil is draining diagonally downwards along the interface between the oil and solvent chamber. Table 7 lists the arithmetic and harmonic mean permeabilities for each of the models. For the quadrant model we also calculated the geometric mean permeability (38.7D) and the effective permeability using renormalization (23.2D, see equation in [57]). In Fig. 13 we plot the observed oil drainage rate against the square root of the arithmetic and the harmonic average permeability for each of the heterogeneous packs. It can be seen that the oil drainage rates seen in the layered models fall very close to the straight line formed by plotting the oil drainage rate versus square root of permeability for the homogeneous models when the arithmetic mean permeability is used. The rates do not fall on this line when the harmonic mean permeability is used. This indicates that, as suggested by Jiang and Butler [39], the arithmetic mean is a more appropriate way to upscale absolute permeability in layered systems. We do note, however, that lower oil drainage rates were obtained when the lower permeability layer was above the higher permeability layer than when the higher permeability layer was on top but both cases have the same arithmetic mean permeability. We also note that the solvent does not reach the top of the model when the low permeability layer is on top. In this case (which also applies to the three layer and quadrant models) the effective permeability must be calculated for the region of the reservoir accessed by the solvent and the drainage height is also reduced. According to Eq. (1) the oil Table 7 The arithmetic and harmonic average calculated permeabilities for the layered and quadrant bead packs.

Model description
k Arithmetic k Harmonic  13. Stabilised measured oil drainage rates for the layered and quadrant models plotted as a function of the square root of the effective permeabilities using arithmetic and harmonic averaging methods and also using a scaling based on drainage height. As can be seen, the correlation of the oil rate with the square root of permeability times scaled drainage height is most consistent with the results from the homogeneous packs.  Comparison of oil drainage rates versus time for the high permeability case, high permeability with 10 cm discontinuous shale and 20 cm discontinuous shale. The standard deviation in the measured oil rates were (r ± 0.32, 0.21 and 0.16 cm 3 / h) for the high permeability without shale, high permeability with 20 cm shale and high permeability with 10 cm shale models respectively. rate would be reduced by where h is the thickness of the high permeability layer and H is the overall reservoir thickness. These results are also plotted Fig. 13 and are most consistent with the results from the homogeneous packs.

Impact of discontinuous shale barriers
Figs. 14 and 15 compares the cumulative oil recovery and oil drainage rates obtained from the bead packs containing a single discontinuous shale with those obtained from the homogeneous pack using the same Grade 6 beads. The averaged stabilised drainage rates for the three systems are given in Table 8. Again this could be due to an edge effect along the shale bottom, as the beads will be less densely packed here, but it may also be due to the way gravity is modelled in the simulator, as discussed earlier.
As expected, the oil rates for models with the shale were slightly lower than those seen in homogenous packs. The interesting result here is the slight improvement in oil drainage rates as a result of increasing the length of the shale. It should be noted, however, that this difference is close to the uncertainty in the measurements, considering the error margin of ±0.21 and ±0.16 cm 3 /h respectively for the two scenarios. It is also interesting to note that the effect of the shale on cumulative oil recovery is much less occurs when there are layers in the system. This is because the only that is bypassed is just above the shale whereas in layered systems complete layers remain unswept, especially if the top layer is a lower permeability than the layer associated with the injection well.
Figs. 16 and 17 shows the solvent-oil distributions seen in the experiments with those predicted by simulations while Table 9   Table 8 The experimentally measured average drainage rates for the cases with a shale barrier.

Model description
Average experimental oil rate (cm 3    compares the average oil drainage rates seen experimentally and obtained from the simulations. The agreement between experiments and simulations is reasonable in both cases although once again the simulations seem to predict less gravity dominated flows than the experiments. It can also be seen that the simulation predicts very little change in oil rate when there is a shale in the system whilst there is a significant drop in oil drainage rate seen in the experiments. The numerical simulation models were able to capture the solvent-oil distributions and the basic flow behaviour seen in the experiments although, as noted in previous sections, the level of gravity segregation was higher in the experiments than in the simulations. Again, the experiments gave a higher oil drainage rate than was predicted by the simulations (see Table 9). This could be due to an ''edge effect'' along the shale bottom, as the beads will be less densely packed here, but it may also be due to the way gravity is modelled in the simulator, as discussed earlier.

Summary and conclusions
We have investigated the impact of simple macroscopic heterogeneities on the performance of VAPEX using a mixture of experiments and very fine grid numerical simulations. The fine grid was needed to capture the development of the vapour chamber and the ensuing drainage of the diluted oil under gravity as well as to ensure that physical diffusion dominated the solvent oil mixing rather than numerical diffusion. The experiments used wellcharacterised glass bead packs and the analogue fluid system of glycerol and ethanol, previously used by [2][3][4][5][6]. These experimentally measured data were assigned in all the simulations, so that the predictions could be compared directly with the experimental results without the need for history matching.
The laboratory measurements of oil rate obtained from homogeneous bead packs were consistently higher than estimates obtained from the Butler-Mokrys analytical expression, although the drainage rate did vary linearly with the square root of permeability as predicted by this expression. Incorporating the effect of convective dispersion (as suggested by earlier works) and using the end-point density difference instead of the concentrationdependent density difference slightly improved the estimated rates. The stabilised oil rates predicted by simulations were closer to those obtained experimentally but still were slightly lower than the measured oil rates in both homogeneous and heterogeneous systems. It appeared that the discrepancy between the predicated and the physically measured rates was more marked for the model with higher permeability (in which flow was more gravity dominated) whilst a better agreement was obtained for low permeability models.
Comparisons between the observed and predicted solvent-oil distributions suggested that flow was less gravity dominated in the simulations than in the experiments. It is possible that the horizontal growth of the vapour chamber in the experiments was due to higher permeabilities at the top of the pack however increasing the permeability in the top few rows of grid blocks in the simulation did not improve the match. It seems possible that this discrepancy is due to the fact that simple upstream weighting does not capture flows properly around the sonic point. Improved numerical methods such as those described by [56] may be needed to improve predictions of gravity dominated flow.
These results indicate that VAPEX is most likely to be adversely affected by layering in reservoirs. This effect of this layering can be approximated by using an arithmetic mean of the layer permeabilities although it does not capture the impact of layer ordering on oil drainage rate and recovery. Upscaling is improved if the arithmetic mean permeability (or absolute permeability of the upper layer does not contribute to the oil recovery) is multiplied by the effective drainage height. It would appear that permeability contrast is most important in determining whether a layer will be swept rather than its absolute permeability. A single discontinuous shale just above the injection well will have a slight impact on oil drainage rate and oil recovery.
Further work needs to be done to determine the exact cause of the discrepancies between the simulation and the experiments. Until this is resolved the quantitative results from even fine grid simulations of the VAPEX process should be treated with caution although the qualitative differences in recovery between different realizations of geological heterogeneity are probably correct. The impact of smaller scale heterogeneities on VAPEX also needs to be investigated, especially those that reduce the vertical permeability of reservoirs as we expect that these will significantly affect the oil drainage rate.