Impacts of Kinetics Scheme Used To Simulate Toe-to-Heel Air Injection (THAI) in Situ Combustion Method for Heavy Oil Upgrading and Production.

While simulating toe-to-heel air injection (THAI), which is a variant of conventional in situ combustion that uses a horizontal producer well to recover mobilized partially upgraded heavy oil, the chemical kinetics is one of the main sources of uncertainty because the hydrocarbon must be represented by the use of oil pseudo-components. There is, however, no study comparing the predictive capability of the different kinetics schemes used to simulate the THAI process. From the literature, it was determined that the thermal cracking kinetics schemes can be broadly divided into two: split and direct conversion schemes. Unlike the former, the latter does not depend on the selected stoichiometric coefficients of the products. It is concluded that by using a direct conversion scheme, the extent of uncertainty imposed by the kinetics is reduced as the stoichiometric coefficients of the products are known with certainty. Three models, P, G, and B, each with their own different kinetics schemes, were successfully validated against a three-dimensional combustion cell experiment. In models P and G, which do not take low-temperature oxidation (LTO) into account, the effect of oil pseudo-component combustion reactions is insignificant. For model B, which included LTO reactions, LTO was also found to be insignificant because only a small fraction of oxygen bypassed the combustion front and the combustion zone was maintained at temperatures of over 600°C. Therefore, in all the models, it is observed that coke deposition was due to the thermal cracking taking place ahead of the combustion zone. During the first phase of the combustion, peak temperature curves of models P, G, and B closely matched the experimental curve, albeit with some deviations by up to 100°C between 90 and 120 min. After the increase in the air injection flux, only the model P curve overlapped the experimental curve. The model P cumulative oil production curve deviated from the experimental one by only a relative error of 4.0% compared to deviations in models G and B by relative errors of 6.0 and 8.3%, respectively. Consequently, it follows that model P provided better predictions of the peak temperature and cumulative oil production. The same conclusion can be drawn with regard to the produced oxygen concentration and combustion front velocity. With regard to American Petroleum Institute (API) gravity, it is found that all the three models predicted very similar trends to the experiment, just like in the case of the oil production rate curves, and therefore, no model, in these two cases, can be singled out as the best. Also, all the models' predictions of the produced CO X concentration prior to the increase in the air flux closely match the experimental curve. There are, however, serious differences, especially by model P, from the reported experimental curve by up to 15% after the increase in the air flux.

Introduction

Conventional in situ combustion (ISC) involves air injection into an oil reservoir to oxidize the immobile carbonaceous fraction of the bitumen/heavy oil in place. An advancing and expanding combustion front is created and sustained with the continuous air injection via a vertical injector well. The heat generated during the combustion results in substantial viscosity reduction and significant upgrading of the heavy oil. The mobilized oil has to travel over several hundreds of meters before reaching the vertical production well. Toe-to-heel air injection (THAI) is a variant of conventional ISC that uses the horizontal well for heavy oil mobilization and production. In THAI, the combination of heat from combustion reactions, mass and momentum transfer, and gravity-assisted drainage is used to mobilize the heavy oil to the surface.[1−7]The combustion front propagates continuously from the toe to the heel of the horizontal producer (HP) well. During the process, typical temperatures in the vicinity of the combustion region range between 500 and 700°C.[1] The high temperature is similar to that reported by Petrobank during their WHITESAND’s field pilot project.[4] Unlike in the case of conventional ISC, the mobilized oil is continuously produced as it gravity-drained and as it needs not bank and/or travel over hundreds of meters. Therefore, THAI automatically falls into the category of the short distance gravity drainage oil recovery and upgrading techniques.[2] The full description of the advantages of the THAI process is given elsewhere.[5−7]

The air injection enhanced oil recovery technique is usually studied by the use of either the combustion tube[8−14] or three-dimensional (3D) combustion cell.[1,15,16] The experiment allows the evaluation of the likely mechanism of the physicochemical processes at the field scale. Numerical models are then developed and validated against the experimental results. The model is then upscaled for field scale evaluations.[17] However, the kinetics of fuel deposition and combustion reactions required to simulate the air injection enhanced oil recovery process remains one of the main sources of uncertainty. Not only that, oil pseudo-components must be used due to the complex nature of hydrocarbons and the impracticality of developing kinetics based on the total number of the individual compounds making up the heavy oil. In addition to that, the main mechanism through which fuel for the sustenance of high-temperature oxidation (HTO) is deposited remains a contentious issue. Some authors consider fuel deposition to be mainly as a result of low-temperature oxidation (LTO) of the heavy oil,[13,18] while others consider it to be due to thermal cracking of the heavier fraction ahead of the combustion zone.[8,9,19] Alexander et al.,[20] however, observed that fuel deposition via LTO is only significant during the start-up period (i.e., prior to coke combustion). Therefore, the aim of this work is twofold: first to present a brief review about the kinetics schemes used when simulating ISC and the second is to present a comparative study of the numerical simulation results of three different kinetics schemes validated against experimental results. It should be noted that the rock–fluid properties could also be one of the main sources of uncertainty when simulating the THAI process. However, because they are not the subject of this work, they are not considered here.

