Numerical Simulation of Non-Isothermal Two-Phase Flow in Oil Reservoirs: Slab Geometry and Heated Fluid Injection

1. Introduction

We can mention energy generation, the environment, health, and biological systems where porous media flow knowledge is significant. Specifically, we have as examples oil and gas exploitation, fuel cells, batteries, solid mechanics, and cardiovascular systems [1]. For example, numerical simulation can be used in the oil industry, which invests significant resources to reduce its costs and risks inherent to production.

We can obtain data about flow patterns by performing experiments using reservoir samples. On the other hand, information can also be collected through pressure tests in injection/producing wells and temperature and pressure data analysis. However, these are expensive and may be unfeasible when we want to study several production scenarios. For this purpose, numerical simulation becomes a valuable alternative tool [2]. In this context, we propose here to use it to solve non-isothermal two-phase flow using a one-equation model derived without assuming the hypothesis of local thermal equilibrium.

Numerical reservoir simulation requires a prior understanding of several areas of knowledge, such as mathematics, physics, numerical methods, and advanced programming [3], and it is widely used in the oil industry [2]. It is due to the inherent difficulties in accurately obtaining reservoir properties throughout its extension and the complexity of determining flow behavior in situ. Therefore, it is an ally when we need to maximize hydrocarbon recovery.

We know that the hydrocarbon production stages can be divided into primary, secondary, and tertiary [4]. In the primary, the oil or gas is obtained due to the difference in pressure between the interior of the reservoir and the surface. In the second, part of the large quantity of hydrocarbons remaining in the reservoir is recovered via the injection of an immiscible fluid, available in abundance and at low cost. Finally, we have the chemical, miscible displacement, biological, and thermal methods [5]. All aim to promote higher mobility of hydrocarbons, decreasing the resistance to their flow and promoting increased production [4].

In the case of conventional oils (API gravity greater than 30), production costs are comparatively lower than those associated with heavier oils [6]. However, due to the current depletion of their reserves, special attention has been given to exploring those with API gravity between 5 and 22. On the other hand, as they have a higher viscosity, there is higher resistance to their flow, and production costs are higher [7]. Under these circumstances, the dominant recovery techniques are thermal, non-thermal, or cold and solvent-based [8].

Therefore, an attractive alternative to increase the recovery of heavy oils is the use of so-called thermal methods [9]. For example, static heaters can be positioned inside the reservoir. As a result, there will be a decrease in its viscosity and a consequent increase in its fluidity [10].

As with other recovery techniques, we can also use numerical simulation for thermal methods [11]. The most commonly used is steam injection [12], which has the disadvantage of requiring a large amount of available water in addition to the energy cost associated with steam generation. In addition, it causes harmful effects on the environment due to the means generally used to generate the necessary energy, causing the emission of greenhouse gases [12]. An alternative, in situ combustion, has been applied with relative success in Canada [8].

We can also use electric energy to heat the reservoir by introducing heating elements that transfer heat to the fluid and the rock [7,13]. Another possibility is electromagnetic heating, enabling localized heating and being able to transfer up to 92% of the energy generated to its surroundings [14]. The number and location of the heating elements can vary to maximize the recovery process.

In the physical-mathematical modeling of porous media flow, compared to isothermal flow, we must consider the energy equation in non-isothermal flow in addition to the continuity and momentum equations. In the literature, one- or two-equation models are commonly used, with the second being the most suitable when there are significant differences in the order of magnitude of the thermal properties of the fluids and the porous matrix [15,16]. The one-equation model for single-phase flow [11,17,18] can be derived with or without assuming the hypothesis of local thermodynamic equilibrium [19]. In the case of the two-equation model [15,16], we use separate equations to obtain the temperature of the fluid and the rock, and the heat transfer between them is accounted for via a source term [20].

In the most general case, we must consider a multiphase flow due to the presence of more than one phase when using advanced recovery techniques. For example, Siavashi et al. [21] simulated the non-isothermal two-phase flow of immiscible fluids in a saturated medium, adopting the hypothesis of local thermal equilibrium, and they simulated the injection of heated water via the streamline method. On the other hand, also assuming local thermal equilibrium, Roy et al. [22] studied the non-isothermal, compositional two-phase flow, considering the presence of chemical reactions. We find other applications in Lindner et al. [23] and Kumar et al. [24]. The first dealt with cooling by phase change and used a mixing model, considering a multi-component fluid, without considering local thermal equilibrium. The last one addressed the drying process in food via microwaves, in which a fluid flow model in a porous medium with phase change was used.

Although phase changes may occur due to the temperature increase inside the reservoir, we did not consider them here. As a preventive measure, we monitored the maximum temperature values to prevent them from occurring. Therefore, the case of a three-phase flow is outside the scope. In this initial work, we will only deal with the problem of injecting a heated fluid into a reservoir with slab geometry. Later, we will attack the heating problem using static heating wells.

2. Non-Isothermal Flow in Porous Media

We present below the governing equations that will provide the saturation and average temperature of the reservoir. In addition to establishing how the properties of the fluids and the porous matrix will be calculated.

We considered the following hypotheses in the model used:

• the rock formation is anisotropic in terms of absolute permeability;
• the rock compressibility is small and constant;
• the fluids are Newtonian;
• no chemical reactions occur;
• the flow is laminar at a very low velocity;
• the fluids have a steady chemical composition;
• the flow is non-isothermal, two-phase and three-dimensional;
• the fluids are slightly compressible;
• the thermal conductivities of the rock and fluids are constant;
• and the hypothesis of local thermal equilibrium is not considered.

For two-phase flow in porous media, we we obtain from the mass balance [3]

$${\displaystyle \frac{\partial }{\partial t}}\left(\varphi {\rho }_{\alpha }{S}_{\alpha }\right)+\nabla .\left({\rho }_{\alpha }{\mathbf{v}}_{\alpha }\right)-q_{m\alpha }=0,$$

(1)

where $\varphi$ is the porosity, $\rho$ is the density, $\mathbf{v}$ is the superficial velocity, S is the saturation, $q_m$ is the source term, and the subscript $\alpha$ indicates the phase in question. Furthermore,

$$\varphi ={\varphi }^0\left[1+c_{\varphi }\left(p_o-p^0\right)-c_{\varphi T}\left(T-T^0\right)\right],$$

(2)

