# Iterative Numerical Scheme for Non-Isothermal Two-Phase Flow in Heterogeneous Porous Media
## 2 Modeling and Mathematical Formulation  
Consider the problem of two-phase flow with convective heat transfer in heterogeneous porous media. The flow is assumed to be incompressible and immiscible without phase transition. The variation in density is caused by differences in temperature and is related to the buoyancy effect which is described by the Boussinesq approximation. The basis of this approximation is that there are flows in which the temperature varies slightly, in these, the density varies slightly, and this buoyancy contributes the motion. Therefore, the variation in density is neglected anywhere except in the buoyancy term. Let \(\rho_{\alpha,0}\) be a reference density of the phase \(\alpha\), the temperature is \(T_{r}\). For small temperature difference we can write  
\[\rho_{\alpha}=\rho_{\alpha,0}(1-\beta_{\alpha}(T-T_{r})),\quad\beta_{\alpha}=- \frac{1}{\rho_{\alpha,0}}\frac{\partial\rho_{\alpha}}{\partial T},\quad\alpha =w,n,\]  
where \(\rho_{\alpha,0}\) is a reference density of the phase \(\alpha\). \(\beta_{\alpha}\) is the coefficient of volume expansion. \(\beta_{\alpha}\) is in the range \(10^{-3}\)-\(10^{-4}\). For a temperature variation of moderate amount, we have  
\[\frac{\delta\rho_{\alpha}}{\rho_{\alpha,0}}=\frac{|\rho_{\alpha}-\rho_{\alpha,0}|}{\rho_{\alpha}}=\beta_{\alpha}|T-T_{r}|<<1.\]  
Thus, the variation in the density and therefore in the specific volume is negligible. The velocity equation of the phase \(\alpha\) is obtained by [PERSON]'s law,  
\[\mathbf{u}_{\alpha}=-\frac{k_{ra}\mathbf{K}}{\mu_{\alpha}}\
abla(p_{\alpha}+g \rho_{\alpha,0}(1-\beta_{\alpha}(T-T_{r}))z),\qquad\alpha=w,n, \tag{1}\]  
where \(k_{ra}\) is the relative permeability of the phase \(\alpha\); \(\mathbf{K}\) is the permeability tensor; \(\mu_{\alpha}\) is the \(\alpha\)-fluid viscosity; \(p_{\alpha}\) is the phase-\(\alpha\) pressure; \(g\) is the gravitational acceleration; \(z\) is the depth; \(q_{\alpha}\) is the external mass flow rate; \(w\) refers to the wetting phase; and \(n\) refers to the non-wetting phase. The temperature difference may cause a convective flow, which is represented by the Boussinesq approximation in the body force term of the previous Equation (1). For measures of central tendency (mean, mode, median), pore-throat sizes (diameters) are generally greater than 2 mm in conventional reservoir rocks, they range from about 2 to 0.03 mm in tight-gas sandstones, and range from 0.1 to 0.005 mm in shales [19]. On the other hand, the ratio between the inertial forces and the viscous forces driving the flow is computed by the Reynolds number, which is used as a criterion to distinguish between the laminar flow, the turbulent flow, and the transition zone. For porous media, the Reynolds number is defined as \(Re=\rho u_{i}d_{i}/\mu\), where \(u_{i}\) is the injection velocity and \(d_{i}\) is the injection width of the inlet. According to [PERSON] [20], [PERSON]'s law, which supposes a laminar flow, is valid for Reynolds number less than 1, but the upper limit can be extended up to 10. This may be because the Reynolds number is calculated based on the discharge (inlet) flow rate, however, the velocity reduces significantly in the porous media as it depends on the permeability. For the current problem, the Reynolds number is about 0.6-1.8 for the injected water depending on the discharge inlet diameter 0.01-0.03 m. Therefore, based on [PERSON]'s classification, [PERSON]'s law is valid for our study. Moreover, according to [PERSON]'s classification, the nonlinear effects are taken into consideration for 10 < Re < 100 [20]. The flow in the current problem is linear, therefore, no need to consider the Brinkman and [PERSON] effects. The saturation equation which represents the mass conservation law can be written as  
\[\phi\frac{\partial s_{\alpha}}{\partial t}+\
abla\cdot\mathbf{u}_{\alpha}=q_{ \alpha},\quad\sum_{\alpha}s_{\alpha}=1,\quad\alpha=w,n, \tag{2}\]  
where \(\phi\) [-] is the porosity, and \(s_{\alpha}\) is the saturation of the phase \(\alpha\).