1. Introduction
It goes without saying that available oil contents are recognized as an unsustainable source of fuel energy in the world, and it will not take a long time for whatever remains in conventional oil reservoirs to run out. As a result, the global demand will be more inclined toward exploiting heavy oil reservoirs, a confident energy source after the conventional one. In contrast to several other reservoir types, such as light shale reserves, where enhanced oil recovery (EOR) applications can be implemented after the primary production stage, Steam-assisted gravity drainage (SAGD) operators tend to start the process after the wells are completed. As has been vastly delineated in the literature, because it has viscosities that are too high, the oil in heavy oil reservoirs cannot move naturally. Thus, researchers have been developing various methods to give companies the ability to practically deal with the problem [
1,
2]. Among the different methods, it appears that thermal methods are more convincing than others. One of the thermal methods, which has a huge reputation in extracting heavy oil, is the SAGD process (or its variations).
Figure 1 shows a simplified demonstration of this process. As can be seen, the reservoir is drilled with at least a pair of horizontal wells for the SAGD process to progress [
3]. Steam is injected from the upper horizontal well (named injector) and the heated oil is summoned up from the well located at the bottom of the reservoir (named producer). The role that steam plays is to transport heat to the viscous, cold areas in the pay zone. There, the conduction mechanism has its foremost effect in heating the untouched oil. While the untouched heavy oil gradually becomes warmer, its viscosity sharply decreases and consequently, thanks to the gravity force, flows down toward the producer [
4].
In general, the simulation of a reservoir in a SAGD process is performed by two methods. The first one originates from numerical methods, which include many equations and sheer conditions [
6]. Almost all commercial software related to petroleum engineering uses this method. It is believed that numerical methods express accurate results in different EOR manners [
6,
7,
8], since they take into account all governing equations, including all transform phenomena expressions and thermodynamic equations of state. However, its most noteworthy flaw that ought to be highlighted is the weaknesses in quickly converging on the ultimate results [
7], which has made researchers seek other methods. That is why petroleum engineers have been developing an alternative to overcome this issue. SAGD analytical models, as a second method, tend to adapt the main in-situ mechanisms and processes that occur during the oil production stage. Researchers have usually made many simplifying assumptions to obtain an equation capable of rapidly predicting the oil production rate in SAGD. In addition to accurate outputs, attempts were made for the developed equation(s) to require fewer reservoir input data to yield the desired outputs [
9].
The first person to model the SAGD process was Butler [
10]. In 1981, according to different observations from some SAGD tests, Butler and Mac Nab [
10] developed the first simple mathematical formula to predict the production rate of the heated oil in their experimental model. Similarly to what is depicted in
Figure 1, they assumed that a steam chamber is formed as superheated steam is injected into the reservoir through the injector. Moreover, they defined a thin layer of oil between the steam chamber and the cold oil region as the interface. As the heated oil flows down toward the producer, the interface grows sideways to reach the walls of the bed. Therefore, Butler and Mac Nab [
10], using Darcy’s law, heat conduction, and also a mass balance in the vicinity of the interface, formulated the sideways expansion of the steam chamber at steady state conditions. It should be mentioned that their outcomes were not very precise in comparison to the experimental data. Again, it was Butler [
11] who came up with a new idea in 1985. To modify his previous work, he changed the steady state conditions that were assumed in the earlier work. As he stated, the major problem in all preceding equations was assuming a steady temperature distribution along the interface. It was quite non-realistic and yielded inaccurate results. Thus, to overcome this shortcoming, he divided the interface into several segments and then, applied the mass and energy balances to each. In his new model, each segment (element) independently moves sideways with its own specific velocity, and the location of interface as well as the oil drainage rate was forecasted more accurately.
Several years later, a linear geometry model was suggested by Reis [
12] to find out more confident results under steady state production rates. Having introduced steam-oil interface as a straight line moving at a velocity
U, Reis combined all governing equations in an element of interface to form a single oil production equation as follows:
which shows the oil drainage rate from one side of the steam-oil interface. Equation (1) has the following parameters:
= oil production rate
= porosity
= initial oil minus residual oil saturation of system
= oil permeability
= reservoir height
= thermal diffusivity of reservoir
= coefficient of velocity
= oil viscosity at steam temperature
= coefficient of viscosity
Nonetheless, Equation (1) shows the independency of the oil production rate from time, which compared with other models, does not look desirable.
Later, it was Pooladi-Darvish [
13] who used Butler’s idea and the HIM (heat integral method) theorem in order to change the energy balance PDE (partial differential equation) over an interface segment into a simple ODE (ordinary differential equation). In this way, Pooladi-Darvish et al. [
13] stated the SAGD process as a moving boundary layer problem and imagined an exponential or polynomial equation for the temperature distribution profile ahead of the interface layer. Eventually, he modified the performance of Butler’s model [
11] in a tremendous way. Similarly to what Pooladi-Darvish et al. [
13] had done in his work, Rabiei et al. [
14] developed a semi-analytical model to estimate the amount of oil drainage in a steam-solvent gravity process. In fact, they involved the solvent and thermal layer in the interface layer in the expanding solvent steam assisted gravity drainage (ES-SAGD) process. HIM was used to find the distribution of the solvent concentration and temperature ahead of the steam/solvent-oil interface. Next, the inverted triangle theory first introduced by Reis [
12] was employed to obtain the oil flow rate. Then, they compared their model with a numerical simulator to prove the accuracy of their model. Finally, it was postulated that the ES-SAGD process had a better potential in oil recovery than the general form of the SAGD process.
Regarding Reis’s linear geometry model, Azad and Chalaturnyk [
14] suggested a circular steam chamber shape for the sideways expansion of the chamber. Like Reis, they defined the flow potential function ∇Φ such that the increase in the chamber’s diameter in each time step influences Darcy’s equation results and the production rate. Aside from that, they considered the oil reservoir as circular slices and could replace the constant relative permeability with varying relative permeability in Darcy’s equation. They assumed three different possible steam chamber locations for their circular model. Consequently, as the first researchers, they could develop a semi-analytical model to simulate the beginning of a SAGD process when the steam chamber was forming. Despite a huge success in production prediction, this model is unable to show the real chamber shape when a SAGD process is in operation. The idea of presenting a model to delineate a conical steam chamber shape is more convincing. This was why Sabeti et al. [
15] proposed an exponential geometry for the interface locations to compensate the obvious problem in Reis’s [
11] and Azad’s [
14] models. More explanation about the exponential geometry approach will be presented in the next section. But concisely, instead of the linear assumption based on Reis’s theory, the interface location was plotted by such the exponential function of
. Finally, Equation (2) was derived as oil production rate related to the SAGD process:
where,
b and
C are dimensionless coefficients and L represents the length of the steam-oil interface layer. Compared with Equation (1), Equation (2) has a second term which implicitly defines the influence of the interface curvature over the amount of the production rate. Compared with the old version of the SAGD models, this exponential geometry had a promising accuracy in both the interface location and the oil production rate. Shi and Okuno [
16] took Reis’s idea and temperature variation along the edge of a steam chamber account to develop an analytical model to estimate oil production rate and steam oil ratio (SOR) in the Surmont SAGD project.
Taking into account all aforementioned studies and following our previous work, the authors of this study decided to move one step further to increase the similarity between the mathematical model and the real mechanism behind the SAGD process. To adapt the exponential geometry model to the real conditions, the assumption of must be obviated from the model. To do so, this paper recommends that the exponential interface, as an analogy for the steam-oil interface, be divided into several segments so that the more appropriate relation, , can be used in the mathematical formula. The present work also uses temperature variation ahead of each segment of the steam-oil interface to find the flowing oil rate.
2. Theory
At the SAGD process, steam is injected from the top well in a reservoir, and the drained oil goes out from the bottom well. Meanwhile, three different scenarios happen respectively. The operation starts by developing a steam chamber, and it may need a thermal circulation process for a while before the producer well can play its own role in the SAGD process. After the steam chamber touches the cap rock, the second stage starts. In this stage, the oil production occurs along with the sideways expansion of the steam chamber. Afterwards, when the steam chamber has completely expanded to the reservoir boundaries, the falling phase begins. In the third stage, the height of oil reservoir begins to gradually descend and at the same time, the oil production rate drops, too. The study described in this work encompasses the second and third stages of a real SAGD process, therefore, it is clear that the equations deployed here are not applicable in the first stage at all. It is worth mentioning that the first stage matters to some researchers, such as Zargar [
17] and Irani [
18], and they presented their own equations.
The flow of the heated oil ahead of the interface layer complies with the Darcy’s law in the porous media. As shown in
Figure 2, for a tiny element ahead of the steam-oil interface, the production rate is formulated as:
in which, the following parameters are included:
= differential oil flow
= distance from interface
= the flow potential function
As mentioned before, unlike the exponential geometry model [
5], the flow potential function is calculated in this way:
where,
θ represents the angle between the steam-oil interface and horizontal line.
To find the oil viscosity ahead of the interface layer, the equation suggested by Butler [
11] is used:
In Equation (5), Butler simply related the oil viscosity at distance
to heat penetration depth,
. Now, combining Equations (3)–(5), the oil production rate is calculated:
In this survey, the steam chamber is depicted as the exponential function of Equation (7) in which C and b are coefficients that determine the curvature of the steam chamber interface during the steam injection:
It was shown that the above exponential function has the ability to locate the interface pertinently with respect to the linear geometry demonstrated by Reis [
11], hence, the authors of this study take advantage of Equation (7) again. Using Equation (7), the amount of
can be defined as:
If Equation (8) is plugged into Equation (6), the final oil production rate with respect to time is obtained:
The only parameter that so far remains unknown is the heat penetration depth. Many studies have been conducted to find
[
11,
13,
14,
20,
21,
22]. As an example, here is an equation derived by Butler [
11], which is used in developing the new flow rate equation presented in this study:
To get Equation (10), an energy balance on the element defined in
Figure 2 is required (see
Appendix A). This differential heat equation asserts that the gradient of penetration depth with respect to time is inversely proportional to its own amount, and that linearly varies with the velocity
perpendicular to the interface. According to Equation (10), it makes sense that the model needs the velocity of each segment to calculate the penetration depth in each time step. By performing a lumped mass balance over the whole area of the reservoir, as is shown in
Figure 3, the maximum frontal velocity
is computed. In this case, the total oil production is equal to the swept area in the reservoir, which covers the same area as the steam chamber inside the reservoir does. Thus,
where
H and
Ws represent the reservoir height and the half-width of the steam chamber, respectively. Equation (11) is differentiated with respect to time as shown in the following line:
Now, by manipulating Equation (12), the required relation for locating the steam-oil interface as well as the penetration depth is obtained:
The maximum velocity, which is believed to occur at the top of the interface, has the following relationship with the velocity perpendicular to it:
Having the maximum velocity (
) available from Equation (14), the location of the top of the interface can be found at each time step:
where
indicates a time interval and n denotes the number of the time intervals which increases with the passing of time. Following the new width of the steam chamber, an interface velocity for each segment can be calculated with the following chain of formulas:
Here, subscript
i indicates each segment number of the interface:
In the previous paper (the exponential geometry model), the heat penetration depth through the interface does not appear in the oil drainage rate equation. But in this work, it clearly does (pay attention to Equation (9)). Please note that, the more accurately the heat penetration depth is used in Equation (9), the more accurately results are obtained through solving the equation.
The procedure for finding the amount of steam required for a SAGD process is the same as that introduced by Sabeti et al. [
5]. The three main sources of energy consumption in SAGD consist of that to keep the steam chamber warm, that to keep the interface heated, and that lost to the overburden. Eventually, after some manipulating, the equation which represents the required steam rate to proceed with the SAGD process is obtained by:
where
and
are the rate of latent heat injection and formation heat capacity of the reservoir, respectively. To get the interface length by using Equation (7), the following function can be defined:
4. Results and Discussion
This section is aimed at validating the new semi-analytical model, and then applying it to a field reservoir under real-world conditions. To do so, one of well-known experimental works recorded by Chung and Butler [
23] was used.
Table 1 represents the physical properties of the oil and reservoir that they considered in their tests. They used Ottawa sand as porous media and saturated it with Cold Lake bitumen at room temperature. The experiment was conducted in two discrete well configurations, A and B. Here, well configuration B, which is in compliance with the present model, was chosen for more investigation. The distance between the producer and injector was set to 20 cm in well configuration B. In SAGD, unlike other EOR processes, the injected fluid does not push the oil back. Thus, before an experiment was started, they applied a steam circulation through the injector for 30 min to form an initial communication path between the two wells.
It must be mentioned that test B does not have the first step (rising of the chamber) in the SAGD process. Therefore, its recorded results are certainly closer to the presented mathematical model in comparison to test A. In
Figure 5, the results of this new model are compared with the recorded data (test B) from Chung and Butler [
23], as well as outputs from the exponential model. It was assumed that the scale coefficient
C in Equation (9) and the heat penetration depth
had the values of 2.25 ×
H and 0.15 ×
H, respectively. The oil production rate, seen in
Figure 5, is approximately constant in both the experimental data and the exponential model from the beginning of the experiment for about 1.5 h, due to the expansion of the steam chamber toward the formation boundaries. After that, the oil production rate follows a falling trend as the reservoir becomes depleted. The oil production rate follows a very similar trend, according to the results of the new model and the experiment. However, at the beginning of the process, the mathematical model predicts a lower production rate, which can be understandably attributed to the first estimation of the heat penetration depth in step 3 of the proposed solution algorithm. Equation (9) includes the three important parameters of heat penetration depth, the angle of interface with respect to horizontal, and the formation height, each of which automatically affects the production rate. Using Equation (10) in the solving procedure, the amount of heat penetration in each segment can be expected to gradually increase. The steady rise of the last segment’s depth on the interface as time goes on is shown in
Figure 6. This increment influences the production rate. However, although the proposed model presents an increase in the depth, as shown in
Figure 7, there is a reduction in the formation height and the interface angle simultaneously. While the penetration depth directly increases the oil production rate, the reductions in formation height as well as the interface angle indirectly lower the production rate. Eventually, these changes lead to a decreasing trend in the production rate (as shown in
Figure 5 after 1.5 h). Please note that Equation (10) takes into account the two extremes (see
Appendix A) and is not exact. Furthermore, using Equation (10), one may notice the unrealistic non-stop rise of the heat penetration depth with time. Thus, taking that into account, the new model cannot give very accurate results in comparison with the experimental data, as illustrated by
Figure 5.
In order to show the importance of the new model, the temperature distribution in the reservoir at various times has been sketched in
Figure 8. Isothermal lines are clear in both the simulation and experimental cases. The contour plots designed by the model fairly agree with Chung and Butler’s [
23] thermal screenshots. It appears that after 1 h, when the production rate reaches its peak, the steam chamber is fully expanded across the reservoir and starts its falling phase. The value of
C in the exponential function (Equation (7)) has the potential to change the curvature of the isothermal lines in
Figure 8. In other words, the more curved the isothermal lines, the lower the value of the scale coefficient.
It is worthwhile to determine if the initial estimation for the value of heat penetration depth (
) has a considerable impact on the oil production rate. Theoretically, it does not hugely affect the results. However, in order to scrutinize its influence, several runs have been taken with different initial estimates. As expected, the higher the initial estimate for the heat penetration depth, the higher the prediction of Equation (10). In the same way, as can be seen in
Figure 9, as the time passes, higher initial estimates create higher values for the heat penetration depth. Apart from that,
Figure 10 represents various outputs coming from these initial estimates. At the first glance, it may be concluded that there are huge discrepancies among different cases at the beginning of the process. However, they gradually demonstrate similar trends after passing their peaks and converge to a production rate of 20 gr/h. As a result, the first estimate does not actually limit the use of this method in other reservoirs with different properties.
As conduction is the main heat transfer mechanism in SAGD, thermal diffusivity is a key parameter in transferring heat to oil. To study the effect of thermal diffusivity on oil production in SAGD, three different values of thermal diffusivity were selected in the model under study. The values are half of, equal to, and twice the value of those mentioned in
Table 1. The resulting oil production rates in the three cases during SAGD production are shown in
Figure 11. As can been seen, the higher thermal diffusivity, the higher oil the production rate, especially at the beginning of the process. Conductive heat transfer occurs at a higher rate across the reservoirs with a higher thermal diffusivity, which in turn raises the temperature ahead of the steam interface. In fact, in a reservoir with a higher thermal diffusivity, the heat transfer is enhanced and higher oil production rates occur.
5. Model Verification Using the UTF Project Phase B
Aside from validating the new model with the pilot test, one might want to know if it is accurate enough to estimate both the oil production rate and the amount of steam required in a real SAGD project. To find out, the properties of the Athabasca reservoir in phase B of the UTF project were used. A summary of reservoir properties is listed in
Table 2. In 1993, phase B of the UTF project was initiated by drilling three pair wells, all of which had a length of 500 m and were located 70 m apart from each other at a depth of 22 m [
24,
25]. The process continued until 1998, when its investors decided to modify the process by adding solvents into the injectors [
24,
26].
Introducing all the data in
Table 2 to the solution algorithm designed in the previous sections produced diagrams for both oil production and steam injection rates, as demonstrated in
Figure 12. Two important facts must be remembered about these diagrams. When the UTF project was started in 1993, the steam chamber began to form its usual conical shape as the first stage of the SAGD process. This lasted almost a year, and during this period, the oil production rate from each well kept growing before reaching a plateau. Firstly, the oil production diagram in
Figure 12 does not take into account the production data before 1994. Secondly, when the new model was employed using the Athabasca data, it was observed that the value of heat penetration depth encroached its real size. Therefore, a maximum value was set as a constraint to prevent its excessive growth. To do so, the maximum thickness
of the interface was assumed to be an unknown adjustable parameter. That value was determined by introducing an objective function, such as the time-weighted root-mean-squared relative error defined by Li et al. [
27] (see
Appendix B). By optimizing that objective function, the value of the maximum depth
was calculated to be about 2 m. Despite these two stipulations,
Figure 12 asserts the potential of the proposed model in predicting SAGD processes in homogenous formations. Not only does the oil production rate show good compatibility with the field data, but also, the amount of steam per day follows a logical and acceptable estimation.