where ${\varphi }^0$ is the reference porosity, determined at the reference pressure ($p^0$) and temperature ($T^0$), $c_{\varphi }$ is the rock compressibility coefficient, $c_{\varphi T}$ represents the rock thermal expansivity coefficient, $p_o$ is the pressure of the non-wetting phase (oil) and T is the mean reservoir temperature (to be defined later) [3].

Regarding density [3], we employ

$${\rho }_{\alpha }={\rho }_{\alpha }^0\left[1+c_{\alpha }\left(p_{\alpha }-p^0\right)-c_{\alpha T}\left(T-T^0\right)\right],$$

(3)

where ${\rho }_{\alpha }^0$ is the density of the phase $\alpha$ under reference conditions, $c_{\alpha }$ is the fluid compressibility coefficient, and $c_{\alpha T}$ is the thermal expansion coefficient of the phase $\alpha$.

Although viscosity varies depending on pressure and temperature, we will only consider its variation with temperature [11],

$${\mu }_n=a_{n}exp\left({\displaystyle \frac{b_n}{T-T_{ref}}}\right),$$

(4)

where $a_n$ and $b_n$ are constants determined for a given type of oil, $T_{ref}$ is the reference temperature for calculating viscosity, and the water viscosity is considered constant.

In two-phase flow, Darcy’s classical law must be modified by introducing the relative permeability $k_{r\alpha }$ [3],

$${\mathbf{v}}_{\alpha }=-{\displaystyle \frac{k_{r\alpha }}{{\mu }_{\alpha }}}\mathbf{k}\left(\nabla p_{\alpha }-{\gamma }_{\alpha }\nabla Z\right),$$

(5)

where $k_{r\alpha }$ is the relative permeability of $\alpha$ phase, ${\mu }_{\alpha }$ is the viscosity, $\mathbf{k}$ is the absolute permeability tensor (diagonal), $p_{\alpha }$ is the pressure, ${\gamma }_{\alpha }={\rho }_{\alpha }g$, where g is the magnitude of the acceleration due to gravity, and Z is the depth.

Since this is a non-isothermal flow, we use a model that takes into account the variation of relative permeability as a function of temperature [25]:

$$k_{rw}=exp\left[-0.08832T+15.38ln\left(T\right)-64.87\right]{\left[{\displaystyle \frac{S_w-2.673·{10}^{-4}T-0.0382}{1-2.673·{10}^{-4}T+0.1036ln\left(T\right)-0.7623}}\right]}^{A_1}$$

(6)

for the wetting phase (water), where $A_1=-5,578·{10}^{-3}T+2,995$ and, for the non-wetting phase,

$$k_{rn}=\left[0.4947ln\left(T\right)-1.6666\right]{\left[1-{\displaystyle \frac{S_w-2.673·{10}^{-4}T-0.0382}{1-2.673·{10}^{-4}\left(T\right)+0.1036ln\left(T\right)-0.7623}}\right]}^{A_2},$$

(7)

where $A_2=-9.885·{10}^{-4}T+1.349$, $S_w$ is the saturation of the wetting phase, and the temperature T should be given in ^oC. These correlations are valid for temperatures ranging from 70 to 220 degrees Celsius.

Combining Equations (1) and (5), we have, for each phase,

$$\nabla ·\left[{\displaystyle \frac{k_{ro}}{{\mu }_o{B}_o}}\mathbf{k}\left(\nabla p_o-{\gamma }_o\nabla Z\right)\right]={\displaystyle \frac{\partial }{\partial t}}\left[{\displaystyle \frac{\varphi S_o}{B_o}}\right]-q_{s{c}_o},$$

(8)

$$\nabla ·\left[{\displaystyle \frac{k_{rw}}{{\mu }_w{B}_w}}\mathbf{k}\left(\nabla p_w-{\gamma }_w\nabla Z\right)\right]={\displaystyle \frac{\partial }{\partial t}}\left[{\displaystyle \frac{\varphi S_w}{B_w}}\right]-q_{s{c}_w},$$

(9)

where $B_{\alpha }={\rho }_{sc\alpha }/{\rho }_{\alpha }$ is the formation-volume factor, $q_{s{c}_w}=q_{mw}/{\rho }_{s{c}_w}$ and $q_{s{c}_o}=q_{mo}/{\rho }_{s{c}_o}$, being ${\rho }_{s{c}_w}$ and ${\rho }_{s{c}_o}$ the densities at standard (constant) conditions of the wetting phase (water), and of the non-wetting phase (oil), respectively.

The reservoir is saturated

$$S_w+S_o=1$$

(10)

and the effects of capillary pressure,

$$p_c=p_o-p_w,$$

(11)

are accounted for [26]

$$p_c\left(S_w\right)=p_{c_{max}}{\left({\displaystyle \frac{1-S_w-S_{ro}}{1-S_{iw}-S_{ro}}}\right)}^{epc},$$

(12)

where $p_{c_{max}}$ is the maximum capillary pressure, the exponent $epc$ must be provided, $S_{iw}$ is the irreducible water saturation, and $S_{ro}$ is the residual oil saturation.

Therefore, by substituting Equations (10) and (11) into Equations (8) and (9) we obtain the two governing equations:

$$\nabla ·\left[{\displaystyle \frac{k_{ro}}{{\mu }_o{B}_o}}\mathbf{k}\left(\nabla p_o-{\gamma }_o\nabla Z\right)\right]={\displaystyle \frac{\partial }{\partial t}}\left[{\displaystyle \frac{\varphi \left(1-S_w\right)}{B_o}}\right]-q_{s{c}_o},$$

(13)

$$\nabla ·\left[{\displaystyle \frac{k_{rw}}{{\mu }_w{B}_w}}\mathbf{k}\left(\nabla p_o-\nabla p_c-{\gamma }_w\nabla Z\right)\right]={\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{\varphi S_w}{B_w}}\right)-q_{s{c}_w}.$$

(14)

As previously stated, it is common to use the assumption of local thermal equilibrium (equality of mean temperatures of fluids and porous matrix) [22,27]. However, in some situations, it may be inappropriate [28]. As alternatives, we have the one-equation model without adopting such an assumption [19], or the two-equation model [15] or three-equations, with phase change [28] or without phase change [29]. Here, we employ the one-equation model (without local thermal equilibrium), whose dependent variable is the mean reservoir temperature (weighted average of the phase and rock temperatures), and the thermal heating is modeled via source terms [11].