Chemical Kinetics

Direct Conversion Thermal Cracking Kinetics

Thermal cracking reaction kinetics schemes, which are very similar in terms of general representation and proposed for implementation when simulating the air injection enhanced oil recovery process, are shown in Figure .[18,19,21−24] In all the schemes except that of Wiehe,[22] coke was considered to be formed from asphaltene cracking after a certain induction period which depends on the initial asphaltene content of the oil.[22,25,26] In addition to the formation from asphaltene, the coke was considered to be formed from thermal cracking of the maltene pseudo-component by Gray et al.[23] Radmanesh et al.[24] also considered the thermal cracking of both the light and heavy oil pseudo-components to result in coke formation (Figure and Table ). The number of pseudo-components considered by Radmanesh et al.[24]as been reduced by lumping the “distillates” and “gas oil” pseudo-components to form “light oil” in this study. In the case of Wiehe,[22] in which thermal cracking kinetics of the Cold Lake vacuum residue was developed, the coke generation was considered to be the result of the formation of supersaturated solution of the “asphaltene core”. The coke formed only when the solubility limit of the “asphaltene core”, which was formed directly from the heavy oil, was reached. Phillips et al.[19] developed two sets of cracking kinetics schemes (referred to “model A” and “model B”) for Athabasca bitumen. In one of the models, (i.e.,“model A” which is shown in Figure ), the formation of light pseudo-components from “asphaltene” was considered insignificant based on the assumption that the light oil pseudo-component is formed from thermal cracking of the heavy oil pseudo-component. The thermal cracking kinetics data fitted excellently in “model A” compared to the fit obtained in the other model (i.e.,“model B”) which had six pseudo-components. These types of thermal cracking kinetics can be described as direct conversion kinetics because the formation of any pseudo-component has a different set of kinetics parameters and the stoichiometric coefficient of the products is easily determinable.

Figure 1

General thermal cracking scheme showing coke formation from asphaltene, maltenes, and light oil.

Table 1

Reaction Scheme Used by the Different Authors as Depicted in Figure

authors	kinetics scheme used	comment
Phillips et al.[19]	1,2,3,4,8	“model A” as defined by the authors
Belgrave et al.[18]	1,6,8	LTO reactions have not been included
Adegbesan et al.[21]	1,8	“resins” and “asphaltene” are lumped
Wiehe[22]	1,2,3,6,8	“asphaltene” and “asphaltene core” are lumped
Gray et al.[23]	2,5,6,8	only liquid phase cracking considered
Radmanesh et al.[24]	2,3,5,7,8	reaction no. 6 has zero stoichiometry

Split Conversion Thermal Cracking Kinetics

Split thermal cracking kinetics can be described as those schemes that consider coke formation to be the result of conversion of the heavy oil pseudo-component into coke and light oil pseudo-components. A single frequency factor, as well as activation energy, for the formation of the two (or more) pseudo-components is associated with this kind of scheme. The scheme does not take into account the induction period before coke formation. A significant number of simulation studies have used this kind of thermal cracking reaction scheme.[9,14,26−31] A typical general representation is shown in Figure . Quite good matches were obtained using this scheme. However, the main disadvantage with such a scheme is its heavy dependence on the selected stoichiometric coefficients of the products.

Figure 2

General thermal cracking scheme showing coke and light oil formation from the heavy component.

Low-Temperature Oxidation

LTO is an oxygen addition reaction which results in an increase in the asphaltene content of heavy oil. The increase in the asphaltene content results in an increase in the overall oil viscosity. Operating in an LTO mode, which usually takes place in the temperature region of less than 380°C,[32−34] has been reported to result in poor combustion propagation because of restricted gas distribution.[3,35] This is why an adequate air injection rate must be maintained during combustion in order to not allow the system to switch to the LTO mode. The kinetics of Athabasca bitumen developed by Belgraveet al.[18] considered coke formation to be due to both thermal cracking (Figure and Table ), as described earlier, and LTO reactions. In the LTO reactions, “maltenes”combined with oxygen to form “asphaltenes”and “asphaltenes” in turn reacted with oxygen to form coke. A typical LTO reaction is shown in Figure . For their model to be accurately accounting for the oxygen atoms that combined with partially oxidized “asphaltenes” to form solid coke (CH1.3), gaseous products must be produced directly from the “asphaltenes”. This was represented by the thermal cracking reaction 1 shown in Figure .

Figure 3

Typical representation of LTO reaction.

