Case Study: Modeling the Phase Behavior of Solvent Diluted Bitumen

In Western Canada, thermal recovery methods such as cyclic steam stimulation and steam assisted gravity drainage are the methods of choice to recover heavy oil and bitumen with viscosities exceeding 10,000 mPa.s. These methods require significant volumes of natural gas and water to generate steam: approximately 34 m³ of natural gas and 0.2 m³ of groundwater (assuming 90 to 95% recycle) to produce one barrel of bitumen^(1). Solvent based and solvent assisted recovery methods are a potential alternative to reduce or replace steam usage. However, potential solvents, such as light n-alkanes⁽²⁾, are expensive relative to heavy oil and the success of process depends on how much solvent can be recovered. Predicting the performance of solvent-based and solvent-assisted processes (including both oil and solvent recovery) is challenging because the introduction of a solvent can lead to complex phase behavior. For any given heavy oil and solvent mixture, it may be necessary to predict the phase boundaries, amounts and compositions for liquid-liquid (LL), vapour-liquid (VL), vapour-liquid-liquid (VLL), and asphaltene precipitation regions.

Cubic equations-of-state (CEoS) are widely used for modeling the phase behavior of crude oil^{(3, 4)}. To use CEoS, the composition, critical properties, molecular weight and density of each component or pseudo-component of the fluid are required. However, a complete compositional analysis of heavy oils and bitumens is not available. Even with high temperature vacuum distillation, only 20 to 40 wt% of the heavy oil is distillable. Simulated distillation assays can extend the characterization to approximately 50 wt% of the heavy oil but the boiling point calculations are based on extrapolations that cannot be validated against measured data. Hence, for heavy oils, a method is required to extrapolate the true boiling point curve to obtain accurate phase behavior predictions.

Castellanos Díaz et al^{(13, 16)} developed such a method to characterize Athabasca bitumen and predict the phase behaviour of mixtures of solvent and bitumen using the Advanced Peng Robinson Equation of state⁽²⁸⁾. They assessed several extrapolations of the true boiling curve including Guassian and Gamma distribution based extrapolations. They then divided the extrapolated curve into number of pseudo-components and assigned the critical properties and accentric factors for the pseudo-components using well known property correlations. Interaction parameters between solvent and bitumen pseudo-components were optimized to fit the experimental saturation pressures of pseudo-binary mixtures of bitumen and propane or carbon dioxide. The model correctly predicted saturation pressures and the boundaries of the vapor-liquid-liquid region for pseudo-ternaries of these components. However, asphaltene precipitation data were not available for these mixtures so model predictions for asphaltene precipitation could not be tested.

While there are considerable data for asphaltene precipitation from solvent diluted dead heavy oils^(17–21) and some data for LL and VL equilibrium of heavy oils and light solvents^(5–13), there are few datasets that include both VLE and asphaltene precipitation data, particularly for live oils^(22–24). This combination of data is required to ensure that the fluid characterization and model are self-consistent.

The following case study is intended to test the characterization method proposed by Castellanos-Diaz et al^{(13, 16)} on dead and live bitumen mixed with n-pentane and a multi-component condensate. For n-pentane, the composition at the onset of asphaltene precipitation and the amount of asphaltene precipitated from the dead oil at different solvent contents were measured at room temperature condition. The saturation pressure of the dead oil with n-pentane was also measured at temperatures 90°C and 120°C. n-Pentane was selected because it is a pure component and, therefore, there is less ambiguity in its contribution to the saturation pressure and asphaltene onset. Also, n-pentane is a liquid at ambient conditions and asphaltene yield data can be readily collected to provide more data for model tuning. Hence, the n-pentane dataset is ideally suited for model tuning.

For the condensate, saturation pressures were measured with the live oil at condensate contents up to 10 wt% and temperatures from 20 to 200°C. The multi-component condensate was selected for practical reasons; it is a cheaper alternative to pure solvents in commercial applications. From a modeling point of view, since it is a multi-component mixture and is combined with live oil over a broad temperature range, it provides a more challenging test of the characterization and modeling methodology.

