Subduction megathrusts develop the largest earthquakes, often close to large population centers. Understanding the dynamics of deformation at subduction zones is therefore important to better assess seismic hazards. Here I develop consistent earthquake cycle simulations that incorporate localized and distributed deformation based on laboratory-derived constitutive laws by combining boundary and volume elements to represent the mechanical coupling between megathrust slip and solid-state flow in the oceanic asthenosphere and in the mantle wedge. The model is simplified, in two dimensions, but may help the interpretation of geodetic data. Megathrust earthquakes and slow-slip events modulate the strain rate in the upper mantle, leading to large variations of effective viscosity in space and time and a complex pattern of surface deformation. While fault slip and flow in the mantle wedge generate surface displacements in the same, that is, seaward, direction, the viscoelastic relaxation in the oceanic asthenosphere generates transient surface deformation in the opposite, that is, landward, direction above the rupture area of the mainshock. Aseismic deformation above the seismogenic zone may be challenging to record, but it may reveal important constraints about the rheology of the subducting plate. Plain Language Summary Mitigation and preparation for seismic hazards rely on numerical simulation of earthquakes grounded in laboratory experiments and field data. Surprisingly, bringing together theory and predictions is still remarkably challenging. In this study, I demonstrate how recent progress in numerical methods now affords an efficient way to simulate how earthquakes come about at subduction zones, including how other regions of the Earth deform during the seismic cycle. The model predicts unexpected results about the Earth's surface deformation following large earthquakes.