Jia et al.[36] developed a kinetics model of thermal cracking and LTO reactions based on Athabasca bitumen. Their model quite accurately predicted the thermal cracking products obtained by Hayashitaniet al.[37] The induction period before coke formation took place was also closely predicted. An experimental study of LTO of Athabasca bitumen performed by Millouret al.[32] showed that more coke was deposited when bitumen was used as the starting material instead of any of the pseudo-components. The coke formation was observed to take place after some induction period. They also observed that at temperatures more than 175°C, the LTO reaction order with respect to oxygen partial pressure tended to zero, implying that once a certain high temperature is reached, LTO becomes insignificant and hence has no effect on the process. In another study carried out by Adegbesanet al.[21] on LTO of Athabasca bitumen, the overall oxygen consumption kinetics together with four different thermal cracking reactions were developed. The first scheme showed that bitumen combined with oxygen in the LTO region to produce products. However, no information about what the products were was given. A good match was however obtained when one of the thermal cracking kinetics, which did not take the formation of coke from the LTO reaction into account, was validated against experimental results.

High-Temperature Oxidation

In HTO, coke is combusted in the presence of oxygen to produce carbon dioxide, carbon monoxide, and water. A typical HTO representation is shown in Figure below. A detailed review about the reactions involved during in situ combustion has been given elsewhere.[38,39]

Figure 4

Typical representation of a HTO reaction.

Another set of reactions often incorporated into the numerical simulations of the in situ combustion process are the individual pseudo-component combustion reactions. Typically, the oil pseudo-component is oxidized to carbon oxides and water. Greaves et al.[29] and Marjerrison and Fassihi[40] have found that the effect of these reactions on the simulation results is negligible. This is because oxygen has to bypass the combustion front before oil is burned. However, authors such as Lin et al.[9] and Anaya et al.[28] have included them in their model and their influence have not been reported.

Overall, what could be deduced from these different studies is that the formation of coke, either due to thermal cracking or as a result of LTO, is not an instantaneous process but requires an induction period. The effect of LTO as a mechanism of fuel deposition as applied to the THAI process has not been investigated. Moreover, in the literature, to the best of my knowledge, two different kinetics schemes were validated against the THAI experiment reported in Xia and Greaves.[1] The first used the split conversion[29] and the second used the direct conversion[5] thermal cracking schemes, respectively. No comparison between the predictions of each kinetics scheme was however made. Given that the kinetics is one of the main sources of uncertainty when simulating the THAI process, it is the aim of this work to investigate how closely each of the three different kinetics schemes will history-match the 3D combustion cell experimental results. The differences and similarities among the predictions by the three different kinetics schemes are discussed. The first numerical simulation model (i.e., model P) used Phillips et al.[19] direct conversion thermal cracking kinetics scheme, which can also be found in Rabiu Ado et al.[5] and Ado,[41,42] the second used modified Greaves et al.[29] split conversion kinetics scheme (i.e., model G), while the third used Belgrave et al.[18] direct conversion kinetics scheme with LTO reactions for fuel deposition and upgrading (i.e., model B). The two former schemes (models P and G) included combustion of oil pseudo-components together with HTO combustion reactions while the latter considered HTO combustion reaction only.

Results and Discussion

After successful validation of the three different models against the 3D combustion cell experimental result, the prediction of peak temperature, oil production rate, cumulative oil production, American Petroleum Institute (API) gravity, produced oxygen mole percent, CO (i.e., combination of CO2 and CO), and combustion front velocity by each model are compared.

Peak Temperature

A close match between the simulated and experimental peak temperature is obtained over the significant portion of the dry combustion period (Figure ). All the peak temperature predictions by the three models are better than that by Greaves et al.[29] However, the closest match is achieved with the model P with the predicted peak temperature overlapping the experimental one during most of the dry combustion period. Models G and B deviated from the experimental peak temperature toward the end of the dry combustion period just like in the case of Greaves et al.[29] The predicted spike in the peak temperature by models P and G, which occurred when the combustion front most likely reached the toe of the HP well and consumed the large amount of coke, only lagged the experimental one by only 20 min. This is a significant improvement when compared to the prediction by the model of Greaves et al.[29] here the discrepancy is more than six times that observed in the three new models (Figure ). Furthermore, it should be noted that model B predicted multiple spikes on the peak temperature (Figure ) from around 170–220 min. This indicated that there were different grid-blocks around the toe of the HP well that contained different large coke amounts.

Figure 5

Peak temperature vs time.

All the models also, dynamically, respond to the increase in the air injection rate, as reflected in the increase in peak temperature at 190 min. Because, during the experiment and as predicted by all the models, the peak temperature is maintained above 600°C and only small concentration of oxygen is observed to bypass the combustion front in the case of model B, LTO reactions were found to be insignificant. Thus, all the fuel deposition as predicted by each model was due to thermal cracking ahead of the combustion front.

Oil Production Rate