For constant heat capacities, the macroscopic equation governing heat transfer is [30]:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial t}}\left[\left(\rho c_p\right)T\right]& -\nabla ·\left[{\rho }_o{h}_o{\displaystyle \frac{k_{ro}}{{\mu }_o}}\mathbf{k}\left(\nabla p_o-{\gamma }_o\nabla Z\right)\right]-\nabla ·\left[{\rho }_w{h}_w{\displaystyle \frac{k_{rw}}{{\mu }_w}}\mathbf{k}\left(\nabla p_o-\nabla p_c-{\gamma }_w{\nabla }_{w}Z\right)\right]\hfill \\ \hfill & =\nabla ·(\mathbf{K}\nabla T)+{\displaystyle \frac{q_H}{V_b}}+{\displaystyle \frac{{\rho }_o{h}_o{q}_{s{c}_o}}{V_b}}+{\displaystyle \frac{{\rho }_w{h}_w{q}_{s{c}_w}}{V_b}},\hfill \end{array}$$

(15)

where $\mathbf{K}$ represents the effective thermal dispersion tensor, $h_o$ and $h_w$ are the enthalpies of oil and water, $q_H$ represents a volumetric thermal source (when there are heating wells), $V_b$ is the total volume, and $q_{s{c}_o}$ and $q_{s{c}_w}$ represent the volumetric mass flow rates when there are injection/producing wells.

The average temperature of the porous medium, considering both phases and the rock, is given by [30]:

$$\left(\rho c_p\right)T=\varphi \left(1-S_w\right){\left(\rho c_p\right)}_o{T}_o+\varphi S_w{\left(\rho c_p\right)}_w{T}_w+\left(1-\varphi \right){\left(\rho c_p\right)}_r{T}_r,$$

(16)

where ${c_p}_o$, ${c_p}_w$, $T_o$, and $T_w$ are respectively the specific heat and temperatures of the oil and water, and T represents the average reservoir temperature. Furthermore [30],

$$\left(\rho c_p\right)=\varphi \left(1-S_w\right){\left(\rho c_p\right)}_o+\varphi S_w{\left(\rho c_p\right)}_w+\left(1-\varphi \right){\left(\rho c_p\right)}_r$$

(17)

represents the average heat capacity.

The enthalpies of the fluids are calculated knowing the average temperature of the reservoir [11]:

$$h_o={c_p}_{o}T,h_w={c_p}_{w}T.$$

(18)

For example, after applying the Method of Volume Averaging [31] to the energy balance equation, we see the appearance of the effective thermal dispersion tensor in the macroscopic equation describing energy transfer at the laboratory scale [19,32]. It contains the contribution of three terms related to thermal conduction, tortuosity, and hydrodynamic dispersion and, in the case of two-phase flow,

$$\mathbf{K}=\left[\varphi \left(1-S_w\right){\kappa }_o+\varphi S_w{\kappa }_w+\left(1-\varphi \right){\kappa }_r\right]\left(\mathbf{I}+\mathbf{\tau }\right)+\mathbf{D},$$

(19)

where ${\kappa }_o$, ${\kappa }_w$ and ${\kappa }_r$ are the thermal conductivities of oil, water and rock, respectively, $\mathbf{I}$ is the identity tensor, $\mathbf{\tau }$ is the tortuosity tensor, and $\mathbf{D}$ is the hydrodynamic dispersion tensor.

We assume that the tortuosity tensor is diagonal and that its components are constant and given by [33]:

$$\tau =1.23{\displaystyle \frac{{\left(1-\varphi \right)}^{4/3}}{\varphi }}.$$

(20)

We take the hydrodynamic dispersion tensor to be diagonal and anisotropic, with its axial $D_{ax}=D_{xx}$ and transversal $D_{tr}=D_{yy}=D_{zz}$ components determined via [32]:

$$D_{ax}=\varphi \left(1-S_w\right){\kappa }_{n}0.7P{e}_o^{1,1}+\varphi S_w{\kappa }_{w}0.7P{e}_w^{1,1},$$

(21)

$$D_{tr}=\varphi \left(1-S_w\right){\kappa }_{n}0.2P{e}_o+\varphi S_w{\kappa }_{w}0.2P{e}_w,$$

(22)

where $P{e}_{\alpha }=R{e}_{\alpha }P{r}_{\alpha }$ is the Péclet number. The Reynolds ($Re$) and Prandtl ($Pr$) numbers are defined for the mobile phases (water and oil) by $R{e}_{\alpha }={\rho }_{\alpha }{v}_{\alpha }{d}_p/{\mu }_{\alpha }$ and $Pr=c_p{\mu }_{\alpha }/{\kappa }_{\alpha }$, where $d_p$ represents the mean particle diameter of the porous medium.

2.1. Initial and Boundary Conditions

We still need to provide the initial and boundary conditions to complete the physical-mathematical model, whose governing equations are (13), (14), and (15).

We provide pressure, temperature, and saturation as initial conditions at $t=t_0$. For pressure, their values consider the hydrostatic gradients [3]. Similarly, for temperature, geothermal energy causes it to vary at a rate $|\nabla G|$ of 25-30 K/km for the region of the Earth’s crust [34].

Regarding boundary conditions, Dirichlet, Neumann, or Robin-type conditions can be employed on the entire boundary $\partial \Omega$ or on part of it [4].

The reservoir has a slab geometry, and the main flow is in the direction of the x-axis. At the inlet boundary, we impose a pressure gradient and the temperature and saturation of the injected fluid. At the right boundary, we establish a pressure value and zero derivatives for the saturation and temperature in the flow direction. At the other boundaries, we consider zero flow conditions and zero derivatives for saturation and temperature in the directions perpendicular to the reservoir face.

3. Numerical Resolution Methodology

In solving governing differential equations, we employ the Finite Volume Method [2,35]. It is a well-known method having conservative properties that make it a good option when solving systems of partial differential equations [35]. As usual, the mean values of the dependent variables are determined at the center of the finite volumes and do not vary inside them (for a given time).

The domain of dimensions $L_x$, $L_y$, and $L_z$ is discretized using $n_x$, $n_y$, and $n_z$ finite volumes (whose dimensions, not necessarily constant, are $\Delta x$, $\Delta y$, and $\Delta z$) in the x-, y-, and z-directions, Figure 1, so that,

$$\sum _{i=1}^{n_x}\Delta x_i=L_x,\sum _{j=1}^{n_y}\Delta y_j=L_y,\sum _{k=1}^{n_z}\Delta z_k=L_z,$$

