# Iterative Numerical Scheme for Non-Isothermal Two-Phase Flow in Heterogeneous Porous Media
## 4 Spatial Discretization  
The cell-centered finite difference (CCFD) is a type of finite volume method and is locally conservative. In the case of rectangular elements, it was proven early on by [PERSON] (2007) and [PERSON] et al. (2016) to be equivalent to the mixed finite element method. In the CCFD method, the velocity is located on the edges, while the pressure and permeability in the cell centers. Now, we apply the CCFD scheme to the system of Equations (13)-(17), to obtain the full implicit iterative discretization. The discretization form of the pressure equation may be given as  
\[ \begin{split}\mathbf{A}_{a}\left(\mathbf{s}_{w}^{n+1,k},\mathbf{T }^{n+1,k}\right)\mathbf{p}_{w}^{n+1,k+1}+\mathbf{A}_{c}\left(\mathbf{s}_{w}^{ n+1,k},\mathbf{T}^{n+1,k}\right)\mathbf{P}_{d}\left(\mathbf{s}_{w}^{n+1,k+1}, \mathbf{T}^{n+1,k}\right)\\ =\mathbf{A}_{T}\left(\mathbf{T}^{n+1,k}\right)+\mathbf{Q}_{ac}^{ n+1}.\end{split}  \tag{18}\]  
In the above equation, the two matrices, \(\mathbf{A}_{a}\) and \(\mathbf{A}_{c}\), depend on the vectors \(\mathbf{s}_{w}\) and \(\mathbf{T}\) from the previous iteration step. On the other hand, the matrix \(\mathbf{A}_{a}\) depends only on the vector \(\mathbf{T}\) from the previous iteration step. Moreover, the capillary pressure approximation (14) can be discretized as  
\[\mathbf{P}_{d}\left(\mathbf{s}_{w}^{n+1,k+1},\mathbf{T}^{n+1,k}\right)=\mathbf{ P}_{c}\left(\mathbf{s}_{w}^{n+1,k},\mathbf{T}^{n+1,k}\right)+\mathbf{d} \mathbf{P}_{s}\left(\mathbf{s}_{w}^{n+1,k},\mathbf{T}^{n+1,k}\right)\left[ \mathbf{s}_{w}^{n+1,k+1}-\mathbf{s}_{w}^{n}\right], \tag{19}\]  
such that the matrix \(\mathbf{d}\mathbf{P}_{s}\) is diagonally defined, for \(N_{c}\) total number of cells, by  
\[\mathbf{d}\mathbf{P}_{s}\left(\mathbf{s}_{w}^{n+1,k},\mathbf{T}^{n+1,k}\right) =\text{diag}\left(\mathbf{d}\mathbf{P}_{s}\left(\mathbf{s}_{w,h}^{n+1,k}, \mathbf{T}^{n+1,k}_{w,h}\right)\right),\quad h=1,2,...,N_{c}. \tag{20}\]  
It is clear that the derivative of \(p_{c}\) is a function of \(p_{c}\) itself, as the saturation at each spatial point varies with time even though it is discontinuously distributed in space. The spatial discretization of the saturation equation to be coupled with the pressure equation given as  
\[\mathbf{M}\frac{\mathbf{s}_{w}^{n+1,k+1}-\mathbf{s}_{w}^{n}}{\Delta t^{n}}+ \mathbf{A}_{w}\left(\mathbf{s}_{w}^{n+1,k},\mathbf{T}^{n+1,k}\right)\mathbf{P }_{w}^{n+1,k+1}=\mathbf{Q}_{w}^{n+1}, \tag{21}\]  
such that \(\mathbf{M}\) is a diagonal matrix appearing from the porosity coefficient. The coupled pressure equation can be obtained from substituting (20) and (21) into (19), i.e.,  
\[\mathbf{A}_{t}\left(\mathbf{s}_{w}^{n+1,k},\mathbf{T}^{n+1,k}\right)\mathbf{ P}_{w}^{n+1,k+1}=\mathbf{Q}_{t}\left(\mathbf{s}_{w}^{n+1,k},\mathbf{T}^{n+1,k} \right), \tag{22}\]  
where  
\[\mathbf{A}_{t}=\mathbf{A}_{a}-\Delta t^{n}\mathbf{A}_{c}\mathbf{d}\mathbf{P}_ {s}\mathbf{M}^{-1}\mathbf{A}_{w},\]  
and  
\[\mathbf{Q}_{t}=\mathbf{A}_{T}+\mathbf{Q}_{ac}^{n+1}-\mathbf{A}_{c}\left[ \mathbf{p}_{c}+\mathbf{d}\mathbf{P}_{s}\left(\mathbf{s}_{w}^{n}-\mathbf{s}_{w} ^{n+1,k}\right)\right]-\Delta t^{n}\mathbf{A}_{c}\mathbf{d}\mathbf{P}_{s} \mathbf{M}^{-1}\mathbf{Q}_{w}^{n+1}.\]  
The upwind scheme is used to treat the convection term of Equation (16) to obtain  
\[\mathbf{M}\frac{\mathbf{s}_{w}^{n+1,k+1}-\mathbf{s}_{w}^{n}}{\Delta t^{n}}+ \mathbf{A}_{s}\left(\mathbf{s}_{w}^{n+1,k},\mathbf{P}_{w}^{n+1,k+1}\right) \mathbf{f}_{w}\left(\mathbf{s}_{w}^{n+1,k}\right)=\mathbf{Q}_{s}^{n+1}. \tag{23}\]Then, update \(\mathbf{f}_{w}^{n+1,k+1}\) i.e., \(\mathbf{u}_{w}^{n+1,k+1}\) to be used in the energy Equation (17),  
\[