As the bitumen has no mobility prior to preheating, oil production begins only after the first 18 min of the pre-ignition heating cycle (PIHC) (Figure ). All the models, including Greaves et al.,[29] over predicted the oil production rate over the 18–30 min period by 12–14 cm3min–1. The over prediction was observed to be due to pressure build-up which is caused by significant thermal cracking of oil around the preheated zone. Prior to the increase of air injection flux by 33% (i.e., to 16 m3 m–2·h–1), the trends in the predicted and experimental oil production rates matched closely. The Greaves et al.[29] oil production rate curve slightly lies above the experimental curve from 190 to 280 min. Thereafter, it lies below the experimental curve up to the end of the dry combustion period. The decrease in the oil production rate corresponds to the period over which oxygen production was predicted to take place, as will be seen in Figure . The oil production rate curves predicted by models G, P,and B lie below the experimental curve from 190 min to the end of the dry combustion period. On the whole, there is, however, a good agreement between the models predictions and the experimental oil production rate, over most part of the combustion period.

Figure 6

Oil production rate vs time.

Figure 9

Produced oxygen concentration vs time.

Cumulative Oil Production

All the models predicted a cumulative oil production of 4% of oil originally in place at the end of the 30 min of the PIHC and thus exactly matching the experimental value (Figure ). It has been noted, however, that the trend in the predicted curves did not follow the experimental one, which is caused by the delay in oil production (Figure ). This is because the mobilized oil has to reach the HP first before being produced on the surface. On the overall general trend, the closest match was obtained with model P, which under predicted the experimental cumulative production by an overall relative error of 4% (Figure ). On the other hand, models G and B under predicted the cumulative oil production at the end of the dry combustion period by overall relative errors of 6.0 and 8.3%, respectively. This means a significant improvement over the prediction by the Greaves et al.[29] model, which deviated from the experimental cumulative oil production by a relative error of 7.2%, is realized by models P and G.

Figure 7

Cumulative oilproduction vs time.

Oil Upgrading

The API gravity gives the measure of the extent to which the bitumen is upgraded. The general trend of the experimental API gravity is closely predicted by all the three models (Figure ). However, a significant deviation by up to 4 API points between the experiment and the predictions can be observed over the time period of 80–140 min. The deviation could be due to the fact that the API gravity measured during the experiment is that of oil collected over every 15 min period (sampling takes place every quarter of an hour during the experiment). Therefore, unlike in the case of the numerical predictions in which the API gravity of produced oil over a fraction of a minute is reported, the variation in the experiment was ironed out because of the low frequency of sampling. After 150 min, all the model predictions closely match the experiment up to the end of the dry combustion period. From Figure , it can be determined that all the models predicted a very similar trend in the API gravity and no particular model, just like in the case of the oil production rate curves (Figure ), can be singled out as the best in terms of API gravity predictions.

Figure 8

API gravity vs time.

Oxygen Concentration

An accurate prediction of oxygen utilization is critical to the safety and economics of the THAI process. The experiment showed that oxygen production started prior to the increase in the air injection flux (Figure ). This is contrary to the explanation provided by Greaves et al.[29] that the oxygen production was due to an increase in the air injection flux. One of the major improvements realized over Greaves et al.[29] model is that the duration when oxygen production began is excellently predicted by models P and G. However, in the case of Greaves et al.[29] prediction, a serious delay by up to 105 min can be observed (Figure ). The trend predicted by model G is very similar to the experimental trend over time periods of 175–275 min. Thereafter, the concentration of the produced oxygen dropped from 1.00 mol% to about 0.86 mol% up to 315 min (Figure ). The increase in oxygen utilization is reflected in the corresponding increase in peak temperature after 250 min. The trend predicted by model P is essentially similar to the experimental one. However, there is the occurrence of regular periodic cycles in the oxygen concentration curve P which is as a result of propagating combustion along the HP in zones of alternating low and high coke concentrations. However, observing Figure closely, it can be seen that all the models predicted periodic cycles in the oxygen concentration curves. Similarly, there are multiple minimum and maximum points on the oxygen concentration experimental curve, which are thought to be replicated by all the models. However, the regular periodic cycles in model P are more pronounced compared to the other models. The difference between the periodic cycles in model P and the experiment could be explained by the fact that the flue gas concentration was sampled discretely during the experiment and as a result, the rate of sampling might not be sufficient to display pronounced periodicity like that shown in model P. Because alternated low and high coke concentrations along the HP well are thought to cause the periodicity in all the models and the experiments, it can be explained that the alternating low and high coke concentrations along the HP are due to the variability (periodic cycles of crest and trough) in the rates of oil production (see Figure ) as the combustion and the mobile oil zone propagated in a toe-to-heel manner. It could be that when there is a decrease in the oil production rate (i.e., trough in the curve), more of the oil is available to be converted to coke and when there is an increase in the oil production rate (i.e., crest in the curve) less of the oil is available to be converted to coke. In the case of model B, it can be observed that the oxygen production began at around 40 min and the concentration remained below 0.25 mol% up to 220 min. The cause of early breakthrough was observed to be as a result of air bypassing the HTO combustion zone. However, despite that, LTO is found to be insignificant. After 220 min, the combustion front at the toe started to propagate along the HP. This is why the concentration of the produced oxygen increased up to a maximum of 0.7 mol%. It follows that model B deviated the most from the experiment.