(23)

where the indices i, j and k indicate the blocks (cells or finite volumes) in the x-, y- and z-directions, respectively. In this system, we reference the boundaries (faces) of the blocks via the notation ($i\pm 1/2,j,k$), ($i,j\pm 1/2,k$), and ($i,j,k\pm 1/2$) [35].

In Figure 2, we indicate the faces of the finite volumes by $(i-1/2,j,k)=w$, $(i+1/2,j,k)=e$, $(i,j-1/2,k)=n$, $(i,j+1/2,k)=s$, $(i,j,k-1/2)=b$ and $(i,j,k+1/2)=a$. On the other hand, we use capital letters when we refer to the nodes positioned at their respective centers: $(i-1,j,k)=W$, $(i-1,j,k)=E$, $(i,j-1,k)=N$, $(i,j+1,k)=S$, $(i,j,k-1)=A$, and $(i,j,k+1)=B$.

Once we have generated the computational mesh, we perform space-time integration of Equations (13), (14), and (15) in each cell of the domain:

$$\begin{array}{cc}\hfill {\int }_t^{t+\Delta t}{\int }_V{\displaystyle \frac{\partial }{\partial t}}\left[{\displaystyle \frac{\varphi \left(1-S_w\right)}{B_o}}\right]dVdt& ={\int }_t^{t+\Delta t}{\int }_V\nabla ·\left[{\displaystyle \frac{k_{ro}}{{\mu }_o{B}_o}}\mathbf{k}\left(\nabla p_o-{\gamma }_o\nabla Z\right)\right]dVdt\hfill \\ \hfill & +{\int }_t^{t+\Delta t}{\int }_V{q}_{s{c}_o}dVdt,\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\int }_t^{t+\Delta t}{\int }_V{\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{\varphi S_w}{B_w}}\right)dVdt& ={\int }_t^{t+\Delta t}{\int }_V\nabla ·\left[{\displaystyle \frac{k_{rw}}{{\mu }_w{B}_w}}\mathbf{k}\left(\nabla p_o-\nabla p_c-{\gamma }_w\nabla Z\right)\right]dVdt\hfill \\ \hfill & +{\int }_t^{t+\Delta t}{\int }_V{q}_{s{c}_w}dVdt\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill {\int }_t^{t+\Delta t}{\int }_V{\displaystyle \frac{\partial }{\partial t}}\left[\left(\rho c_p\right)T\right]dVdt& ={\int }_t^{t+\Delta t}{\int }_V\nabla ·\left[{\rho }_o{h}_o{\displaystyle \frac{k_{ro}}{{\mu }_o}}\mathbf{k}\left(\nabla p_o-{\gamma }_o\nabla Z\right)\right]dVdt\hfill \\ \hfill & +{\int }_t^{t+\Delta t}{\int }_V\nabla ·\left[{\rho }_w{h}_w{\displaystyle \frac{k_{rw}}{{\mu }_w}}\mathbf{k}\left(\nabla p_o-\nabla p_c-{\gamma }_w\nabla Z\right)\right]dVdt\hfill \\ \hfill & +{\int }_t^{t+\Delta t}{\int }_V\nabla ·(\mathbf{K}\nabla T)dVdt+{\int }_t^{t+\Delta t}{\int }_V{\displaystyle \frac{q_H}{V_b}}dVdt\hfill \\ \hfill & +{\int }_t^{t+\Delta t}{\int }_V\left({\displaystyle \frac{{\rho }_o{h}_o{q}_{s{c}_o}}{V_b}}+{\displaystyle \frac{{\rho }_w{h}_w{q}_{s{c}_w}}{V_b}}\right)dVdt,\hfill \end{array}$$

(26)

where $\Delta t$ is the time increment and the subscript V indicates that the integration domain is the control volume.

Since we chose a widely used method, we will omit the details concerning space-time integration of the governing equations [30]. We opted for a fully time implicit formulation and used space-centered difference approximations.

In the discretized equations, we introduced a new variable, $\mathbb{T}$, which represents the transmissibility evaluated on the faces of the finite volumes

$${{\mathbb{T}}_{\alpha }}_e={\left({\displaystyle \frac{k_x{A}_x{k}_{r\alpha }}{{\mu }_{\alpha }{B}_{\alpha }\Delta x}}\right)}_e,{{\mathbb{T}}_{\alpha }}_w={\left({\displaystyle \frac{k_x{A}_x{k}_{r\alpha }}{{\mu }_{\alpha }{B}_{\alpha }\Delta x}}\right)}_w,$$

(27)

where the subscripts e and w refer to the values that must be evaluated on the $east$ and $west$ faces, $\Delta x_e$ and $\Delta x_w$ represent the distances between the centers of the cells neighboring the finite volume considered, and $A_x=\Delta x_P\Delta y_P$ is the area of the face of the finite volume. Analogous expressions are obtained for the other faces in the y- and z-directions.

Another point of fundamental importance is the approximation of transmissibilities at the interfaces of the control volumes. These terms can be highly nonlinear due to their dependence on the unknowns of the problem. This fact can generate problems with convergence, as reported in [36], and the stability of the numerical method. Knowing this context and remembering that we define the properties at the centroid of the finite volume, we must introduce appropriate approximations to overcome these issues. First, we calculate $\mathbb{T}$ from the product between a geometric term, G, and two dependent on the dependent variables, $F_{(p,T)}{F}_S$. Then,

$$\mathbb{T}=G{F}_{(p,T)}{F}_S.$$

(28)

For example, the geometric terms appearing in the three balance equations are approximated on the east face e by

$$G_{xe}={\left({\displaystyle \frac{k_x{A}_x}{\Delta x}}\right)}_e={\displaystyle \frac{k_{xe}{A}_{xe}}{\Delta x_e}},$$

(29)

where the permeability, in turn, is determined via a harmonic mean

$$k_{xe}={\displaystyle \frac{k_{xP}{k}_{xE}\left(\Delta x_P+\Delta x_E\right)}{k_{xP}\Delta x_E+k_{xE}\Delta x_P}},$$

(30)

and the same type of approximation is used for the thermal conductivities that appear in the diffusion term of the energy equation.

The F terms can be computed in two different ways. For example, we used a centered difference scheme (uniform mesh) for those that are pressure- and temperature-dependent (weakly nonlinear)

$${\left(F_{(p,T)}\right)}_e={\displaystyle \frac{{\left(F_{(p,T)}\right)}_E+{\left(F_{(p,T)}\right)}_P}{2}},$$

