Description: Numerically modelling the dynamics of a self-consistently subducting lithosphere is a challenging task because of the decoupling problems of the slab from the free surface. We address this problem with a benchmark comparison between various numerical codes (Eulerian and Lagrangian, Finite Element and Finite Difference, with and without markers) as well as a laboratory experiment. The benchmark test consists of three prescribed setups of viscous flow, driven by compositional buoyancy, and with a low viscosity, zero-density top layer to approximate a free surface. Alternatively, a fully free surface is assumed. Our results with a weak top layer indicate that the convergence of the subduction behaviour with increasing resolution strongly depends on the averaging scheme for viscosity near moving rheological boundaries. Harmonic means result in fastest subduction, arithmetic means produces slow subduction and geometric mean results in intermediate behaviour. A few cases with the infinite norm scheme have been tested and result in convergence behaviour between that of arithmetic and geometric averaging. Satisfactory convergence of results is only reached in one case with a very strong slab, while for the other cases complete convergence appears mostly beyond presently feasible grid resolution. Analysing the behaviour of the weak zero-density top layer reveals that this problem is caused by the entrainment of the weak material into a lubrication layer on top of the subducting slab whose thickness turns out to be smaller than even the finest grid resolution. Agreement between the free surface runs and the weak top layer models is satisfactory only if both approaches use high resolution. Comparison of numerical models with a free surface laboratory experiment shows that (1) Lagrangian-based free surface numerical models can closely reproduce the laboratory experiments provided that sufficient numerical resolution is employed and (2) Eulerian-based codes with a weak surface layer reproduce the experiment if harmonic or geometric averaging of viscosity is used. The harmonic mean is also preferred if circular high viscosity bodies with or without a lubrication layer are considered. We conclude that modelling the free surface of subduction by a weak zero-density layer gives good results for highest resolutions, but otherwise care has to be taken in (1) handling the associated entrainment and formation of a lubrication layer and (2) choosing the appropriate averaging scheme for viscosity at rheological boundaries.