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.

Bayesian-based data assimilation methods integrate observational data into geophysical forward models to obtain the temporal evolution of an improved state vector, including its uncertainties. We explore the potential of a variant, a particle method, to estimate mechanical parameters of the overriding plate during the interseismic period. Here we assimilate vertical surface displacements into an elementary flexural model to estimate the elastic thickness of the overriding plate, and the locations and magnitudes of line loads acting on the overriding plate to produce flexure. Assimilation of synthetic observations sampled from a different forward model than is used in the particle method, reveal that synthetic seafloor data within 150 km from the trench are required to properly constrain parameters for long wavelength solutions of the upper plate (i.e. wavelength ∼500 km). Assimilation of synthetic observations sampled from the same flexural model used in the particle method shows remarkable convergence towards the true parameters with synthetic on-land data only for short to intermediate wavelength solutions (i.e. wavelengths between ∼100 and 300 km). In real-data assimilation experiments we assign representation errors due to discrepancies between our incorrect or incomplete physical model and the data. When assimilating continental data prior to the 2011 Mw Tohoku-Oki earthquake (1997-2000), an unrealistically low effective elastic plate thickness for Tohoku of ∼5-7 km is estimated. Our synthetic experiments suggest that improvements to the physical forward model, such as the inclusion of a slab, a megathrust interface and viscoelasticity of the mantle, including accurate seafloor data, and additional geodetic observations, may refine our estimates of the effective elastic plate thickness. Overall, we demonstrate the potential of using the particle method to constrain geodynamic parameters by providing constraints on parameters and corresponding uncertainty values. Using the particle method, we provide insights into the data network sensitivity and identify parameter trade-offs.