During melting under gravity in the presence of a horizontal thermal gradient, buoyancy-driven convection in the liquid phase affects significantly the evolution of the liquid-solid interface. Due to the obvious engineering interest in understanding and controlling melting processes, fluid dynamicists and applied mathematicians have spent many efforts to model and simulate them numerically. Their endeavors concentrated in the twenty-five years period between the publication of the paper by Brent, Voller & Reid (1988) and that by Mansutti & Bucchignani (2011).

In this paper, an algorithm for the numerical evaluation of hypersingular finite-part integrals with rapidly oscillating kernels is proposed. The method is based on an interpolatory procedure at zeros of the orthogonal polynomials with respect to the first kind Chebyshev weight. Bounds of the error and of the amplification factor are also provided. Numerically stable procedure are obtained and the corresponding algorithms can be implemented in a fast way.

Recent literature confirms the crucial influence of non-viscous deformations together with temperature impact on glacier and rock glacier flow numerical simulation. Along this line, supported by the successful test on a one-dimensional set-up developed by two of the author, we propose the numerical solution of a two-dimensional rock-glacier flow model based on an ice constitutive law of second grade differential type .