Bitumen was obtained from the Peace River area. The sample was distilled and centrifuged at the Alberta Research Council to remove water and solids. Vaporized hydrocarbons were condensed and recombined with the bitumen although some light ends were lost. The residual water content was less than 1 wt%. The gas used to make live-oil samples was provided by PraxAir Canada and had the following composition: 37.80 wt% Methane, 61.11 wt% Carbon Dioxide, 1.09 wt% Nitrogen.

Asphaltenes were precipitated from the bitumen using a standard procedure⁽³²⁾. n-Peptane was added to the crude oil at a 40:1 (cm3/g) ratio and sonicated for 45 min. The mixture was left to equilibrate for 24 h. After settling, the supernatant was filtered without disturbing the whole solution. At this point approximately 20% of the original mixture remained unfiltered. Additional n-pentane was added to this solution at a 4:1 ratio of n-pentane to original bitumen mass, sonicated again for 45 minutes and left to equilibrate for 16 hrs. Then the mixture was filtered using the same filter paper. The filter cake (asphaltenes/solids) was washed with approximately 30 ml of n-pentane for 5 days after which effluent color was almost colorless. The filter cake was then dried for three days in a fume hood and for an additional day under vacuum condition to remove remaining n-pentane.

The asphaltenes were fractionated into solubility cuts by dissolving whole asphaltenes in toluene and then adding n-heptane to precipitate a desired fraction of the asphaltenes. A given mass of asphaltenes was dissolved in a specified amount of toluene at 23°C and sonicated for 60–90 minutes to ensure complete dissolution. A specified amount of heptane was added and left to equilibrate for 24 hrs. Then the solutions were centrifuged at 3500 rpm for 5 min. The supernatant was decanted and the solvent was evaporated off to obtain a "light" asphaltene cut. The precipitate was dried to obtain a "heavy" asphaltene cut. A series of fractionations were performed at different heptane-to-toluene ratios to obtain different cuts.

C30+ Compositional Analysis of the Condensate was carried out by Core Lab. A vacuum distillation ASTM D1160 was performed on the bitumen and converted to NBP data using the Edmister method. A modified ASTM D2007 Clay-Gel Adsorption Chromatography⁽²⁵⁾ method (SARA analysis) was also performed. An ASTM D5307 simulated distillation (SimDist) was measured by Core Lab Calgary.

All densities were measured in an Anton Paar DMA HPM density meter and were repeatable within 0.5%. Note, the density of the resin and asphaltene fractions were measured in toluene at different concentrations and the resin or asphaltene density determined by extrapolation⁽¹⁸⁾. Molecular weights were measured in toluene at 23°C using a Jupiter Instrument Model 833 Vapour Pressure Osmometer and were repeatable within 15%.

Live oil samples were prepared using a contactor which consisted of a 600 cm³ horizontal 5 cm diameter cylinder equipped with a piston at each end and a perforated plate placed in the middle of the cylinder in order to enhance mixing of the bitumen and solution gas. Dead oil and solution gas were injected through two injection ports mounted in the middle of the cylinder. Heating tape was used to control the temperature. Two live oil samples were prepared.

To prepare Live Oil 1, 430 g of dead bitumen and 600 cm³ of gas (37.80 wt% Methane, 61.11 wt% Carbon Dioxide, 1.09 wt% Nitrogen) at 17 °C and 4000 kPa were injected in the mixing cylinder. The temperature of the cylinder was maintained at 50°C. The dead bitumen and gas mixture was moved back and forth through the mixing plate for 2–3 days and then the pressure was reduced to 2100±50 kPa and the system was equilibrated at 18±1 °C for 10 days. After equilibration, any excess gas was displaced slowly at constant pressure. After removing the excess gas, the live oil sample was transferred to a sample cylinder. A small amount of the live was taken in pycnometer and single stage flash at 17 °C, 101 kPa was carried out in a Jefri Gasometer to determine its GOR, Table 1.

To prepare Live Oil 2, the same amounts of dead bitumen and solution were mixed using the same procedure except that the system was equilibrated at 18±1 °C and 2600±50 kPa. The GOR and composition data are given in Table 1. Note that Live Oil 1 was erroneously equilibrated at the lower pressure but the data are still suitable for modeling.

