## 3 Numerical Approach and Benchmarking of the Code  
The dimensionless governing equations are written in the scalar form as:  
\[\omega_{1}=-\Big{[}\frac{\partial^{2}\Psi_{1}}{\partial\kappa^{2}}+\frac{\partial ^{2}\Psi_{1}}{\partial\gamma^{2}}+\frac{\partial^{2}\Psi_{1}}{\partial\alpha^{2 }}\Big{]},\omega_{2}=-\Big{[}\frac{\partial^{2}\Psi_{2}}{\partial\kappa^{2}}+ \frac{\partial^{2}\Psi_{2}}{\partial\gamma^{2}}+\frac{\partial^{2}\Psi_{2}}{ \partial\alpha^{2}}\Big{]}\text{ and }\omega_{3}=-\Big{[}\frac{\partial^{2}\Psi_{3}}{\partial\kappa^{2}}+ \frac{\partial^{2}\Psi_{3}}{\partial\gamma^{2}}+\frac{\partial^{2}\Psi_{3}}{ \partial\alpha^{2}}\Big{]} \tag{12}\]  
\[\frac{\partial\omega_{1}}{\partial\alpha}+\frac{\partial}{\partial \kappa}\Big{[}\mathrm{U}_{1}\omega_{1}-\Gamma_{\omega}\frac{\partial\omega_{1 }}{\partial\kappa}\Big{]}+\frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{2} \omega_{1}-\Gamma_{\omega}\frac{\partial\omega_{1}}{\partial\kappa}\Big{]}+ \frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{3}-\Gamma_{\omega}\frac{ \partial\omega_{1}}{\partial\kappa}\Big{]}=\omega_{1}\frac{\partial\omega_{1 }}{\partial\kappa}+\omega_{2}\frac{\partial\Psi_{1}}{\partial\gamma}+\omega_{ 3}\frac{\partial\Omega_{1}}{\partial\alpha} \tag{13}\] \[-\mathrm{Pr}\cdot\mathrm{K}\Big{(}\frac{\rho_{1}}{\rho_{\rm{d}}} \Big{)}\Big{(}\frac{\partial\Psi_{1}}{\partial\kappa^{2}}+\frac{\partial\Psi_ {1}}{\partial\gamma^{2}}+\frac{\partial\Psi_{1}}{\partial\alpha^{2}}\Big{)}+ \mathrm{R}\mathrm{Pr}\Big{(}\frac{\rho_{1}\rho_{\rm{d}}}{\left(\rho\right) \Psi_{1}}\Big{)}\Big{(}\frac{\partial\tau}{\partial\kappa}-\mathrm{N}\frac{ \partial\Gamma}{\delta x}\Big{)}+\mathrm{Pr}\,\mathrm{Ha}^{2}\big{(}\frac{\rho _{1}}{\rho_{\rm{d}}}\Big{)}\Big{(}\frac{\sigma_{1}}{\sigma_{1}}\Big{)}\Big{(} \frac{\partial\upsilon_{1}}{\partial\kappa}\Big{)}\] (14) \[-\mathrm{Pr}\cdot\mathrm{K}\Big{(}\frac{\rho_{1}}{\rho_{\rm{d}}} \Big{)}\Big{(}\frac{\partial\Psi_{1}}{\partial\kappa^{2}}+\frac{\partial\Psi_ {1}}{\partial\gamma^{2}}+\frac{\partial\Psi_{1}}{\partial\alpha^{2}}\Big{)}+ \mathrm{Pr}\,\mathrm{Ha}^{2}\Big{(}\frac{\rho_{1}}{\rho_{\rm{d}}}\Big{)}\Big{(} \frac{\partial\tau}{\partial\kappa}\Big{)}\] \[\frac{\partial\omega_{2}}{\partial\tau}+\frac{\partial}{\partial \kappa}\Big{[}\mathrm{U}_{1}\omega_{3}-\Gamma_{\omega}\frac{\partial\omega_{3 }}{\partial\kappa}\Big{]}+\frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{2} \omega_{1}-\Gamma_{\omega}\frac{\partial\omega_{2}}{\partial\alpha}\Big{]}= \omega_{1}\frac{\partial\upsilon_{1}}{\partial\kappa}+\omega_{2}\frac{\partial \Psi_{1}}{\partial\gamma}+\omega_{3}\frac{\partial\Omega_{1}}{\partial\alpha}\] (15) \[-\mathrm{Pr}\cdot\mathrm{K}\Big{(}\frac{\rho_{1}}{\rho_{\rm{d}}} \Big{)}\Big{(}\frac{\partial\Psi_{1}}{\partial\kappa^{2}}+\frac{\partial\Psi_ {1}}{\partial\gamma}+\frac{\partial\Psi_{1}}{\partial\alpha^{2}}\Big{)}+ \mathrm{R}\mathrm{Pr}\Big{(}\frac{\rho\Psi_{1}}{\left(\rho\right)\Psi_{1}} \Big{)}\Big{(}\frac{\rho_{1}}{\rho_{\rm{d}}}\Big{)}\Big{(}\frac{\pi\tau}{ \alpha}+\mathrm{N}\frac{\partial\chi}{\Delta\kappa}\Big{)}+\mathrm{Pr}\,\mathrm{ Ha}^{2}\big{(}\frac{\rho}{\rho_{\rm{d}}}\Big{)}\Big{(}\frac{\partial \Psi_{1}}{\alpha\kappa}\Big{)}\] (16) \[\frac{\partial\mathrm{C}}{\partial\tau}+\frac{\partial}{\partial \chi}\Big{[}\mathrm{U}_{1}-\Gamma_{\mathrm{C}}\frac{\partial\mathrm{C}}{ \partial\chi}\Big{]}+\frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{2}- \Gamma_{\mathrm{C}}\frac{\partial\mathrm{C}}{\partial\gamma}\Big{]}+\frac{ \partial}{\partial\alpha}\Big{[}\mathrm{U}_{3}-\Gamma_{\mathrm{C}}\frac{ \partial\mathrm{C}}{\partial\chi}\Big{]}=0\] (17) \[\frac{\partial\mathrm{h}_{1}}{\partial\tau}+\frac{\partial}{ \partial\chi}\Big{[}\mathrm{U}_{1}-\Gamma_{\mathrm{h}}\frac{\partial\mathrm{h}_ {1}}{\partial\kappa}\Big{]}+\frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{2} -\Gamma_{\mathrm{h}}\frac{\partial\mathrm{h}_{1}}{\partial\gamma}\Big{]}+\frac{ \partial}{\partial\chi}\Big{[}\mathrm{U}_{3}-\Gamma_{\mathrm{h}}\frac{ \partial\mathrm{h}_{1}}{\partial\chi}\Big{]}=\mathrm{Pr}\cdot\mathrm{K}\Big{(} \frac{\rho_{\rm{f}}}{\rho_{\rm{nf}}}\Big{)}\big{(}\omega_{1}-\mathrm{h}_{1}\big{)}\] (18) \[\frac{\partial\mathrm{h}_{2}}{\partial\tau}+\frac{\partial}{\partial \chi}\Big{[}\mathrm{U}_{1}-\Gamma_{\mathrm{h}}\frac{\partial\mathrm{h}_{2}}{ \partial\chi}\Big{]}+\frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{2}- \Gamma_{\mathrm{h}}\frac{\partial\mathrm{h}_{2}}{\partial\gamma}\Big{]}+\frac{ \partial}{\partial\alpha}\Big{[}\mathrm{U}_{3}-\Gamma_{\mathrm{h}}\frac{ \partial\mathrm{h}_{2}}{\partial\chi}\Big{]}=\mathrm{Pr}\cdot\mathrm{K}\Big{(} \frac{\rho_{\rm{f}}}{\rho_{\rm{nf}}}\Big{)}\big{(}\omega_{2}-\mathrm{h}_{2} \big{)} \tag{19}\]  
One notes that the above equations can be written in the following general form:  
\[\frac{\partial\Phi}{\partial\tau}+\frac{\partial}{\partial\chi}\Big{[}\mathrm{U}_ {1}\Phi-\Gamma_{\Phi}\frac{\partial\Phi}{\partial\chi}\Big{]}+\frac{\partial}{ \partial\gamma}\Big{[}\mathrm{U}_{2}\Phi-\Gamma_{\Phi}\frac{\partial\Phi}{ \partial\gamma}\Big{]}+\frac{\partial}{\partial\chi}\Big{[}\mathrm{U}_{3}\Phi- \Gamma_{\Phi}\frac{\partial\Phi}{\partial\chi}\Big{]}=\mathrm{S}_{\Phi} \tag{20}\]  
With:  
\[\frac{\partial\omega_{1}}{\partial\alpha}+\frac{\partial}{\partial \kappa}\Big{[}\mathrm{U}_{1}\omega_{1}-\Gamma_{\omega}\frac{\partial\omega_{1 }}{\partial\kappa}\Big{]}+\frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{2} \omega_{1}-\Gamma_{\omega}\frac{\partial\omega_{1}}{\partial\kappa}\Big{]}+ \frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{3}-\Gamma_{\omega}\frac{ \partial\omega_{1}}{\partial\kappa}\Big{]}=\omega_{1}\frac{\partial\omega_{1 }}{\partial\kappa}+\omega_{2}\frac{\partial\Psi_{1}}{\partial\gamma}+\omega_{ 3}\frac{\partial\Omega_{1}}{\partial\alpha} \tag{13}\]0  
\(\Gamma_{\Phi}\text{: dimensionless coefficient}\)  
\(\mathrm{S}_{\Phi}\text{: source term}\)  
Using the following expressions:  
\[\frac{\partial\omega_{1}}{\partial\alpha}+\frac{\partial}{\partial \kappa}\Big{[}\mathrm{U}_{1}\omega_{1}-\Gamma_{\omega}\frac{\partial\omega_{1 }}{\partial\kappa}\Big{]}+\frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{2} \omega_{1}-\Gamma_{\omega}\frac{\partial\omega_{1}}{\partial\kappa}\Big{]}+ \frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{3}-\Gamma_{\omega}\frac{ \partial\omega_{1}}{\partial\kappa}\Big{]}=\omega_{1}\frac{\partial\omega_{1 }}{\partial\kappa}+\omega_{2}\frac{\partial\Psi_{1}}{\partial\gamma}+\omega_{ 3}\frac{\partial\Omega_{1}}{\partial\alpha} \tag{13}\]1  
The general equation can be written as:  
\[\frac{\partial\omega_{1}}{\partial\alpha}+\frac{\partial}{\partial \kappa}\Big{[}\mathrm{U}_{1}\omega_{1}-\Gamma_{\omega}\frac{\partial\omega_{1 }}{\partial\kappa}\Big{]}+\frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{2} \omega_{1}-\Gamma_{\omega}\frac{\partial\omega_{1}}{\partial\kappa}\Big{]}+ \frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{3}-\Gamma_{\omega}\frac{ \partial\omega_{1}}{\partial\kappa}\Big{]}=\omega_{1}\frac{\partial\omega_{1 }}{\partial\kappa}+\omega_{2}\frac{\partial\Psi_{1}}{\partial\gamma}+\omega_{ 3}\frac{\partial\Omega_{1}}{\partial\alpha} \tag{13}\]2  
The integration of the previous equation on a control volume leads to (Patankar [47]):  
\[\frac{\partial\omega_{1}}{\partial\alpha}+\frac{\partial}{\partial \kappa}\Big{[}\mathrm{U}_{1}\omega_{1}-\Gamma_{\omega}\frac{\partial\omega_{1 }}{\partial\kappa}\Big{]}+\frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{2} \omega_{1}-\Gamma_{\omega}\frac{\partial\omega_{1}}{\partial\kappa}\Big{]}+ \frac{\partial}{\partial\gamma}\Big{[}\mathrm{U}_{3}-\Gamma_{\omega}\frac{ \partial\omega_{1}}{\partial\kappa}\Big{]}=\omega_{1}\frac{\partial\omega_{1 }}{\partial\kappa}+\omega_{2}\frac{\partial\Psi_{1}}{\partial\gamma}+\omega_{ 3}\frac{\partial\Omega_{1}}{\partial\alpha} \tag{13}\]3  
\(e\), \(w\), \(n\), \(s\), \(b\) and \(f\) represent the faces of control volume centered in \(P\).