Produced CO Concentration

During the first phase of the dry combustion when the air flux was 12 m3 m–2·h–1, an excellent agreement between the models predictions and the experimental CO composition is achieved (Figure ). However, as the air flux was increased to 16 m3 m–2·h–1, all the models predictions remained between 17 and 18 mol% which is a deviation from the experimental curve by a relative error of 10–15%. This is attributed to the fact that the CO2/(CO2 + CO) ratio and thus the stoichiometric coefficients of the combustion reactions were kept constant in all the models. Therefore, it is acknowledged that all the models are incapable of replicating the dynamic increase in the CO concentration after the increase in the air flux.

Figure 10

Produced COconcentration vs time.

Combustion Front Velocity

Figure shows the graph of the distance against time for the three different models and for the two different air injection fluxes of 12 and 16 m3 m–2·h–1. At an air flux of 12 m3 m–2·h–1, models P and G predicted essentially the same combustion front velocity of 0.026 m·h–1 while model B predicted 0.035 m·h–1. As the air flux was increased to 16 m3 m–2·h–1, models P and B showed an increase in the combustion front velocity as reflected by the increase in the slope of the distance–time graph. However, because of high fuel availability, no increase in the combustion front velocity is observed in model G. The average combustion front velocity predicted by model P is 0.039 m·h–1 which deviated from the experimental one of 0.035 m·h–1 by a relative error of 12%. Models B and G predicted an average of 0.041 and 0.026 m·h–1 which imply deviations from the experiment by relative errors of 17 and 25%, respectively. This shows that model P provides a better prediction of the experimental result.

Figure 11

Distance–time graph for the determination of the combustion front velocity along the vertical middle plane and at the top of the combustion cell as predicted by the three models.

Conclusions

From the brief literature review carried out, it was determined that the thermal cracking kinetics schemes can be broadly divided into two: the split conversion thermal cracking scheme which heavily depends on the selected stoichiometric coefficients of the products and direct conversion thermal cracking scheme which does not depend on the stoichiometric coefficients of the products. The split conversion thermal cracking kinetics scheme depends heavily on the stoichiometric coefficients of the reaction’s products and is of the formsuch that, there are infinite numbers of possible combinations of x and y, and hence, the scheme is highly indeterminable, while the direct conversion thermal cracking kinetics scheme does not depende on the stoichiometric coefficient of the reaction’s products and is of the formsuch that x has a single value and is uniquely determined; similarly y has a single value and is uniquely determined. It is therefore concluded that by using the direction conversion thermal cracking kinetics scheme, the extent of uncertainty imposed by the kinetics is reduced as the stoichiometric coefficients of the products are known with certainty. It was also determined that the mechanism through which fuel deposition takes place was a contentious issue as a section of authors considered it to be as a result of thermal cracking taking place ahead of the combustion front while the other section considered it to be due to LTO reactions.

Three different kinetics schemes have been successfully validated against a 3D combustion cell experiment. In models P and G which do not take LTO into account, it is found that the effect of the individual oil pseudo-component combustion reaction is insignificant. For model B which included LTO reactions, it is shown that LTO was insignificant because only a small fraction of oxygen bypassed the combustion front and the combustion zone was maintained at temperatures of over 600°C. Therefore, in all the models, it is observed that fuel deposition was as a result of thermal cracking taking place ahead of the combustion zone. This is a similar conclusion drawn by Yang et al.[13]

During the first phase of the combustion, peak temperature curves of models P, G, and B closely matched the experimental curve, albeit with some deviations by up to 100°C between 90 and 120 min. After the increase in the air injection flux, only the model P curve overlapped the experimental curve during most of the second phase of combustion. It is therefore concluded that model P provided the best prediction of peak temperature compared to models G and B.

The model P cumulative oil production curve deviated from the experimental curve by only a relative error of 4.0% compared to deviations in model G and B by relative errors of 6.0 and 8.3%, respectively. Consequently, it follows that model P provided a better prediction of the cumulative oil production. The same conclusions can be drawn with regard to the produced oxygen concentration and combustion front velocity.

With regard to API gravity, it is found that all the three models predicted very similar trends to the experiment, just like in the case of the oil production rate curves, and therefore, no model, in these two cases, can be singled out as the best. Also, all the models’ predictions of the produced CO concentration prior to the increase in the air flux closely match the experimental curve. There are, however, serious differences, especially by model P, from the reported experimental curve by up to 15% after the increase in the air flux.

Methodology

