## 3 Numerical Model  
### Thermal Lattice Boltzmann Equations  
The TLBM consists of simulating the statistical behavior of a set of particles on a lattice with finite velocities. It stems from the discrete Boltzmann equation and allows providing macroscopic fluid properties, such as density, velocity, pressure, _etc._, through weighted averages, or moments, of the particle distribution for all discrete lattice velocities. The SRT-LBM (also called the Bhatnagar-Gross-Krook (BGK) model) for incompressible thermal flows builds on two distribution functions (DFs), \(f_{{}_{i}}\) and \(g_{{}_{i}}\), and their corresponding evolution equations [18,19]. As such, these are given here (Equations (13) and (14)) in the form of two steps, which are collision and streaming (advection) processes:  
\[f_{{}_{i}}\left(x+e_{{}_{i}},t+1\right)=f_{{}_{i}}\left(x,t\right)-\tau_{{}_{f }}^{-1}\left(f_{{}_{i}}\left(x,t\right)-f_{{}_{i}}^{eq}\left(x,t\right)\right) +\mathrm{Boussinesq\,\,force\,\,effect} \tag{13}\]  
\[g_{{}_{i}}\left(x+e_{{}_{i}},t+1\right)=g_{{}_{i}}\left(x,t\right)-\tau_{{}_{h }}^{-1}\left(g_{{}_{i}}\left(x,t\right)-g_{{}_{i}}^{eq}\left(x,t\right)\right) +S_{{}_{h}}^  
are, respectively, the weight coefficient and the velocity vector of the D2Q9 model; \(\overline{u}\left(u,v\right)\) is the macroscopic velocity, with \(u\) and \(v\) representing velocities in the \(x\) - and \(y\) -directions, respectively. Note that the relaxation times \(\tau_{f}\) and \(\tau_{h}\) can be determined via \(\
u_{lattice}=c_{s}^{2}\Delta t\left(\tau_{f}-0.5\right)\) and \(\alpha_{lattice}=c_{s}^{2}\Delta t\left(\tau_{h}-0.5\right),\ c_{s}\) being the lattice sound speed. It should be noted that, the lattice viscosity and diffusivity are selected so as to conform to the intended Prandtl number \(\Pr\left(=\
u_{lattice}/\alpha_{lattice}\right).\) Likewise, the additional force term related to Boussinesq force \(F_{b}\) is incorporated in the model by shifting the velocity field by a term of \(F_{b}\tau_{f}/\rho\), as proposed by [PERSON] and [PERSON] [20], where \(F_{b}\ \left(=-\beta g\left(T-T_{0}\right)/\left(T_{h}-T_{c}\right)\right).\) By this treatment, there is no need to add a force term to the collision operator. On the other hand, in the \(g\) distribution function, the source term is treated as per the method proposed by [PERSON] [21]. Hence, the resulting force in the LBM frame will be: \(S_{h}^{g}=-\omega_{i}S_{h}\) with \(S_{h}\left(=Ste^{-1}\partial\varepsilon/\partial t\right)\) being the source (or sink) term that handles the phase change.  
The liquid fraction \(\left(\varepsilon\right)\) is computed from [15; 16]:  
\[\varepsilon=\begin{cases}0&h<h_{s}=c_{p}T_{m}\\ \dfrac{h-h_{s}}{h_{l}-h_{s}}&h_{s}\leq h\leq h_{l}=h_{s}+L_{f}\\ 1&h>h_{l}\end{cases} \tag{18}\]  
where \(h\) is the local enthalpy defined, at a time step \(j\), by: \(h^{j}=c_{p}T^{j}+\left(\varepsilon^{j}-\varepsilon^{j-1}\right)L_{f}.\)  
It is useful to recall that in the current model, the solid and liquid phases are defined according to the liquid fraction value. Therefore, a mushy zone state is assigned when the value of \(\varepsilon\) is between zero and one. In this case, the velocity field is partially bounced back and the macroscopic velocity is modified as in Equation (19) [22]. The procedure of the implementation of the partial bounce-back approach is described in [23]. Note that, to implement BCs in the current method, we are led to convert them, at the mesoscopic level, in terms of the distribution function.