The numerical modelling of plate subduction requires solving a coupled thermo-mechanical highly-nonlinear transient problem. The mechanical description of the phenomenon results in a multiphase quasi-static Stokes flow, ...