A 3D model based on the THAI experiment conducted by Xia and Greaves[1] was developed using Computer Modelling Group’s thermal reservoir simulator (CMG’s STARS). The commercial thermal simulator allows the incorporation of a discretized wellbore as a representation of a HP. This allows the transient nature of the multiphase flow and the heat transport to be modeled by discretization. The resulting algebraic equations, because of discretization, were then coupled with the reservoir reactive transport equations and solved using PARASOL, which is a parallel processing solver available in STARS. The computer used to run the models has two parallel processors each with 8 cores (i.e., 16 cores and, thus, 32 threads in total). However, only 25% of the CPU is used as the maximum number of threads that can be specified to PARASOL is 8.

The physical model comprised the 0.6 m × 0.4 m × 0.1 m combustion cell packed with virgin Athabasca bitumen (Figure ). An electrical heater was used to preheat the inlet zone around the horizontal injector (HI), prior to air injection. During the physical experiment, air was injected at a rate of 8000 Scm3 min–1 (i.e., a flux of 12 Sm3 m–2·h–1). After 190 min, the air injection flux was increased by 33%, to 16 Sm3 m–2·h–1, and maintained up to the end of the dry combustion period (i.e., 320 min).

Figure 12

Well arrangement and dimensions of the 3D combustion cell model.

Pressure, Volume, and Temperature Data

The pressure, volume, and temperature (PVT) data used in models P, G, and B are shown in Tables , 3, and 4, respectively. LC, MC, and IC, which are shown in Table , represent the light, mobile, and immobile pseudo-components, respectively. LITE and HEAV oil, which are shown in Table , represent the light and heavy pseudo-components, respectively. MALT and ASPH, which are shown in Table , represent the maltene and asphaltene pseudo-components, respectively. It should be noted that the PVT data shown in Table are not exactly the same as those provided by Belgrave et al.[18] The change was necessary in order to obtain the best history match possible. The viscosity of the Athabasca bitumen as function of temperature is shown in Figure a. The phase equilibrium K-values for the individual pseudo-components, which are required to account for phase change, were estimated using the Wilson equation.[43] They are presented as a function of temperature and pressure as used in models P, G, and B and are shown in Figure b–d, respectively.

Table 2

PVT Data as Used in Model P as Presented in Rabiu Ado et al.[5]a

components	split(mol%)	RMM (g/mol)	Pc(kPa)	Tc(°C)	ρ (kg/m3)	accentricity	TB(°C)
LC	42.50	210.82	1682.88	464.68	828.24	0.62	281.47
MC	23.91	496.81	1038.46	698.53	961.66	1.18	549.67
IC	33.59	1017.01	729.22	940.36	1088.04	1.44	785.78

Reprinted with Permission from Rabiu Ado et al.[5]

Table 3

PVT Data as Used in the Modified Greaves et al.[29] Model Ga

components	split(mol%)	RMM (g/mol)	Pc(kPa)	Tc(°C)	ρ (kg/m3)	accentricity	TB(°C)
LITE oil	36.47	170.00	2305.95	425.16	903.80	0.48	246.60
HEAV oil	63.53	878.00	1031.29	780.00	1012.07	1.45	711.00

Reprinted with Permission from Greaves et al.[29]

Table 4

PVT Data as Used in Model B Using Modified Belgrave et al.[18] Dataa

components	split(mol%)	RMM (g/mol)	Pc(kPa)	Tc(°C)	ρ (kg/m3)	accentricity	TB(°C)
MALT	70.53	406.70	1350.87	612.88	911.53	0.90	430.00
ASPH	29.47	892.98	818.84	898.42	1139.00	1.59	742.59

Reprinted with Permission from Belgrave et al.[18]

Figure 13

Model input data as a function of temperature: (a) bitumen viscosity[44] used in each of the three models and (b) model P, (c) model G, and (d) model B VLE K-values for the individual pseudo-components.

Thermal Properties of the Oil Pseudo-Components

The PVT data were used in Aspen HYSYS in conjunction with the Peng–Robinson equation of state to obtain the thermal properties of the oil pseudo-components for each model. The equation to obtain both the thermal conductivity α(W/m·K) and the specific heat capacity C (kJ/kmol·K) using the coefficients shown in Tables –10 iswhere T is the temperature in K.

Table 5

Liquid Thermal Conductivity Coefficients for Model P

	coefficients
components	A	B	C	D	E
LC	2.311 × 10–1	–3.484 × 10–4	2.959 × 10–7	–6.768 × 10–11	–2.149 × 10–13
MC	2.140 × 10–1	–1.384 × 10–4	–1.767 × 10–7	3.891 × 10–10	–2.577 × 10–13
IC	1.869 × 10–1	–5.378× 10–5	–2.042 × 10–7	2.760 × 10–10	–1.286 × 10–13

Table 10

Liquid Specific Heat Capacity Coefficients for Model B

	coefficients
Components	A	B	C	D	E
MALT	59.50	2.704	–1.263 × 10–3	–1.243 × 10–14	5.879 × 10–18
ASPH	–155.2	5.825	–2.719 × 10–3	–1.058 × 10–14	3.560 × 10–18