(31)

while for the saturation-dependent terms, which are strongly nonlinear, it is done using a first-order upwind scheme. For example, for relative permeability

$$k_{r{\alpha }_e}=\begin{array}{cc}\hfill k_{r{\alpha }_P},& \mathrm{if}{v}_e\ge 0\\ \\ \hfill k_{r{\alpha }_E},& \mathrm{if}{v}_e<0\end{array}.$$

(32)

That is, it is the calculation of the relative permeabilities on the faces of the finite volumes. The capillary pressures on the faces are calculated analogously.

After carrying out all the steps recommended by the Finite Volume Method, we obtained the discretized form for the non-wetting phase (oil) [30]

$$\begin{array}{cc}\hfill & \left(\sum _{d\in {\Omega }_{face}}{\mathbb{T}}_{o_d}+M_{o_P}\right)p_{o_P}^{n+1}-{\left({\mathbb{T}}_o{p}_o\right)}_{Tot}^{n+1}-I_{o_{grav}}\hfill \\ \hfill =& M_{o_P}{p}_{o_P}^n-M_{o_T}\left(T^{n+1}-T^n\right)-\left(S_w^{n+1}-S_w^n\right)M_{o_S}+{\left(q_{s{c}_o}\right)}_P^{n+1}\hfill \end{array}$$

(33)

where ${\Omega }_{face}$ refers to the set of faces of the finite volume, the term ${\left({\mathbb{T}}_o{p}_o\right)}_{Tot}^{n+1}$ is given by

$$\begin{array}{cc}\hfill {\left({\mathbb{T}}_o{p}_o\right)}_{Tot}^{n+1}& ={\left({\mathbb{T}}_{o_e}{p}_{o_E}\right)}^{n+1}+{\left({\mathbb{T}}_{o_w}{p}_{o_W}\right)}^{n+1}+{\left({\mathbb{T}}_{o_s}{p}_{o_S}\right)}^{n+1}\hfill \\ \hfill & +{\left({\mathbb{T}}_{o_n}{p}_{o_N}\right)}^{n+1}+{\left({\mathbb{T}}_{o_b}{p}_{o_B}\right)}^{n+1}+{\left({\mathbb{T}}_{o_a}{p}_{o_A}\right)}^{n+1}\hfill \end{array}$$

(34)

and with respect to gravitational effects,

$$I_{o_{grav}}={\mathbb{T}}_{o_a}^{n+1}{\gamma }_{o_a}^{n+1}\left(Z_P-Z_A\right)+{\mathbb{T}}_{o_b}^{n+1}{\gamma }_{o_b}^{n+1}\left(Z_B-Z_P\right).$$

(35)

The terms coming from the transient term are given by

$$M_{o_P}={\displaystyle \frac{V_b}{\Delta t}}{\left[\left(1-S_w^n\right)\left({\displaystyle \frac{c_{\varphi }{\varphi }^0}{B_o^n}}+{\displaystyle \frac{{\varphi }^{n+1}{c}_n}{B_o^0}}\right)\right]}_P,$$

(36)

$$M_{o_T}=-{\displaystyle \frac{V_b}{\Delta t}}{\left[\left(1-S_w^n\right)\left({\displaystyle \frac{c_{\varphi T}{\varphi }^0}{B_o^n}}+{\displaystyle \frac{{\varphi }^{n+1}{c}_{oT}}{B_o^0}}\right)\right]}_P,$$

(37)

$$M_{o_S}=-{\displaystyle \frac{V_b}{\Delta t}}{\left({\displaystyle \frac{\varphi }{B_o}}\right)}_P^{n+1},$$

(38)

where $V_b=\Delta x_p\Delta y_P\Delta z_P$.

Continuing, we present the discretized equation for the wetting phase (water) [30]

$$\begin{array}{cc}\hfill & {\left(\sum _{d\in {\Omega }_{face}}{\mathbb{T}}_{w_d}{\displaystyle \frac{d{p}_{c_d}}{d{S}_{w_d}}}-M_{w_S}\right)}^{n+1}{S}_{w_P}^{n+1}-{\left(\sum _{\begin{array}{c}d\in {\Omega }_{face}\\ D\in {\Omega }_{vizinho}\end{array}}{S}_{w_D}{\mathbb{T}}_{w_d}{\displaystyle \frac{d{p}_{c_d}}{d{S}_{w_d}}}\right)}^{n+1}=M_{w_S}{S}_w^n\hfill \\ \hfill +& {\left(q_{s{c}_w}\right)}_P^{n+1}-{\left(\sum _{\begin{array}{c}d\in {\Omega }_{face}\end{array}}{\mathbb{T}}_{w_d}\right)}^{n+1}{p}_{o_P}^{n+1}+{\left(\sum _{\begin{array}{c}d\in {\Omega }_{face}\\ D\in {\Omega }_{neighbor}\end{array}}{\mathbb{T}}_{w_d}{p}_{o_D}\right)}^{n+1}-\left(T^{n+1}-T^n\right)M_{w_T}\hfill \\ \hfill -& \left(p_o^{n+1}-p_o^n\right)M_{w_P}+I_{w_{grav}},\hfill \end{array}$$

(39)

where ${\sum }_{\begin{array}{c}d\in {\Omega }_{face}\\ D\in {\Omega }_{vizinho}\end{array}}$ refers to the summation on each face d and, at the same time, on each centroid of the neighboring cell D, for the same spatial direction. Furthermore, the gravity term is given by

$$I_{w_{grav}}={\mathbb{T}}_{w_b}^{n+1}{\gamma }_{w_b}^{n+1}\left(Z_B-Z_P\right)+{\mathbb{T}}_{w_a}^{n+1}{\gamma }_{w_a}^{n+1}\left(Z_A-Z_P\right),$$

(40)

and the other ones are

$$M_{w_P}={\displaystyle \frac{V_b}{\Delta t}}{\left[S_w^n\left({\displaystyle \frac{c_{\varphi }{\varphi }^0}{B_w^n}}+{\displaystyle \frac{{\varphi }^{n+1}{c}_w}{B_w^0}}\right)\right]}_P,$$

(41)

$$M_{w_T}=-{\displaystyle \frac{V_b}{\Delta t}}{\left[S_w^n\left({\displaystyle \frac{c_{\varphi T}{\varphi }^0}{B_w^n}}+{\displaystyle \frac{{\varphi }^{n+1}{c}_{wT}}{B_w^0}}\right)\right]}_P,$$

