## 2 Governing Equations  
The _Stefan problem_ of melting/solidification in absence of material transport is considered. For the case of the melting/solidification of a two-phase flow, see Reference [20]. Let \(\Omega\) be the domain shown in Figure 1 and \(\partial\Omega\) its boundary with outward normal unit vector \(\mathbf{n}\). The domain is composed of a solid region \(\Omega_{\text{sol}}\) and a liquid region \(\Omega_{\text{liq}}=\Omega\setminus\overline{\Omega}_{\text{sol}}\). We denote by \(\Gamma_{\text{tran}}\) the interface between the two regions and by \(\mathbf{n}_{\text{sol}}\) the outward normal from \(\Omega_{\text{sol}}\). The temperature in the point \(\mathbf{x}\in\Omega\) at time \(t\in\left[0,\overline{t}\right]\) is \(T\left(\mathbf{x},t\right)\). The model accounts for the heat transferred by conduction only; therefore, convection, thermal expansion, and buoyancy are not considered. Since deformations are excluded, the mass density \(\rho\) is the same in the solid and liquid phases. The energy conservation equation in the domain reads  
\[\rho\frac{\partial}{\partial t}\left(cT\right)=-\operatorname{div}\mathbf{q} +s\, \tag{1}\]  
where \(c\) is the heat capacity, \(\mathbf{q}\) the heat flux vector, and \(s\) the volumetric heat source. The constitutive equation which characterises the conductivity within the material is given by Fourier law,  
\[\mathbf{q}=-\mathbf{k}\
abla T\, \tag{2}\]  
where \(\mathbf{k}\) is the thermal conductivity second-order tensor. In the following, we assume isotropy with respect to heat conduction, so that \(\mathbf{k}=k\mathbf{1}\), where \(\mathbf{1}\) denotes the identity tensor and \(k\) the scalar thermal conductivity. \(k\) and \(c\) are taken as constants in each phase and their values in the solid and in the liquid are indicated by subscripts'sol' and 'liq' respectively. The essential and natural boundary conditions are  
Figure 1: Schematic representation of the domain. \[T=\hat{T}\quad\text{ on }\Gamma_{T}\,, \tag{3}\] \[-\mathbf{k}\
abla T\cdot\mathbf{n}=\hat{q}\quad\text{ on }\Gamma_{q}\,,\]  
where \(\hat{T}\big{(}\mathbf{x},t\big{)}\) and \(\hat{q}\big{(}\mathbf{x},t\big{)}\) are assigned functions defined over \(\Gamma_{T}\times\big{[}0,\overline{t}\big{]}\) and \(\Gamma_{q}\times\big{[}0,\overline{t}\big{]}\) respectively, with \(\Gamma_{T}\cup\Gamma_{q}=\partial\Omega\) and \(\Gamma_{T}\cap\Gamma_{q}=\emptyset\). The phase-interface is characterised by the constant temperature \(T_{\text{\tiny tran}}\), so that the condition  
\[T=T_{\text{\tiny tran}}\quad\text{ on }\Gamma_{\text{\tiny tran}} \tag{4}\]  
holds. It is noted that the temperature field must be continuous through the interface; this requirement must be included (or enforced) in the numerical models as discussed in Section 4. It is also possible to consider more complicate interface conditions which involve the velocity of the interface or its geometrical properties. An example is Gibbs-Thomson relation, which describes the unstable dendritic growth of crystals into an undercooled melt, and modifies Eq. (4) as \(T=T_{\text{\tiny tran}}+\gamma_{{}_{\kappa}}\kappa+\gamma_{{}_{\
u}}v\), where \(\gamma_{{}_{\kappa}}\) and \(\gamma_{{}_{\
u}}\) are the surface tension and the kinetic mobility coefficients, \(\kappa=\operatorname{div}\mathbf{n}_{\text{\tiny sol}}\) is the mean curvature of the interface, and \(v\) its normal speed (negative for melting and positive for solidification).