Petrophysical Parameters

The initial fluid saturations, porosity, and absolute permeability used in this work are shown in Table . The producer back pressure is kept constant as shown in Table . The relative permeability curves, which can also be found in Rabiu Ado et al.,[5] are given in Figure . All the three models were run with these parameters.

Table 11

Fluid Saturation, Porosity, and Absolute Permeability as Used in Each Model

initial oil saturation, So	0.85
initial water saturation, Sw	0.15
initial gas saturation, Sg	0.00
porosity	0.34
vertical permeability (mD)	3450
horizontal permeability (mD)	11,500
producer back pressure (kPa)	200

Figure 14

(a) Oil/water and (b) gas/oil relative permeability curves for the Athabasca bitumen.

THAI Kinetics

The thermal cracking kinetics scheme and the corresponding combustion reactions used in model P are taken from Rabiu Ado et al.[5] and are shown in Table . In this model, the combustion of the oil pseudo-component was found to be insignificant. This is because at the high preheating temperature (up to 900°C), there is no oil presence around the preheated zone and oxygen did not bypass the HTO zone. Throughout the combustion period, the combustion zone peak temperature remained over 600°C implying that only HTO reaction is present.

Table 12

Direct Conversion Cracking Kinetics Scheme and Combustion Reactions as Used in Model P[5]a

thermal cracking reactions	frequency factor (min–1)	activation energy (kJ/mol)	heat of reaction (kJ/mol)
IC → 2.0471 MC	3.822 × 1020	239.01	0.00
MC → 0.4885 IC	3.366 × 1018	215.82	0.00
MC → 2.3567 LC	1.132 × 1015	184.88	0.00
LC → 0.4243 MC	1.524× 1015	180.45	0.00
IC → 77.4563 COKE	2.320 × 1015	180.88	0.00
Combustion Reactions
IC + 98.869 O2→77.456 COX + 46.904 H2O	1.812 × 108 kPa–1	138.00	4.00 × 104
MC + 49.069 O2→37.075 COX + 25.953 H2O	1.812 × 109 kPa–1	138.00	1.60 × 104
LC + 32.025 O2→4.600 COX + 35.623 H2O	1.812 × 1010 kPa–1	138.00	4.91 × 102
COKE + 1.22 O2→COX + 0.565 H2O	1.000 × 1010 kPa–1	123.00	3.90 × 102

Reprinted with Permission from Rabiu Ado et al.[5]

Model G is based on the modified Greaves et al.[29] model in which the stoichiometric coefficients and the kinetics parameters of the cracking reaction were adjusted (Table ) in order to obtain better prediction of the fuel availability and produced oxygen concentration. In this model, the oil combustion reactions were also observed to be insignificant for the same reasons given above.

Table 13

Split Conversion Cracking and Combustion Kinetics Based on Modified Greaves et al.[29] Model and as Used in Model Ga

thermal cracking reaction	frequency factor (min–1)	activation energy (kJ/mol)	heat of reaction (kJ/mol)
HEAV → 1.6 LITE + 46.6 Coke	1.50 × 109	99.00	0.00
Combustion Reactions
HEAV +80 O2→ 26.7 H2O + 68.7 COX	1.81 × 1011kPa–1	138.00	5.91 × 102
LITE + 19 O2→ 14.5 H2O + 11.8 COX	1.81 × 1012 kPa–1	138.00	4.91 × 102
CH + 1.2 O2→ 0.5 H2O + COX	8.60 × 107 kPa–1	123.00	3.65 × 102

Reprinted with Permission from Greaves et al.[29]

Model B is based on Belgrave et al.[18] kinetics schemes (Table ). After running this model with the default Belgrave et al.[18]PVT data and kinetics parameters, a significant deviation between the experimental and simulated peak temperature, oil production rate, cumulative oil, and so forth was observed. This is similar to the observation made by Yang and Gates.[13] The LTO reactions were also found to be insignificant despite the presence of oil around the preheated zone. Therefore, to better represent the Athabasca bitumen in question, Aspen HYSYS was used to generate the PVT data for the two oil pseudo-components. Thus, it became necessary to tune the kinetics parameters in order to obtain a better history match. After running the model with slightly tuned parameters, it became obvious that the model is most sensitive to the level of fuel available and because of that to the kinetics parameters of the coke formation reaction (i.e., thermal cracking). With peak temperature reaching up to 900°C during the PIHC, all the oil present around the inlet zone of the sandpack, where the heater was embedded, was displaced by the end of the PIHC. As a consequence, only coke was observed to be present around the HI. This means that only the kinetics parameters of coke formation was tuned because the effect of LTO reactions was negligible.

Table 14

Direct Conversion Cracking, LTO, and HTO Reactions Based on the Belgrave et al.[18] Scheme as Used in Model Ba