Since the live oils were prepared with an excess of solution gas and the solubility of methane and carbon dioxide in the bitumen are different, the composition of the dissolved solution gas is not the same as the original gas composition. Therefore, a GC analysis was performed by Core Labs Calgary on the Live Oil 2 solution gas recovered from the Gasometer, Table 1. The sample for Live Oil 1 was contaminated when sent for testing and could not be analyzed. However, since both live oils were recombined at the same temperature and similar pressures, the solution gas compositions are expected to be the same.

Apparatus: Saturation pressures measurements were carried out in a DB Robinson Jefri PVT cell placed in an air-bath which controls the temperature to within ±0.1 °C. The maximum pressure rating for the PVT Cell is 69 MPa and it can operate at temperatures from −15 to 200°C. The maximum capacity of the PVT cell is 100 cm³. The PVT cell includes a sight glass and is equipped with magnetically coupled impeller mixer mounted on the bottom end cap. The volume of the sample fluid in the cell is determined using a calibrated cathetometer which is precise to ±10⁻⁵ m. The volume of cell, and hence the pressure of the sample fluid under investigation, were controlled by a variable volume computer-controlled positive displacement pump which allowed for the injection or removal of hydraulic oil.

Procedure: A detailed procedure is provided elsewhere⁽⁹⁾ and is summarized below. Before sample injection, the PVT cell was cleaned and vacuumed and the temperature of the air bath was set to the test temperature and allowed to equilibrate. Heavy oil was injected into PVT cell as a liquid phase well above the saturation pressure. If required for the experiment, solvent was then injected also as a liquid. The volumes of sample injected were determined from cathetometer readings plus the dead volume of the PVT cell. The cathetometer readings were verified with the measured pump displacements which agreed within the error of measurement (±0.2 cc). The mass of fluid injected was determined from the volume and the known density at the injection pressure and temperature.

To measure a saturation pressure, the air bath was first equilibrated at the test temperature and the sample fluid was compressed and equilibrated at a pressure well above the bubble point. Then, the pressure was decreased with step-wise volume expansions. After each step, sufficient time was given for the fluid to equilibrate until the pressure was stable. The magnetic stirrer was used to accelerate the equilibration but was turned off for the last 2 hours of equilibration. In general, the equilibration time required in the single liquid phase region and two-phase region were approximately 8 and 24 hours, respectively at room temperature condition while at higher temperature it is 2 and 6 hours, respectively.

The volume was measured using cathetometer after each step and the specific volume was determined using the known mass of the fluid in the cell. Then, the bubble point can be identified from the change in slope in the pressure-volume plot, Figure 1.

The onset of asphaltene precipitation is here defined as the composition (wt% solvent) at which asphaltenes first precipitate. For dead oils and low volatility solvents at ambient conditions, the asphaltene onset condition can be determined by measuring asphaltene yield from bitumen diluted with solvent⁽¹⁷⁾ and extrapolating the yield curve to zero yield. A series of solutions of solvent and bitumen with different solvent contents were prepared in centrifuge tubes from specified masses of bitumen and solvent. The solutions were sonicated for one hour and left to settle for a day to ensure complete precipitation. The mixture was then centrifuged at 3500 rpm for 6 minutes. The supernatant was decanted and the residual asphaltenes were washed with the same solvent and dried. The mass of dry asphaltenes was measured to determine the yield of precipitate (mass of precipitate divided by mass of bitumen) for each solution.

For live oils or volatile solvents where the pressure is above ambient pressure, the onset of the asphaltene was measured using a high pressure microscope (HPM) apparatus.

Apparatus: The apparatus includes the HPM Cell, two mixing cylinders equipped with floating pistons and magnetic stirrers, all placed in an air bath, Figure 2. The HPM cell consists of two sapphire windows with a gap between the windows that can be adjusted with the range of 100 to 400 microns. On one side of the HPM cell there is light source and on the other side there is high focal length camera which is connected to the computer so that digital image and video can be captured in the computer terminal. The HPM Cell can operate at temperatures up to 200°C and pressures up to 20,000 psi (138 MPa). The mixing cylinders are connected to the pump and back pressure assembly which allows the sample to be moved back and forth through the microscope.