(42)

$$M_{w_S}={\displaystyle \frac{V_b}{\Delta t}}{\left({\displaystyle \frac{\varphi }{B_w}}\right)}_P^{n+1}.$$

(43)

Finally, we present the discretized form of the energy equation [30]

$$\begin{array}{cc}\hfill & {\left[\left(\rho c_p\right){\displaystyle \frac{V_b}{\Delta t}}+{\Gamma }^{*}\right]}^{n+1}{T}_P^{n+1}-{\Gamma }_e^{n+1}{T}_E^{n+1}-{\Gamma }_w^{n+1}{T}_W^{n+1}-{\Gamma }_s^{n+1}{T}_S^{n+1}-{\Gamma }_n^{n+1}{T}_N^{n+1}\hfill \\ \hfill -& {\Gamma }_b^{n+1}{T}_B^{n+1}-{\Gamma }_a^{n+1}{T}_A^{n+1}=I_{advec}^{n+1}-I_{grav}^{n+1}+I_{source}^{n+1}+{\displaystyle \frac{V_b}{\Delta t}}{\left(\rho c_p\right)}^n{T}_P^n\hfill \end{array}$$

(44)

where

$${\Gamma }^{*}=\sum _{d\in {\Omega }_{face}}{\Gamma }_d,{\Gamma }_d={\displaystyle \frac{A_d{K}_d}{\Delta x_d}},$$

(45)

and the advective, gravitational and source terms are given, respectively, by

$$\begin{array}{cc}\hfill I_{advec}^{n+1}& ={\left[\left(\sum _{d\in {\Omega }_{face}}{\rho }_{s{c}_w}{h}_{w_d}{\mathbb{T}}_{w_d}\right)p_{c_P}\right]}^{n+1}-{\left[\left(\sum _{d\in {\Omega }_{face}}{\rho }_{s{c}_o}{h}_{o_d}{\mathbb{T}}_{o_d}+{\rho }_{s{c}_w}{h}_{w_d}{\mathbb{T}}_{w_d}\right)p_{o_P}\right]}^{n+1}\hfill \\ \hfill & +{\left[\sum _{\begin{array}{c}d\in {\Omega }_{face}\\ D\in {\Omega }_{neighbor}\end{array}}\left({\rho }_{s{c}_o}{h}_{o_d}{\mathbb{T}}_{o_d}+{\rho }_{s{c}_w}{h}_{w_d}{\mathbb{T}}_{w_d}\right)p_{n_D}\right]}^{n+1}-{\left({\rho }_{s{c}_w}{h}_{w_d}{\mathbb{T}}_{w_d}{p}_{c_P}\right)}^{n+1},\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill I_{grav}^{n+1}& ={\left({\gamma }_{n_b}{\mathbb{T}}_{n_b}{\rho }_{s{c}_n}+{\gamma }_{w_b}{\mathbb{T}}_{w_b}{\rho }_{s{c}_w}\right)}^{n+1}\left(Z_B-Z_P\right)\hfill \\ \hfill & +{\left({\gamma }_{n_a}{\mathbb{T}}_{n_a}{\rho }_{s{c}_n}+{\gamma }_{w_a}{\mathbb{T}}_{w_a}{\rho }_{s{c}_a}\right)}^{n+1}\left(Z_A-Z_P\right),\hfill \end{array}$$

(47)

$$I_{source}^{n+1}=Q_H+{\left({\rho }_o{h}_o{Q}_{s{c}_o}\right)}^{n+1}+{\left({\rho }_w{h}_w{Q}_{s{c}_w}\right)}^{n+1}.$$

(48)

It is worth noting that the nonlinear system formed by these three equations is coupled since the oil pressure and water saturation depend on the average temperature value.

3.1. Operator Splitting and Linearization

Since we decided to use linear algebra methods, we must linearize the system of discretized equations. To do so, we opted for the Picard Method [37],

$${\mathbb{T}}^{n+1}\approx {\mathbb{T}}^{n+1,v},$$

(49)

where ${\mathbb{T}}^{n+1,v}$ represents, for example, the transmissibility evaluated at time $n+1$, but at the previous iterative level v. On the other hand, the dependent variables are calculated at $n+1,v+1$, where $v+1$ is the subsequent iterative level. Therefore, we have an external iterative process linked to the Picard Method and an internal one for the numerical resolution of the system of discretized equations.

However, in the present work, we employed an Operator Splitting Method [16,38,39] to enable the unknowns of the system to be determined by solving the equations sequentially. Thus, it was possible to arbitrate the most appropriate algebraic methods for each of them [36].

We established the following resolution order [30]:

• firstly, the subsystem that will provide the pressure of the non-wetting phase is solved (Equation (33));
• next, the saturation of the wetting phase is obtained (Equation (39));
• finally, the average temperature of the reservoir is determined (Equation (44)).

3.2. Resolution of Linear Algebraic Systems

Due to the strategy employed in solving the governing equations, the pressure and average temperature are determined via the Preconditioned Conjugate Gradient Method [4], while the saturation is obtained using the Preconditioned Stabilized Biconjugate Gradient Method [40].

The choice of preconditioned methods was based on the possibility of accelerating the convergence process of iterative methods. To this end, preconditioners constructed from the incomplete LU Factorization process were used [41].

4. Numerical Results

This section aims to present and discuss the results obtained, which are the values calculated for the variables that depend on oil pressure, water saturation, and average reservoir temperature. Before starting the discussion, we first validated the numerical method.

4.1. Numerical Verification

We focused on a non-isothermal one-dimensional flow problem to verify the accuracy of the results for two-phase flow. To this end, we used the analytical solution of a simplified dimensionless problem, the soil remediation process, which consists of injecting liquid water into a porous medium to remove the oil in the soil. The fluids are immiscible, and we did not consider the effect of capillarity on the flow.

The injection of heated water begins, and the effects of capillary pressure are disregarded. The fluids have approximately the same thermal capacity, and we neglected energy transfer by diffusion. Therefore, after simplifications [42]

$${\displaystyle \frac{\partial T}{\partial t}}+{\displaystyle \frac{\partial T}{\partial x}}=0,-\infty <x<\infty ,$$

(50)

subject to a Riemann-type initial condition

$$T(x,0)=\begin{array}{cc}\hfill 1,& x<0\\ \hfill 0,& x>0\end{array}.$$

