## 1 Introduction  
In its original formulation, the _Stefan problem_ is a representation of (temperature-driven) first-order phase transitions in matter through a boundary value problem for PDEs in which a discontinuity surface, internal to the domain, can move with time. Indeed, the problem was originally formulated in relation to ice formation and melting. However, the same mathematical framework can be applied to a variety of problems arising in different areas of science and technology, such as the diffusion of gasses in biologic tissues, biofilm and hydrogel growth [1; 2], the penetration of solvents in polymers, the flow in porous media, filtration problems, free surface flows, etching, shock propagation, and financial mathematics (e.g. [3]). For the formulation details of some of the physical problems mentioned above, and for further examples, see Reference [4]. Since the topic of this paper is the numerical treatment of the problem, we select a specific physical case, that is the solid-liquid phase transition (e.g. [5]). Applicative examples are the casting of metals, the freezing and thawing of foods, the production of ice, or ice formation on pipes. As in all the first-order phase transitions, a discontinuity in the first derivative of the free energy with respect to a thermodynamic variable is exhibited and a latent heat is involved. In a solid-liquid transition at constant volume, the pressure, i.e. the first derivative of the free energy with respect to the specific volume, is discontinuous. During the transition, the system either absorbs or releases a fixed amount of energy at a constant temperature. Since the heat cannot be exchanged instantly, a front of solidification or melting is present and moves according to the heat transport across the front itself. From the mathematical point of view, the solid and liquid regions are subdomains in which the coefficients (representing the physical properties of the medium in each phase) of the underlying PDEs are continuous and differentiable up to the order of the PDEs. The coefficients are discontinuous across the surfaces that separate the adjacent phases and the PDEs are not valid there, so that additional equations are needed for closure. These are derived from energy conservation, namely by the Stefan condition that expresses the local velocity of the moving freezing/melting front as a function of the heat flux evaluated at both sides of the phase boundary. Therefore, the solution of the heat-conduction equations are required within unknown subdomains (the liquid and the solid regions) and their interface must be determined as part of the solution. Existence and uniqueness of the solution to Stefan problem are demonstrated in References [6] and [7], respectively. However, analytical solutions are available in a close form only for a restricted number of simple, particular cases, all characterised by a high degree of symmetry in the geometry and in the boundary and initial conditions. For all the other cases a numerical treatment is required. For developing a numerical scheme for Stefan problem, two issues must be solved. Firstly, the evolving geometry of the phase-interface must be suitably described with a discrete model; secondly, the temperature field must be discretised in such a way as the jump in the heat flux, i.e. in the temperature gradient, at the phase interface can be reproduced. The description of the phase-interface can be done explicitly or implicitly. In the explicit methods, a Lagrangian approach is followed and some marker points (finite in number) on the phase interface are selected and tracked explicitly. This is the case, for instance, in which the temperature field is discretised by standard Finite Element Method (FEM); the mesh is generated at the initial time instant in such a way as the element boundaries lay on the phase interface. The finite element nodes on the interface are used as marker points and, at each time step, are displaced according to the solution of the heat-transfer problem within each phase.