## 3 Numerical Solution
#### Numerical Algorithm  
To solve Equations (6)-(8) subject to the boundary conditions given in Equation (10) and the initial conditions given in Equation (9), we employ finite difference numerical techniques. In particular, the computational algorithms are based on the semi-implicit finite difference methods. Details of our numerical algorithms are illustrated, say, in [21; 22]. Taking \(N\) as the current time level, \(N+1\) as the subsequent time level, and \(N+\xi\) with \(0\leq\xi\leq 1\) as an intermediate time level, the semi-implicit scheme for the \(u\)-velocity component reads,  
\[\begin{split}\frac{u^{(N+1)}-u^{(N)}}{\Delta t}+\left[u\frac{ \partial u}{\partial x}+v\frac{\partial u}{\partial y}\right]^{(N)}& =\frac{1}{\text{Re}}\left[\eta(\dot{\gamma},\Theta)\right]^{(N)} \left[\frac{\partial^{2}u}{\partial y^{2}}\right]^{(N+\xi)}\\ &+\frac{1}{\text{Re}}\left[\frac{\partial}{\partial y}\eta(\dot{ \gamma},\Theta)\,\frac{\partial u}{\partial y}\right]^{(N)}+\frac{\text{Gr}}{ \text{Re}^{2}}\,\Theta^{(N)}.\end{split} \tag{11}\]  
The equation for \(u^{(N+1)}\) then becomes:  
\[-r_{1}\,u_{i,j-1}^{(N+1)}+\left(1+2r_{1}\right)u_{i,j}^{(N+1)}-r_{1}\,u_{i,j+ 1}^{(N+1)}=\text{explicit terms}, \tag{12}\]where  
\[r_{1}=\frac{\xi}{\text{Re}}\left[\eta(\gamma,\Theta)\right]_{i,j}^{(N)}\frac{ \Delta t}{\Delta y^{2}}.\]  
The solution procedure for \(u^{(N+1)}\) reduces to efficient linear algebraic processes of inversion of nonsingular, diagonally dominant, and tri-diagonal matrices. The semi-implicit scheme for the temperature equation is similarly approached; unmixed second partial derivatives of the temperature are treated implicitly,  
\[\begin{split}\frac{\Theta^{(N+1)}-\Theta^{(N)}}{\Delta t}+\left[u \frac{\partial\Theta}{\partial x}+v\frac{\partial\Theta}{\partial y}\right]^{( N)}&=\frac{1}{\text{Re}\,\text{Pr}}\left[\frac{\partial^{2}\Theta}{ \partial x^{2}}+\frac{\partial^{2}\Theta}{\partial y^{2}}\right]^{(N+\xi)}\\ &+\frac{1}{\text{Re}\,\text{Pr}}\left[\eta(\gamma,\Theta)\left( \frac{\partial u}{\partial y}\right)^{2}\right]^{(N)}.\end{split} \tag{13}\]  
The algorithm (13) allows for the decoupling of variables, which represents a significant advantage over a fully implicit scheme. The implicit terms are,  
\[1-\frac{1}{\text{Re}\,\text{Pr}}\,\xi\,\Delta t\left(\frac{\partial^{2}}{ \partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}\right)\Theta^{(N+1)}.\]  
The operator on \(\Theta^{(N+1)}\) factorizes with an associated error term of order,  
\[\text{factorization error}=O\Bigg{(}\left(\frac{\xi\,\Delta t}{\text{Re}\, \text{Pr}}\right)^{2}\Bigg{)}. \tag{14}\]  
Small \(\Delta t\) and relatively large \(\text{Re}\,\text{Pr}\) lead to very small factorization errors. The semi-implicit algorithm for the temperature equation, therefore, reduces to,  
\[\left(1-\frac{1}{\text{Re}\,\text{Pr}}\,\xi\,\Delta t\,\frac{\partial^{2}}{ \partial x^{2}}\right)\left(1-\frac{1}{\text{Re}\,\text{Pr}}\,\xi\,\Delta t \,\frac{\partial^{2}}{\partial y^{2}}\right)\Theta^{(N+1)}=\text{explicit terms}. \tag{15}\]  
The solution procedure again reduces to the efficient linear algebraic processes of inversion of nonsingular, diagonally dominant, and tri-diagonal matrices.