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.