Numerical methods to solve problems involving discontinuities (jumps, kinks or singularities) on moving internal boundaries have received much attention over the last decade. Among them, the most suitable is probably the extended finite element method (XFEM) in tandem with the level-set technique due to its ability to take into account these discontinuities without matching meshes [1]. The present contribution aims to elaborate a numerical approach to model interfacial discontinu- ities within a meshless context. This approach couples the constrained natural element method (C-NEM) [2] and the level-set technique. In the former, the natural neighbours interpolation, based on a Voronoi diagram, is locally enriched through the partition of unity concept. This enrichment is built from level-set functions that represent and track implicitly discontinuities inside the domain [3]. Like in XFEM, key features of the proposed approach is (i) to determine the intersection between Voronoi cells and discontinuities and (ii) to integrate numerically the weak form over cells containing discontinuities. After testing the proposed method on classical benchmarks, both accuracy and efficiency are examined on the two phase Stefan problem that deals with heat flow involving a solid-liquid phase boundary on which a jump condition must be satisfied [4]. [1] Chessa J., Smolinski P. and Belytschko T. The extended finite element method (XFEM) for solidification problems. Int. J. Numer. Meth. Engng 53:1959–1977 (2002). [2] Yvonnet J., Chinesta F., Lorong P. and Ryckelynck D. The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces. Int. J. Therm. Sci. 44:559–569 (2005). [3] LiuJ.T.,GuS.T.,MonteiroE.andHeQ.C.Aversatileinterfacemodelforthermalconduc- tion phenomena and its numerical implementation by XFEM. Comp. Mech. 53:825–843 (2014). [4] Carslaw H.S. and Jaeger J.C. Conduction of Heat in Solids. 2th Edition, Clarendon Press, (1959)