thermal cracking reaction	frequency factor (min–1)	activation energy (kJ/mol)	heat of reaction (kJ/mol)
MALT → 0.46 ASPH	5.46 × 1014	234.70	0.00
ASPH → 68.69 COKE	6.96 × 1017	174.20	0.00
ASPH → 20.75 COX	8.17 × 1010	176.30	0.00
LTO Reactions
MALT + 3.43 O2→ 0.58 ASPH	7.69 × 106 kPa–0.43	86.73	1.30 × 103
ASPH + 7.51 O2→ 87.18 COKE	2.49 × 106 kPa–4.76	185.60	2.86 × 103
HTO Reactions
CH + 1.22 O2→0.50 H2O + COX	8.60 × 106 kPa–1	123.00	3.65 × 102

Reprinted with Permission from Belgrave et al.[18]

In each model, the stoichiometric coefficient of CO is obtained using the combination of mole and mass balance. Each oil pseudo-component is assigned an atomic hydrogen to carbon (H/C) ratio. From that, the molecular weight of hydrogen contribution toward the overall molecular weight of the oil pseudo-component was calculated and hence that of the carbon was by default obtained. For the combustion of each of the oil pseudo-component, hydrogen and carbon balances were used to obtain the stoichiometric coefficient of the reactant oxygen. Because the ratio of CO2 over (CO + CO2) was reported from the experiment, it automatically means that the CO2 to CO ratio is known. Additionally, the mol% of CO2 was reported from the experiment, and thus, the mol% of CO is also known. Through oxygen and carbon atom balances, the values of X and the stoichiometry of CO are, respectively, obtained such that each reaction was balanced in terms of mass.

Boundary Conditions

To simulate the combustion cell, no flow boundary condition is assumed all over the cell boundary except via the HI and HP. The HI well is flow-controlled with the air injection rate at the early stage set to 8000 Scm3 min–1 before it was increased by one-third to12,000 Scm3 min–1 at 190 min and maintained until 320 min. A bottom hole pressure of 200 kPa and a total liquid production rate of 25 cm3 min–1 are, respectively, specified as the primary and secondary boundary conditions for the HP well. This allows the simulator to enforce either the pressure or the flow, depending on which one is violated, as the primary constrain. In all the three models, the HP well is assumed to contain no oil as an initial condition. This is because at a typical reservoir temperature of 5–15°C,[45] the Athabasca bitumen is virtually immobile and preheating must be performed in order for the mobility to be established. Also, it is assumed that, during the experiment, the heat loss only occurred from both the overburden and underburden. Therefore, the heat loss parameters were selected via trial and error until all the experimentally measured performance parameters of the THAI process are matched. It should be noted that because there is no material flow in or out of the overburden and underburden, the heat loss only occurred via conduction.

Grid Sensitivity Study

In all the three models, the sensitivity of the simulation results to the grid size was investigated and an optimum number of grid blocks (GBs) was chosen. In this case, 38,000 GBs after refinement were used in each model with the dynamic grid refinement (DynaGrid) option enabled. The DynaGrid option works by derefining the mesh (i.e., restoring the parent GB size) when the temperature difference in the child grids is less than or equal to 30°C and the global mole fraction of any component is less than or equal to 3%. For the result of the sensitivity study, see Rabiu Ado et al.[5]

Table 6

Liquid Specific Heat Capacity Coefficients for Model P

	coefficients
components	A	B	C	D	E
LC	65.23	1.370	–6.404 × 10–4	2.959 × 10–14	–1.738 × 10–17
MC	32.45	3.364	–1.571 × 10–3	–3.072 × 10–14	1.264 × 10–17
IC	–98.16	6.910	–3.226 × 10–3	–1.558 × 10–14	5.003 × 10–18

Table 7

Liquid Thermal Conductivity Coefficients for Model G

	coefficients
components	A	B	C	D	E
LITE oil	2.410 × 10–1	–3.914 × 10–4	3.784 × 10–7	–1.440 × 10–10	–2.389 × 10–13
HEAV oil	1.465 × 10–1	2.129 × 10–4	–9.322 × 10–7	1.110 × 10–9	–4.890 × 10–13

Table 8

Liquid Specific Heat Capacity Coefficients for Model G

	coefficients
Components	A	B	C	D	E
LITE oil	31.09	1.037	–4.843 × 10–4	3.664 × 10–14	–2.287 × 10–17
HEAV oil	–0.7038	6.094	–2.846 × 10–3	2.982 × 10–16	–8.314 × 10–20

Table 9

Liquid Thermal Conductivity Coefficients for Model B

	coefficients
Components	A	B	C	D	E
MALT	2.094 × 10–1	–2.268 × 10–4	5.921 × 10–8	1.414 × 10–10	–1.848 × 10–13
ASPH	1.930 × 10–1	–6.933 × 10–5	–1.964 × 10–7	2.870 × 10–10	–1.420 × 10–13