Sample Preparation: Even though it had been cleaned and dewatered, the bitumen sample still contained very fine solid particles (Silica, clay, inorganic solids) and approximately micron diameter water droplets.

In order to detect the onset of precipitation as exactly as possible, the solids and water were removed using a multistage filtration process at 50°C. The bitumen was displaced through 5 micron, 1, 2 micron and then 0.2 micron silver membrane (4.7 cm diameter) filters at a flow rate of 0.01–0.0001 cc/min. The filtering process removed almost all of the observable particles.

Procedure: Before bitumen injection, the HPM apparatus was cleaned and vacuumed and the temperature in the air bath equilibrated at the test condition. Bitumen was injected into the mixing cylinders and the volume of bitumen injected was determined from the pump displacement reading. The pump displacement was confirmed from the volume of hydraulic oil volume displaced through the back pressure assembly. Then, almost all bitumen was moved into one mixing cylinder. Solvent was then injected from the PVT cell into the same cylinder in a stepwise process at a flow rate of 1 – 2 cc/hr. Note that a high rate of solvent injection must be avoided because it creates a high local solvent concentration above the onset of asphaltene precipitation. The volume of injected solvent was determined at each step from the cathetometer readings verified with pump displacements. The magnetic stirrer was on during injection but turned off for HPM measurements. After injection, the fluid was moved back and forth between the two cylinders until a uniform mixture was observed in the HPM cell. If asphaltene precipitation did not occur, solvent was again injected from the PVT cell and whole process was repeated. An image was captured after each volume step and in some cases a movie was captured as well.

Measured phase behavior data was modeled using the Advanced Peng-Robinson Equation of State (APR EoS)^(28):

where P is the pressure, T is the absolute temperature, R is the universal gas constant, V* is the molar volume with volume translation, a and b are constants related to the attractive and repulsive forces, and α is given by:

where

for w<0.5:

for w>0.5:

and T_R is the reduced temperature (T/T_c), T_c is the critical temperature, ω is the acentric factor.

For pure components, a and b are defined as follows:

where Pc is the critical pressure. The parameters a and b are determined from the following mixing rules:

where x_i is the mole fraction of component i and a_iJ and b_iJ are given by:

where k_ij is a concentration independent interaction parameters between components i and j.

To use an equation of state model for a petroleum fluid, the fluid must be divided into components and pseudo-components with defined mole fractions, densities, critical properties, and acentric factors. The model is tuned to fit experimental data by adjusting the binary interaction parameters. In this work, the interaction parameters were determined from the following correlation:

For pairs of pseudo-components the Gao et al. correlation⁽²⁹⁾ was used, in which the exponent n in Eq. 10 is set to Z_Cij, the critical compressibility factor. For solvent and pseudo-component pairs, n was adjusted to tune the APR EoS to fit the experimental phase behavior data.

Data:

Table 2 is the GC assay of the condensate. The average molecular weight and specific gravity of the C6+ fraction are 111 g/mol and 0.744 respectively. The amount of the C30+ fraction was neglible. Figure 3 shows the normal boiling point (NBP) curve of the Bitumen and Table 4 presents the SARA analysis. The density of the bitumen was 1014.9 kg/m³ at 19.7 °C and the densities of the SARA fractions are reported in Table 4. The molecular weight of the maltenes was 500±50 g/mol.

To aid in the characterization, the asphaltenes were fractionated into the solubility cuts given in Table 3. The densities of each cut at 23°C and atmospheric pressure are also provided in Table 3. The molecular weight of each cut was measured as a function of asphaltene concentration and all but the lightest cut self-associated; that is, the molecular weight increased with concentration. The molecular weight of the lightest cut was 840 ±50 g/mol

