# Iterative Numerical Scheme for Non-Isothermal Two-Phase Flow in Heterogeneous Porous Media
## 3 Iterative Method  
The governing Equations (3)-(6), are strongly coupled and highly nonlinear. In order to treat such a type of complexity, we may use an iterative method [24]. The time interval \([0,T]\) is discretized into a number of \(N_{T}\) time steps, \(\Delta t^{n}=t^{n+1}-t^{n}\), \(n+1\) refers to the current time step, and \(n\) refers to the previous time step. Apply iterations of \(N_{I}\) iterations for each time step, such that \(k+1\) refers to the current iteration step and \(k\) refers to the previous iteration step. The time discretization of time derivatives is achieved by using the backward Euler. The discretized pressure equation takes the form  
\[-\
abla\cdot\lambda_{t}\left(s_{w}^{n+1,k}\right)\mathbf{K} \
abla(p_{w}^{n+1,k+1}+g\rho_{w,0}(1-\beta_{w}(T^{n+1,k}-T_{r}))\
abla z)- \
abla\cdot\] \[\lambda_{n}\left(s_{w}^{n+1,k}\right)\mathbf{K}\
abla\left(p_{c} ^{n+1,k+1}+g\Delta\rho_{0}(1-\beta_{w}(T^{n+1,k}-T_{r}))\
abla z\right)=q_{t}. \tag{13}\]  
The capillary pressure \(p_{c}^{n+1,k+1}\) is linearized as  
\[p_{c}\left(s_{w}^{n+1,k+1},T^{n+1,k}\right)\cong p_{c}\left(s_{w}^{n+1,k},T^{n +1,k}\right)+\frac{\partial p_{c}}{\partial s_{w}}|_{s_{w}^{n+1,k},T^{n+1,k}} \left(s_{w}^{n+1,k+1}-s_{w}^{n}\right). \tag{14}\]  
Using the saturation equation to find  
\[s_{w}^{n+1,k+1}-s_{w}^{n}=-\frac{\Delta t^{n}}{\phi}\left[\
abla\cdot\lambda_{ w}^{n+1,k}\mathbf{K}\
abla(p_{w}^{n+1,k+1}+g\rho_{w,0}(1-\beta_{w}(T^{n+1,k}-T_{r})) \
abla z)\right]. \tag{15}\]  
The above three equations are solved to compute the wetting phase pressure. Then, saturation is updated using the explicit equation  
\[\phi\frac{s_{w}^{n+1,k+1}-s_{w}^{n}}{\Delta t^{n}}+\alpha_{r}\
abla\cdot\left( f_{w}^{n+1,k}\mathbf{u}_{a}^{n+1,k+1}\right)=q_{w}^{n+1}, \tag{16}\]  
where  
\[\mathbf{u}_{a}^{n+1,k+1}=-\lambda_{t}^{n+1,k}\mathbf{K}\
abla p_{w}^{n+1,k+1 }+g\rho_{w,0}(1-\beta_{w}(T^{n+1,k}-T_{r}))\
abla z.\]  
In the above equation, the advection term is multiplied with the relaxation factor, \(\alpha_{r}\). In order to ensure the stability of the proposed scheme, we employ the stability Courant-Friedrichs-Lewy (CFL) condition (CFL \(<\) 1) to estimate a relaxation factor. The convection terms in the discretized saturation and temperature equations are multiplied by the following relaxation factor:  
\[\alpha_{r}=\min(\frac{c_{1}}{\max(\text{CFL})},c_{2}),\quad 0\leq c_{1},c_{2} \leq 1,\quad\alpha_{r}\text{CFL}<1.\]  
The energy equation is computed implicitly by the following scheme:  
\[\frac{1}{\Delta t^{n}}\left[\left((1-\phi)\rho_{s}c_{p,s}+\phi \sum_{\alpha}s_{\alpha}^{n+1,k+1}\rho_{\alpha}(T^{n+1,k})c_{p,\alpha}\right) T^{n+1,k+1}-\right.\] \[\left.\left((1-\phi)\rho_{s}c_{p,s}+\phi\sum_{\alpha}s_{\alpha}^ {n}\rho_{\alpha}(T^{n})c_{p,\alpha}\right)T^{n}\right]+\alpha_{r}\sum_{\alpha} \rho_{\alpha}\mathbf{u}_{a}^{n+1,k+1}\cdot\
abla T^{n+1,k+1}=\
abla\cdot \tag{17}\] \[\left((1-\phi)h_{s}+\phi\sum_{\alpha}\tau s_{\alpha}^{n+1,k+1}h_{ \alpha}^{n+1,k+1}\right)\
abla T^{n+1,k+1}+Q_{T}^{n+1},\quad\alpha=w,n,\]where  
\[h_{\alpha}^{n+1,k+1}=(\phi s_{w}^{n+1,k+1}\tau h_{\alpha,d}+d_{t,w}|\mathbf{u}_{w }|^{n+1,k+1})\mathbf{I}+(d_{l,w}-d_{t,w})\frac{\mathbf{u}_{w}^{n+1,k+1}\mathbf{u }_{w}^{n+1,k+1}}{|\mathbf{u}_{w}|^{n+1,k+1}}.\]  
The advection term in the energy Equation (17) is also multiplied by the relaxation factor, \(\alpha_{r}\). Finally, the variables, \(\lambda_{w}\),\(\lambda_{n}\),\(\lambda_{t}\), \(f_{w}\), and \(\rho_{\alpha}\), are updated. This procedure is repeated until approaching convergence.