# Iterative Numerical Scheme for Non-Isothermal Two-Phase Flow in Heterogeneous Porous Media
## 2 Modeling and Mathematical Formulation
If we sum the two equations of saturation, (2), namely, \(\alpha=w\) and \(\alpha=n\), whilst implementing the constraint \(s_{w}+s_{n}=1\), we reach the following pressure equation  
\[\
abla\cdot\mathbf{u}_{t}=-\
abla\cdot\lambda_{t}\mathbf{K}\
abla \Big{(}p_{w}+g\rho_{w,0}\left(1-\beta_{w}(T-T_{r})\right)z\Big{)}-\
abla \cdot\lambda_{n}\mathbf{K}\] \[\
abla\Big{(}p_{c}+g\Delta\rho_{0}\left(1-\beta_{n}(T-T_{r}) \right)z\Big{)}=q_{t}, \tag{3}\]  
where \(\mathbf{u}_{t}=\mathbf{u}_{w}+\mathbf{u}_{n}\) is the total velocity; \(p_{c}=p_{n}-p_{w}\) is the capillary pressure; \(\Delta\rho_{0}=\rho_{n,0}-\rho_{w,0}\) is the density difference; \(\lambda_{\alpha}=k_{\alpha}/\mu_{a}\) is the phase mobility; \(\lambda_{t}=\lambda_{w}+\lambda_{n}\) is the total mobility; and \(q_{t}=q_{w}+q_{n}\) is the total source mass transfer. Equation (3) with the following saturation equation are solved simultaneously to compute the pressure of water phase:  
\[\phi\frac{\partial s_{w}}{\partial t}=-\
abla\cdot\lambda_{w}\mathbf{K}\
abla (p_{w}+g\rho_{w,0}(1-\beta_{w}(T-T_{r}))\
abla z)+q_{w}. \tag{4}\]  
Therefore, we can update the saturation using the following equation:  
\[\phi\frac{\partial s_{w}}{\partial t}-q_{w}=-\
abla\cdot(f_{w}\mathbf{u}_{a}), \tag{5}\]  
where \(f_{w}=\lambda_{w}/\lambda_{t}\) is the flow fraction, \(\mathbf{u}_{w}=f_{w}\mathbf{u}_{a}\) and \(\mathbf{u}_{a}=-\lambda_{t}\mathbf{K}\
abla p_{w}+g\rho_{w,0}(1-\beta_{w}(T-T _{r}))\
abla z\). The relative permeability of the two phases is defined as  
\[k_{rw}=k_{rw}^{0}\xi^{2},\quad k_{rn}=k_{rn}^{0}\left(1-S\right)^{2},\]  
where \(S=\frac{s_{w}-s_{wr}}{1-s_{wr}-s_{wr}}\) is the normalized wetting phase saturation, such that, \(s_{wr}\) is the residual non-wetting saturation, \(s_{wr}\) is the irreducible saturation, \(0\leq S\leq 1\). In addition, \(k_{rw}^{0}=k_{rw}\left(S=1\right),k_{rw}^{0}=k_{rw}\left(S=1\right)\). Assume that the two fluid phases and the solid phase are of local thermal equilibrium of low flow velocities and relatively small grain sizes (cf. [PERSON] [21]):  
\[T_{matrix}=T_{fluids}=T.\]  
This yields the conservation equation of energy for a two-phase system, consisting of a solid phase plus fluid phases [16; 22]:  
\[\frac{\delta\rho_{\alpha}}{\rho_{\alpha,0}}=\frac{|\rho_{\alpha}-\rho_{\alpha,0}|}{\rho_{\alpha}}=\beta_{\alpha}|T-T_{r}|<<1.\]0 \[\frac{\delta\rho_{\alpha}}{\rho_{\alpha,0}}=\frac{|\rho_{\alpha}-\rho_{\alpha,0}|}{\rho_{\alpha}}=\beta_{\alpha}|T-T_{r}|<<1.\]1  
where \(\rho_{s}\) [kg/m\({}^{3}\)] is density of the solid phase; \(c_{p,s}\) [kg/m\({}^{3}\)] is the heat capacity of the solid phase; \(c_{p,\alpha}\) [J/Kg\(\cdot\)K] is the heat capacity of the fluid phase \(\alpha\); \(T\) [K] is the temperature; \(Q_{T}\) [m\({}^{3}\)/s] is the heat source term; \(\tau\) is the tortuosity of the water phase; and \(h_{s}\) [J/K\(\cdot\)m\(\cdot\)s] is the thermal conductivity of the solid phase. The thermal conductivity of the phase \(\alpha\) is expressed as  
\[\frac{\delta\rho_{\alpha}}{\rho_{\alpha,0}}=\frac{|\rho_{\alpha}-\rho_{\alpha,0}|}{\rho_{\alpha}}=\beta_{\alpha}|T-T_{r}|<<1.\]2  
where \(h_{\alpha,d}\) is the deterministic thermal conductivity of the phase \(\alpha\), while the mechanical thermal-dispersion tensor, \(h^{\text{disp}}\), may be written as  
\[\frac{\delta\rho_{\alpha}}{\rho_{\alpha,0}}=\frac{|\rho_{\alpha}-\rho_{\alpha,0}|}{\rho_{\alpha}}=\beta_{\alpha}|T-T_{r}|<<1.\]3  
Therefore,  
\[\frac{\delta\rho_{\alpha}}{\rho_{\alpha,0}}=\frac{|\rho_{\alpha}-\rho_{\alpha,0}|}{\rho_{\alpha}}=\beta_{\alpha}|T-T_{r}|<<1.\]4  
where \(d_{l,w}\) is the longitudinal dispersion coefficient and \(d_{t,w}\) is the transverse dispersion coefficient.