The dead bitumen was characterized as recommended by Castellanos-Diaz et al.^{(13, 16)}. Maltenes and asphaltenes are characterized separately. For maltenes, the distillation data of Figure 3 was extrapolated over the maltene fraction (first 80.4 wt% of the NBP curve) using a Gaussian distribution of boiling points. The NBP curve of maltene was divided into ten pseudo-components: six for a light oil section (from 200 to 427 °C), three for a medium oil section (from 427 to 580 °C), and one for a heavy oil section (580 to 600 °C). The molecular weight and specific gravity of each pseudo-component were calculated using Lee-Kesler and Katz Firoozabadi correlations, respectively. The calculated molecular weight and specific gravity of each pseudo-component were adjusted by a multiplier to fit the average values of the maltene (495 g/mol and 982 kg/m³).

Asphaltenes self-associate and a Gamma distribution of the aggregate size distribution has been found to give good predictions for asphaltene yields with regular solution based models^{(17, 18)}. The Gamma distribution function is given by:

where MW_m, and are the monomer molar mass and average ‘associated’ molar mass of asphaltenes, and β a parameter which determines the shape of the distribution. The lightest cut from the asphaltene fractionation did not appear to self-associate and its molecular weight was taken as the average monomer molar mass. It measured molecular weight was 840±50 g/mol. A value of 820 g/mol along with the measured density of 1090 kg/m³ give an NBP of 585°C using the Lee-Kesler correlation which is consistent with the end point of the extrapolated maltene NBP. Once the monomer properties were fixed, β was set to 2.5 and MW_Avg was set to 1400 g/mol as was done by Castellanos-Diaz et al.^{(13, 16)}. The asphaltenes were then divided into 6 pseudo-components of equal weight fraction.

The specific gravity of the asphaltenes was determined from the density data presented in Table 3. The following distribution was found to fit the data:

where W is cumulative weight fraction of the asphaltenes and SG is the specific gravity. The distribution was integrated to determine the density of each asphaltene pseudo-component. The average boiling point of each pseudo-component was then calculated from the Lee-Kesler correlation. The complete boiling point curve is shown in Figure 3. Finally, the critical properties and accentric factor of all of the pseudo components were determined using Lee-Kesler correlation. The pseudo-component properties are summarized in Table 5.

SimDist was also collected for the sample. Figure 3 shows that the SimDist is similar to the NBP from vacuum distillation for the first 10 wt% of the sample but deviates for heavier cuts. Extrapolation of SimDist data over maltene fraction using Gaussian distribution gives a very high boiling point value for asphaltene monomer, Figure 3. This characterization leads to the incorrect prediction of two liquid phases for bitumen at ambient conditions. Therefore, SimDist data were discarded in this study.

The preparation of Live Oil 2 was simulated using VMGSim. The exponent n of Eq. 10 for the binary interaction between methane, CO2, and nitrogen with the pseudo-components was adjusted in order to match the solution gas composition and GOR given in Table 1. The fitted values of n are −0.3, 0.75, and 0.1 for methane/pseudos, CO₂/pseudos, and N₂/pseudos. The negative n value for methane/pseudos gives negative binary interaction parameters, consistent with literature⁽³⁰⁾.

The preparation of Live Oil 1 was modeled without further adjustment and the predicted GOR was within experimental error of the measured value. As expected, the predicted composition of the Live Oil 2 solution gas was almost identical to the composition of the Live Oil 1 solution gas.

Based on the boiling point range from the C30+ analysis in Table 2, an NBP curve of the C6+ fraction was generated. Then, the NBP curve was divided into pure components up to pentane and eight pseudo-components for the C6+ fraction: four for a light oil section (from 37 to 90 °C); two for a medium oil section (from 90 to 200 °C), and two for a heavy oil section (200 to 450 °C). A minimum of eight pseudo components was required to accurately represent the TBP curve. The molecular weight and specific gravity of the pseudo-components were calculated using Lee-Kesler and Katz Firoozabadi correlations, respectively, and then adjusted by a multiplier to fit the average values of the C6+ fraction. The critical properties and accentric factor of the each pseudo component were calculated using Lee-Kesler correlations. The properties of pseudo and pure components of the condensate are summarized in Table 7. Binary interaction parameters between CO2 and N2 with the condensate pseudo-components were calculated using Nishiumi-Arai-Takeuchi correlation⁽³¹⁾, while those between methane and the pseudo-components were calculated using the Gao et. al.⁽²⁹⁾ correlation.

