## 1 Introduction
Eventually, the initial FE mesh will become too distorted and a new mesh must be generated in the solid and liquid regions. In the implicit methods, on the other hand, the computational grid is usually fixed and the position of the interface is obtained indirectly from some field defined over the whole domain (and suitably discretised). Examples of implicit methods are the enthalpy method, the phase field method, and the level set method. In the _enthalpy method_ (e.g. [8, 9]), the enthalpy field is considered and the position of the interface is determined in a smeared way by a jump of such field. In the _phase field method_ (e.g. [10]) a phase function is defined over the domain; it assumes fixed values in each of the phases (e.g. 1 in solid and -1 in liquid) and varies smoothly between these values in the interface region. In order to allow such smooth variation, an artificial 'interface thickness' must be added to the model. Hence, the interface position is not defined exactly, although it can be chosen conventionally as the surface where the phase field assumes an average value (e.g. zero) or as some kind of 'average surface' of the interface region; hence, a sufficiently fine mesh is required in order to resolve the interface zone. Of course, it is necessary to provide additional governing equation for the phase field evolution in a thermodynamically-consistent way. Finally, in the _level set method_ a function is defined all over the domain and the position of the phase interface is obtained as a level set of this function. In comparing the last two methods, we see that, although the level set function is artificial, the physics of the problem is fully respected; on the other hand, the phase field has a physical meaning, but it necessitates of an artificial 'thickness' of the interface. Some comparative analyses in Reference [11] lead to the conclusion that methods based on the level set are the most promising and general. For these reasons, and since it allows a general, elegant formulation of the extended finite element method, the level set method is chosen here for the numerical analyses and it is reviewed in Section 3. The second issue in the numerical treatment of the Stefan problem is the discretisation of the temperature field. The crucial feature, in order to account for the evolution of the phase interface, is the discontinuity in the temperature gradient across the interface itself. In the aforementioned enthalpy method the problem is reformulated in terms of enthalpy instead of temperature. Although this approach has the advantage of offering a straightforward discretisation of the problem, it has two main disadvantages, as the discrete set of equation in terms of enthalpy becomes nonlinear and the phase interface is smeared throughout a finite thickness. In the so-called _boundary immobilisation method_, a new set of coordinates is defined in which the phase-boundary is fixed. [PERSON] and [PERSON], in their comparative analysis [12], highlight the effectiveness of this method for one-dimensional problems, but discard it for planar or spatial cases because of the difficulties in defining the coordinate transform with respect to a reference coordinate system fixed in time. Other methods, which are not discussed in depth here due to their lack of generality, are the _perturbation method_, which requires hints from a partial analytical solution, the _nodal integral method_, which only applies to plane symmetric geometries, and the _heat balance integral method_, which is restricted to time-independent boundary conditions. For details on these methods, we remand to the last quoted reference and further references therein. More general approaches to the temperature field discretisation are based on standard techniques, such as the _finite difference method_ with a moving grid (e.g. [11]), the _boundary element method_, and, above all, the _finite element method_ with adaptive mesh (for references on the last two methods applied to Stefan problem, see in Reference [13]).