We seek to quantify bulk viscoelastic flow, afterslip, and locking, within a rheological framework that is consistent over the entire earthquake cycle. We address this using an ensemble smoother. We construct a 2D finite element seismic cycle model with a power-law rheology in the asthenosphere. A priori information, such as a realistic temperature field and a coseismic slip distribution, is integrated into the model. Model pre-stresses are initialized during repeated earthquake cycles wherein the accumulated slip deficit is released entirely. We tailor the last earthquake to match the observed co-seismic slip of the 2011 Tohoku earthquake. The heterogeneous rheology structure is derived from the temperature field and experimental flow laws. Additionally, we simulate afterslip using a thin viscoelastic shear zone. We focus on constraining power-law flow parameters for the asthenosphere and the shear zone. We assimilate 3D GEONET GNSS displacement time series acquired before and after the 2011 Tohoku earthquake. Power-law viscosity parameters are successfully retrieved for all domains. The data require separate viscoelastic domains in the mantle wedge above and below ~50 km depth. The sub-slab asthenosphere has viscoelastic properties that are distinctly different from the mantle wedge. The trade-off between the power-law activation energy and water fugacity hinders their individual estimation. The wedge viscosity is >10^19 Pa·s during the interseismic phase. Postseismic afterslip and bulk viscoelastic relaxation can be individually resolved from the surface deformation data. Afterslip is substantial between 40-50 km depth and extends to 80 km depth. Bulk viscoelastic relaxation in the wedge concentrates above 150 km depth with viscosities <10^18 Pa·s. Landward motion of the near-trench region occurs during the early postseismic period without the need for a separate low-viscosity channel below the slab.