(51)

Therefore, the solution in terms of temperature is a discontinuity that advances through the porous region. We highlight that viscosity influences the flow pattern since it is a temperature-dependent property [42]

$${\mu }_n\left(T\right)=\left(1-\delta \right)exp\left(T-1\right)+\delta ,$$

(52)

where ${\mu }_n$ is the viscosity of the non-wetting phase, $\delta$ is equal to 1/523, and the viscosity of the wetting phase is assumed to be 1 [42].

To determine the saturation (S), we must consider the establishment of two regions separated by the temperature discontinuity interface, one on the left (higher temperature) and the other on the right (lower temperature). Therefore [42],

$${\displaystyle \frac{\partial S_1}{\partial t}}+{\displaystyle \frac{\partial f_1}{\partial x}}=0,-\infty <x<t,$$

(53)

for the first region, subscript 1, and

$${\displaystyle \frac{\partial S_2}{\partial t}}+{\displaystyle \frac{\partial f_2}{\partial x}}=0,t<x<\infty ,$$

(54)

for the second region, subscript 2, and they must be solved considering the following initial condition

$$S(x,0)=\begin{array}{cc}\hfill 1,& x<0\\ \hfill 0,& x>0\end{array}.$$

(55)

In Equations (53) and (54), $f_1$ and $f_2$ represent the respective flux functions that, for the adopted assumptions [42], are given by

$$f_1={\displaystyle \frac{S^3}{S^3+\delta {\left(1-S\right)}^3}},x<t$$

(56)

and

$$f_2={\displaystyle \frac{S^3}{S^3+\delta exp\left(1\right){\left(1-S\right)}^3}},x>t.$$

(57)

The saturation and temperature profiles can be observed in Figure 3 and Figure 3, respectively, obtained with a grid containing 1,000 cells and a time increment equal to 0.0025. The figures also represent the analytical solutions determined by Moyles [42] for t=0.0625.

For temperature, a comparison between numerical and analytical results shows that the simulator correctly captured the jump that occurs at $x=t$. Nevertheless, due to numerical diffusion, there is indeed a slight harmful noticeable effect. However, spurious oscillations are not present, which, according to Pletcher et al. [43], would be undesirable when solving a problem containing discontinuities in its solution, which must be reproduced numerically as accurately as possible. The temperature values are independent of those of saturation, but they modify them since the viscosity varies with the temperature.

Regarding the saturation profile, it is possible to see that, in general, we observed the same behavior for both the numerical and analytical values. In the first rarefaction, they are superimposed. However, at the end of it, it is noticed that there is a deviation, and the numerical results are below the predicted value, similar to what occurs numerically in Moyles [42]. After the first rarefaction, the saturation value remains constant until the appearance of the first shock. In the case of the analytical solution, we can see the discontinuity, but it is smoothed out due to numerical diffusion in our results. It is worth noting that the appearance of the shock is due to the change in the fluid viscosity, showing the existence of coupling of the equations. After the formation of the first shock, a second rarefaction begins, followed by the appearance of the last shock.

Despite the influence of numerical diffusion, we consider that the numerical results adequately reproduced those predicted by theory.

4.2. Sensitivity analysis

Continuing, we move on to the study of the thermal oil recovery process. We considered the problem of injecting heated water into a reservoir with a slab-type geometry. This procedure aims to sweep the oil that is still in the reservoir towards the production region, besides reducing its viscosity and favoring the flow of fluids in this region. It is a fundamental case for sensitivity analysis. It has significant relevance in the hydrocarbon recovery process, highlighting the contribution of the effects of heating in the recovery of, for example, heavy oil.

This study allows us to learn what happens in a complex process of injecting heated water into a heavy oil reservoir. Therefore, to better understand it, we focus on understanding the influence of the effect of the variation of the water injection temperature ($T_{inj}$), oil viscosity, and the water injection rate (${\partial p/\partial x|}_{x=0}$).

We used the listed aspects as an integral part of the sensitivity analysis process. The standard model for the simulations includes a reservoir containing water and oil, with a three-dimensional geometry in the form of a parallelepiped, whose dimensions are defined by the values of its length, width, and height, respectively $L_x$, $L_y$, and $L_z$ provided in Table 1. In this table, $k_x$, $k_y$, and $k_z$ represent the permeability values in the directions of the x, y, and z axes, respectively, and ${\varphi }_{ini}$ is the initial value of the reservoir porosity.

Table 2, Table 3 and Table 4 also contain several parameters used in the simulations, the default case, besides the properties of the fluids and the rock. Unless explicitly mentioned, these were the values used in all simulations. In Table 2, $\Delta t_i$ is the initial time step, $\Delta t_f$ is the final time step, and $F_{{\Delta }_t}$ is the rate of increase of the time increment ($\Delta t^{n+1}=F_{{\Delta }_t}\Delta t^n$). Furthermore, $p_{ini}$ indicates the initial pressure at the top of the reservoir, $S_w^0$ the initial water saturation, $t_{max}$ the maximum production time, and $T_{r_{ini}}$ the initial temperature of the rock. The properties and parameters used in the simulations were based on those available in the works of the authors de Freitas [41] and dos Santos Heringer et al. [39].

4.2.1. Variation of Injected Fluid Temperature

As mentioned throughout the text, the increase in the temperature of heavy oil causes its viscosity to be reduced and the saturation curve to change. Therefore, we carried out simulations with different values of the water injection temperature to verify the behavior of pressure, saturation, temperature, and viscosity inside the reservoir. Figure 4 contains the pressure values along the x axis, obtained by injecting water at distinct temperatures and for a final production time ($t_{max}$) of 7000 days.

In Figure 4, it is possible to see that the higher the injection temperature, the higher the pressure in the reservoir, and this happens due to the reduction in oil viscosity, allowing an introduction of a greater quantity of water, which causes the reservoir pressure to increase. It is worth noting that in the inlet region on the left, the pressure on the injection face will have very similar values due to the same saturation value, which means the same equivalent relative permeability for each temperature value. At the same time, it is possible to see that there is still a difference between these values, with slightly higher pressures for lower temperatures. It is due to a higher resistance to the entry of water, resulting in a consequent increase in the inlet pressure so that the pressure gradient is maintained.

Regarding temperature, as expected, as water is injected at higher temperatures, the reservoir temperature will increase in all regions, although decreasing as a function of the distance to the injection region, leading to an increase of around 12 K between the minimum and maximum values for regions that are 100 m away from the injection boundary.

