An Improved Mathematical Model for Accurate Prediction of the Heavy Oil Production Rate during the SAGD Process

Steam-assisted gravity drainage (SAGD) is one of the most successful thermal enhanced oil recovery (EOR) methods for cold viscose oils. Several analytical and semi-analytical models have been theorized, yet the process needs more studies to be conducted to improve quick production rate predictions. Following the exponential geometry theory developed for finding the oil production rate, an upgraded predictive model is presented in this study. Unlike the exponential model, the current model divides the steam-oil interface into several segments, and then the heat and mass balances are applied to each of the segments. By manipulating the basic equations, the required formulas for estimating the oil drainage rate, location of interface, heat penetration depth of steam ahead of the interface, and the steam required for the operation are obtained theoretically. The output of the proposed theory, afterwards, is validated with experimental data, and then finalized with data from the real SAGD process in phase B of the underground test facility (UTF) project. According to the results, the model with a suitable heat penetration depth correlation can produce fairly accurate outputs, so the idea of using this model in field operations is convincing.


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. 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: q o = oil production rate ϕ = porosity ∆S o = initial oil minus residual oil saturation of system k o = oil permeability H = reservoir height α = thermal diffusivity of reservoir a = coefficient of velocity υ os = oil viscosity at steam temperature m = 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 f (x) = C(1 − e −bx ). 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 ∇Φ = g H L 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, ∇Φ = g sin θ, 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.

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: dq = differential oil flow dξ = 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. 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), (4) and (5), the oil production rate is calculated: (6) 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  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 sin θ can be defined as: If Equation (8) is plugged into Equation (6), the final oil production rate with respect to time is obtained: Processes 2020, 8, 235 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 U v 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 U m 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 W s 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: 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: 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 (U m ) available from Equation (14), the location of the top of the interface can be found at each time step: where ∆t 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 Q inj and M R 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:

Computation Procedure
In order to find the production rate, interface location, interface velocity, temperature distribution ahead of interface, and the amount of steam required to complete a SAGD process, the following sequence needs to be implemented. Nonetheless, before that, the height of the reservoir should be divided into some evenly distributed segments, as displayed in Figure 4.
(1) Using Equation (2), the amount of oil production rate is calculated as a first supposition.
(2) The value of the maximum horizontal velocity is found by Equation (13).
(3) Equations (15)- (17) are used to find the location of the interface. Furthermore, the magnitude of horizontal velocity at each segment is obtained by Equation (18). (4) Having U v and γ and applying Equation (10), the value of the heat penetration depth γ for the next time step is calculated. The first value of the heat penetration depth γ at time t = 0 for each segment of the steam-oil interface needs to be estimated. (5) To get the value of the oil production rate, Equation (9), with its required parameters (formerly defined in the previous steps), is used. (6) The calculated γ should not have a significant discrepancy with its previous value. To fulfill this, the new γ should be recalculated with the new values from step 2. (7) In the end, regarding an accurate oil production rate, Equations (19) and (20) are employed to get the energy required for the SAGD process.
Repeat steps 2 through 7 until the final time step is reached.
respectively. To get the interface length by using Equation (7), the following function can be defined:

Computation Procedure
In order to find the production rate, interface location, interface velocity, temperature distribution ahead of interface, and the amount of steam required to complete a SAGD process, the following sequence needs to be implemented. Nonetheless, before that, the height of the reservoir should be divided into some evenly distributed segments, as displayed in Figure 4. (2), the amount of oil production rate is calculated as a first supposition.

1) Using Equation
2) The value of the maximum horizontal velocity is found by Equation (13).
3) Equations (15)- (17) are used to find the location of the interface. Furthermore, the magnitude of horizontal velocity at each segment is obtained by Equation (18).

4) Having
Uv and  and applying Equation (10), the value of the heat penetration depth  for the next time step is calculated. The first value of the heat penetration depth  at time t = 0 for each segment of the steam-oil interface needs to be estimated. 5) To get the value of the oil production rate, Equation (9), with its required parameters (formerly defined in the previous steps), is used.
6) The calculated  should not have a significant discrepancy with its previous value. To fulfill this, the new  should be recalculated with the new values from step 2.

