Building an advanced numerical model of subduction requires choosing values for various geometrical parameters and material properties, among others, the initial lithosphere thicknesses, representative lithological types and their mechanical and thermal properties, rheologies, initial temperature profiles as well as model boundary conditions. Some of these can be constrained by observations of present-day subduction zones, such as lithosphere age and thickness, whereas others have a larger uncertainty and require critical evaluation of geophysical and experimental results. To test the model response to systematic variations in input parameters, numerical studies often start from a 'reference' subduction model. However, the reference model often varies between different numerical studies, making it difficult to compare results directly. We aim therefore to define a numerical reference model for thermomechanical subduction. This reference setup will facilitate comparisons of a series of numerical models that focus on different aspects of subduction, such as the effects of elasticity on the stress distribution, the energetic impact of phase transformations or the influence of devolatilization reactions. Our reference model represents ocean-ocean convergence and describes initial geometries and lithological stratification for a three-layered subducting slab and overriding plate along with their respective flow laws and chemical composition. It also includes kinematic and thermal boundary conditions, and initial temperature distribution. We will show results of the evolution and dynamics of the subduction reference model using different numerical codes: a finite element code, SULEC, and two finite difference codes, YACC and FDcon.