Geodetic data spanning different phases of the earthquake cycle offer insights into the spatiotemporal interplay between processes driving surface deformation, such as viscoelastic relaxation, afterslip, and (re)locking. However, quantifying their contributions and explaining pre- and post-earthquake displacements with a single set of rheological parameters is challenging. We set up a 2-D earthquake cycle finite element model that simulates the mantle and a thin low-viscosity shear zone with a temperature-dependent linear Maxwell or nonlinear power-law rheology. We use the ensemble smoother with multiple data assimilation to estimate ensembles of parameters describing the rheological makeup of the subduction zone. We assimilate onshore and offshore displacement time series acquired before and after the 2011 Tohoku-Oki earthquake. Our models provide a unique, robust solution using a temperature-dependent power-law rheology. The estimated creep parameters for the mantle wedge deeper than ∼50 km and sub-slab mantle align with laboratory experiments. However, different creep parameters are required for the shallow part of the mantle wedge than the deeper part to explain the observed postseismic response—highlighting the need for shallow viscoelastic relaxation. The trade-off between water fugacity and activation energy hinders their individual estimation but yields a well-constrained viscosity structure. The spatial distribution of vertical displacements as well as the temporal signature of early postseismic horizontal displacements are required to estimate individual parameters for afterslip and viscoelastic relaxation. Afterslip occurs downdip of the coseismic rupture. Near-trench landward motion during the early postseismic period is driven by elastic stress release beneath the oceanic plate and sub-slab asthenospheric flow.