7)
In the end, regarding an accurate oil production rate, Equations (19) and (20) are employed to get the energy required for the SAGD process.
Repeat steps 2 through 7 until the final time step is reached.

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. 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.  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. 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. 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.  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. 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. Figure 11. Study on the effects on thermal diffusivity variation on oil production rate in the present SAGD model. 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. 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. Figure 11. Study on the effects on thermal diffusivity variation on oil production rate in the present SAGD model. Figure 11. Study on the effects on thermal diffusivity variation on oil production rate in the present SAGD model.

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 γ max 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 γ max 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. Processes 2020, 8, x FOR PEER REVIEW 15 of 19 Figure 12. The oil production and steam injection rates from the proposed model, as well as field data, and phase B of the UTF project.

Conclusions
In this study, an effort was made to introduce a fast, reliable, and valid model for the SAGD process in order to accurately predict the oil production and steam injection rates, as well as the location of the interface at different times. The new model was derived by holding onto most of the assumptions in the latest version of the SAGD model (exponential geometry). It was indicated that the model needed the steam-oil interface to be divided into equal segments for its semi-analytical approach to obtain an equation for the oil production rate. The outcomes of the model were then evaluated by the experimental data from Chung and Butler [23], and results showed an acceptable coincidence between them. It was further found that the rate of oil production is affected by three main parameters: heat penetration depth ahead of the interface, the angle of interface with respect to the horizontal line, and height of the formation. When the heat penetration depth gradually increases, the formation height, as well as the interface angle, declines with time. Thus, their combined effect makes a downward trend in the diagram of the heavy oil production rate versus time. Moreover, having a reliable correlation for heat penetration depth was proven to be vital. Hence, an erroneous estimate of this would negatively affect the accuracy of the oil production rate, steam injection rate, and the interface location in the SAGD process. Finally, the proposed model was run using the data from phase B of the UTF project. The outputs matched fairly well with the field data.

Conclusions
In this study, an effort was made to introduce a fast, reliable, and valid model for the SAGD process in order to accurately predict the oil production and steam injection rates, as well as the location of the interface at different times. The new model was derived by holding onto most of the assumptions in the latest version of the SAGD model (exponential geometry). It was indicated that the model needed the steam-oil interface to be divided into equal segments for its semi-analytical approach to obtain an equation for the oil production rate. The outcomes of the model were then evaluated by the experimental data from Chung and Butler [23], and results showed an acceptable coincidence between them. It was further found that the rate of oil production is affected by three main parameters: heat penetration depth ahead of the interface, the angle of interface with respect to the horizontal line, and height of the formation. When the heat penetration depth gradually increases, the formation height, as well as the interface angle, declines with time. Thus, their combined effect makes a downward trend in the diagram of the heavy oil production rate versus time. Moreover, having a reliable correlation for heat penetration depth was proven to be vital. Hence, an erroneous estimate of this would negatively affect the accuracy of the oil production rate, steam injection rate, and the interface location in the SAGD process. Finally, the proposed model was run using the data from phase B of the UTF project. The outputs matched fairly well with the field data.   Taking a lump heat balance over the element shown in Figure 2.
where, q r and K h are the amounts of heat stored ahead of the interface and thermal conductivity of the reservoir, respectively. Equation (A1) consists of three terms: on the left hand side of the equation it denotes the rate of heat accumulated ahead of the interface; and on the right hand side it represents the heat of the steam transported with a forward flux as well as the heat left behind in the reservoir because of the forward motion of the interface.
Also taking ϕS o ρC p (T s − T r ) as a constant value with time, it is clear that, Besides, the dimensionless temperature is defined by, Then, the temperature gradient on the steam-oil interface where ζ = 0 may be defined as follows: Combining Equations (A1), (A2) and (A4), a very simple heat penetration depth equation is formed.
Similar to what Butler mentioned in his work, the temperature gradient term does not have a constant value in Equation (A5) and depends on the varying heat penetration depth and the advancing interface velocity U v . Equation (A5) may be solved by considering two extreme conditions for ∂T * ∂ζ ζ=0 .
(1) Under the steady state condition, the value of heat penetration depth remains constant, so dγ dt = 0 and then, (2) Under a stationary case when U v equals zero, T * | ζ=0 = 1 and T * | γ,t=0 = 0. The dimensionless temperate has the following relation, The corresponding temperature gradient on the interface layer is then given by: Equation (A6) from the first case and Equation (A8) from the latter have similar general forms. In both, the temperature gradient is inversely proportional to the heat penetration. Therefore, Butler assumed that the real equation for heat penetration depth lies between these two cases, so this premise leads to the following equation: Now, substituting Equation (A9) into Equation (A5), a simplified relation for the heat penetration depth used in this paper is obtained.
If ϕS o equals unity [28], the final expression is obtained.