Then, from the analysis of the results shown in the figure, we can see that as the temperature of the injection fluid increases, a greater quantity of oil will be recovered, and water will be present in higher amounts in the region close to the injection boundary. We know it is a direct consequence of the decrease in oil viscosity due to the temperature increases because of the transfer of energy from the wetting phase. Ultimately, this is the practical result of thermal recovery via heated fluid injection into the reservoir. It is precisely here that we can attest to the importance of using thermal methods in heavy oil recovery. The reduction in oil viscosity provided by the heating increases sweep efficiency.

4.2.2. Oil Viscosity Variation

From the variation in oil viscosity as a function of the changes in the temperature of the injected fluid (Figure 5), we could confirm the assertions made previously. The heated fluid will cause a decrease in the viscosity of the oil due to the increase in its temperature. Since the injected water will be confined to a region that is less than 20 m from the injection boundary, the decrease in the viscosity of the oil will also be higher in this region. As a consequence, it will increase progressively as the distance increases. The diminished viscosity will increase the mobility of the oil within the reservoir. We assumed that the viscosity of the water remains constant even with the introduction of temperature variation. So, we could disregard its changes concerning those of the oil. The parameter $a_n$ of Equation (4) directly affects the oil viscosity. Consequently, we understood that it was worth studying its influence in determining the pressure, saturation, temperature, and oil viscosity. Unlike the previous case, where viscosity changes were confined in the heated region, we analyzed the effect of the variation in viscosity throughout the entire reservoir. We chose four different values of $a_n$, and the same trends were observed compared to the previous results.

We begin analyzing the pressure curves. It can be seen that the increase in viscosity causes a considerable increase in pressure in the injection region as well as in the rest of the reservoir. This can be explained by the perception that fluid mobility is reduced with the increase in viscosity, requiring an increase in injection pressure to maintain the same pressure gradient. Therefore, the lower the oil viscosity, the lower the flow resistance and the pressure values.

Regarding the temperature results, Figure 5 shows that the higher the viscosity, the lower the temperature of the medium. We can explain this fact due to the reduction of fluid mobility. Consequently, this leads to decreasing velocity and Péclet number. Thus, the effect of heat transfer by advection is reduced in the case of flows of more viscous fluids compared to lower viscosity cases.

In addition, the results show a higher amount of residual oil when the oil is more viscous, which is explained by the higher flow resistance experienced by the more viscous fluid compared to that of lighter oils. Still, regarding the analysis of the variation in saturation, it can be observed, as expected, that for the cases of the flow of more viscous fluids, the advanced front of the injected fluid does not progress as much as in the case of less viscous fluids, advancing more slowly for more viscous oils. On the other hand, for less viscous oils, most distant regions will be reached by the injected fluid in the same period of time.

The results also confirm the effects of the parameter $a_n$ on the oil viscosity. It is clear that the higher its value, the higher the fluid viscosity throughout the entire reservoir. However, due to the heating of the injected fluid near the injection region, they will assume lower values and increase as the fluid is further away from the left border of the reservoir.

It is relevant to highlight that viscosity varies differently as a function of temperature for each specific oil. Therefore, unlike the injection temperature of water, we can not change the way viscosity varies unless the physical properties of the oil are modified by injecting specific products for this purpose.

4.2.3. Injection Rate Variation

The water injection rate is another fundamental parameter for this study since it defines the volume of water injected into the reservoir. The respective curves obtained by changing the injection rate can be found in Figure 6, for the variation of pressure, temperature, saturation, and viscosity.

When considering the behavior of the pressure curves for different values of the pressure gradient imposed on the injection boundary, it is clear that the curves tend to get closer to each other as the injection rate increases. Therefore, it can be seen that due to the values of the reservoir and fluid properties, such as absolute permeability, porosity, and viscosity, there is a limit for the reservoir to produce under the prevailing conditions (initial pressure and temperature, fluid, and rock formation properties), remembering that the pressure is prescribed and kept constant at the production boundary. However, an influence is noted in the region close to the injection boundary due to the increase in the injection rate, which ceases to be relevant as the rate increases.

From Figure 6, regarding the temperature curves, it is not difficult to conclude that the variation in the injection rate of the heated fluid did not cause significant changes in the average temperature values inside the reservoir. Therefore, we should not expect effective changes in the oil viscosity values.

When analyzing the results regarding the saturation curves, the beneficial effect of the initial increase in the injection rate becomes evident, with the consequent increase in the swept efficiency. However, as mentioned, this effect is seen from a gradient corresponding to 15 kPa/m, and the saturation curves become increasingly closer. Therefore, the volume of residual oil in the swept region tends to have the same value for higher injection rates.

For the cases evaluated, based on what has already been explained about the variation in temperature inside the reservoir, we verify that the viscosity curves are practically the same regardless of the injection rate values.

5. Conclusions

In this work, we developed a numerical simulator to solve the problem of non-isothermal three-dimensional two-phase flow in an oil reservoir, and we solved the subsystems for oil pressure, water saturation, and average temperature sequentially. The coefficients of the equations that appear in these subsystems are linearized using the Picard method, and we calculated implicitly pressure and the water saturation. When calculating the average temperature, we did not consider the hypothesis of local thermal equilibrium, and we also used an implicit formulation in time. The effects of tortuosity and hydrodynamic dispersion are taken into account.

The numerical code was validated based on a problem already known in the literature and that has an analytical solution. From its resolution, it was possible to verify the numerical convergence and that the numerical results reproduced those provided by theory. However, we observed a small smoothing due to numerical diffusion caused by the first-order upwind scheme.

Unquestionably, the variation of the injection temperature proved to be the more efficient technique with regard to oil recovery for the cases we tested in the sensitivity analysis, taking into account the value of water saturation at the end of the rarefaction curve. For higher temperatures, there is a reduction in viscosity and the amount of oil found in the region swept by the heated water.

In the other cases, we did not obtain comparable practical gains in swept efficiency despite the reduction in oil viscosity. In the first test case, variation of the injected fluid temperature, the water saturation immediately before the discontinuity, which defines the position of the advancing front, went from approximately 0.55 to 0.65 for the lowest and highest values of water injection temperature, an increase of about 18%.

Overall, the sequential method and the methodologies for discretization, linearization, and solution of algebraic systems proved to be effective when applied to the problem of non-isothermal two-phase flow.