[1] The underestimation of the size of recent megathrust earthquakes illustrates our limited understanding of their spatiotemporal occurrence and governing physics. To unravel their relation to associated subduction dynamics and long‐term deformation, we developed a 2‐D continuum viscoelastoplastic model that uses an Eulerian‐Lagrangian finite difference framework with similar on‐ and off‐fault physics. We extend the validation of this numerical tool to a realistic subduction zone setting that resembles Southern Chile. The resulting quasi‐periodic pattern of quasi‐characteristic M8–M9 megathrust events compares quantitatively with observed recurrence and earthquake source parameters, albeit at very slow coseismic speeds. Without any data fitting, surface displacements agree with GPS data recorded before and during the 2010 M8.8 Maule earthquake, including the presence of a second‐order flexural bulge. These surface displacements show cycle‐to‐cycle variations of slip deficits, which overall accommodate ∼5% of permanent internal shortening. We find that thermally (and stress) driven creep governs a spontaneous conditionally stable downdip transition zone between temperatures of ∼350°C and ∼450°C. Ruptures initiate above it (and below the forearc Moho), propagate within it, interspersed by small intermittent events, and arrest below it as ductile shearing relaxes stresses. Ruptures typically propagate upward along lithological boundaries and widen as pressures drop. The main thrust is constrained to be weak due to fluid‐induced weakening required to sustain regular subduction and to generate events with natural characteristics (fluid pressures of ∼75–99% of solid pressures). The agreement with a range of seismological, geodetic, and geological observations demonstrates the validity and strength of this physically consistent seismo‐thermo‐mechanical approach.