We introduced a novel approach to produce a 3D heterogeneous viscosity distribution model of the upper mantle below Northeast Japan. The viscoelastic relaxation constitutes a significant portion of the postseismic deformation associated with a megathrust earthquake. Our viscosity distribution was based on a 3D simulation of thermomechanical modeling with realistic 3D oceanic plate geometry and subduction velocity of the Pacific plate in the upper mantle. After obtaining the thermal and strain-rate distributions obtained from the thermomechanical modeling, 3D viscosity distribution was calculated. We used the resulting viscosity distribution in a viscoelastic relaxation model to calculate the viscoelastic response associated with the MW 9.0 2011 Tohoku-Oki earthquake, using a finite element method. We found it necessary to include an elastic cold nose in the mantle wedge and a low-viscosity layer at the subducted lithosphere-asthenosphere boundary to reproduce the observed onshore uplifts and landward horizontal displacements of the seafloor, respectively. We then divided the calculated viscoelastic displacements into 7 one-year intervals and subtracted them from the corresponding observed postseismic displacements to obtain yearly residual displacements. Assuming the residual displacements were the result of afterslip and subduction fault locking on the plate interface, we employed a time-dependent inversion to estimate the spatiotemporal distributions of afterslip and subduction fault locking. Our results showed decrease in the afterslip amount and extent over time, which primarily occurred during the first two years following the 2011 event, which is followed by gradual extension of the subduction fault locking region. We concluded that the plate interface off Northeast Japan is gradually entering an interseismic stage of interplate coupling and elastic strain accumulation within several years after the occurrence of the Tohoku-Oki earthquake.