Figure 4 shows the saturation pressures of a mixture of 11 wt% pentane and 89 wt% bitumen at temperatures of 90 °C and 120 °C. Figure 5 shows that the onset of asphaltene precipitation at 23°C and atmospheric pressure is between 46 to 48 wt% of n-pentane. The onset composition from the High Pressure Microscope is in good agreement with yield data which extrapolate to an onset of 45 ±5 wt% n-pentane, Figure 6.

The APR EoS was fitted to both the saturation pressures and the onsets of asphaltene precipitation simply by adjusting the exponent n in Eq. 10 to 0.62 for the interaction parameters of n-pentane and the pseudo-components. Figures 4 and 6 demonstrate that this single oil characterization allows us to fit both the saturation pressures and the asphaltene precipitation onsets within the accuracy of the data. However, the model under-predicts asphaltene yields at high pentane contents. Measured asphaltene yields increased monotonically while the model predicted that the asphaltenes become soluble at high dilution. The predicted increase in solubility at high dilution occurs because the model over-predicts the amount of non-asphaltene components that partition to the asphaltene-rich phase, Figure 6. At high dilutions, so much n-pentane is predicted to partition to the incipient asphaltene-rich phase that phase separation does not occur. Asymmetric mixing rules that minimize the amount of solvent that can partition to the asphaltene-rich phase could provide a more accurate model.

Figure 7a shows the saturation pressures of Live Oil 1 at temperatures from 22 to 180°C. Figure 7b shows the saturation pressures of Live Oil 2 with 6 and 10 wt% condensate. Figure 7b shows that the saturation pressure of the live oil decreases upon addition of the condensate. The addition of condensate decreases the concentration of the methane, carbon dioxide, and nitrogen which have a significantly higher vapour pressure than the condensate. Hence, the net effect is to decrease the saturation pressure of the mixture.

The saturation pressures were modeled with the same interaction parameters used to model the live oil preparation (n=−0.3 for methane-pseudos and n=0.75 for CO₂-pseudos). As shown by the dotted lines in Figures 7a and 7b, the model over-predicts the saturation pressures. Therefore, the following temperature dependent n values were used to fit the saturation data for Live Oil 1:

where T is in °C. The solid lines on Figure 7a show the fit to the data.

Then, the data for the mixtures of Live Oil 2 and condensate were predicted with the same n values. The solids lines in Figure 7b show that the model predicted the effect of the added condensate within the accuracy of the data up to 120°C. The model over-predicted the saturation pressures at higher temperatures. It is likely that the light components in the condensate, such as n-pentane, also require temperature dependent n-values. Saturation pressure data for pure component/bitumen mixtures over a range of temperatures are required to confirm this hypothesis.

A characterization based on extrapolated distillation data for maltenes and a separate molecular weight based characterization for asphaltenes was used with the Advanced Peng Robinson equation of state to model both saturation pressures and the onset of asphaltene precipitation. The model under-predicted asphaltene yields but this deficiency may be overcome with asymmetric mixing rules.

The same approach was used to model saturation pressure of live bitumens and mixtures of live bitumen with condensate. Temperature dependent interaction parameters were required to fit the data. Saturation pressure data for pure hydrocarbons and bitumen over a wide range of temperatures are required to verify the temperature dependence.

Characterizations based on extrapolated SimDist data led to unrealistic phase behaviour predictions because the boiling points extrapolated to a high end point creating a self-incompatible bitumen model.

The authors are grateful to the sponsors of the NSERC Industrial Research Chair in Heavy Oil Properties and Processing, including NSERC, Shell Energy Ltd., Schlumberger, and Petrobras for financial support. We thank Diana Ortiz for carrying out SARA analysis, vacuum distillation, and molecular weight measurements; Diana Barrera for measuring the density of asphaltene fraction; and Hamed Reza Motahhari for preparing the live oil samples and measuring their density. We also thank Virtual Materials Group for providing simulation software.