We study the evolution of a melting front between the solid and liquid phases of a pure incompressible material where fluid motions are driven by unstable temperature gradients. In a plane-layer geometry, this can be seen as classical Rayleigh–Bénard convection where the upper solid boundary is allowed to melt due to the heat flux brought by the fluid underneath. This free-boundary problem is studied numerically in two dimensions using a phase-field approach, classically used to study the melting and solidification of alloys, which we dynamically couple with the Navier–Stokes equations in the Boussinesq approximation. The advantage of this approach is that it requires only moderate modifications of classical numerical methods. We focus on the case where the solid is initially nearly isothermal, so that the evolution of the topography is related to the inhomogeneous heat flux from thermal convection, and does not depend on the conduction problem in the solid. From a very thin stable layer of fluid, convection cells appear as the depth – and therefore the effective Rayleigh number – of the layer increases. The continuous melting of the solid leads to dynamical transitions between different convection cell sizes and topography amplitudes. The Nusselt number can be larger than its value for a planar upper boundary, due to the feedback of the topography on the flow, which can stabilize large-scale laminar convection cells.