Understanding the buildup of stress in space and time at tectonic plate interfaces is a key component of explaining sequences of earthquakes and aseismic slip (SEAS) at plate boundaries. In turn, successful modeling of coseismic slip and post- and interseismic fault creep can provide insights into mechanical properties of the fault zone material, as well as into potential seismic hazards. In the past, incomplete modeling of geodetic data to infer plate interface coupling has led to false assumptions of relative seismic quiescence. The issue of stress shadows in particular has led to the realization that improving our knowledge of the stress deficit on any given plate interface requires modeling the influence of past earthquakes (e.g., Hetland & Simons, 2010).This study aims to model the known historical sequence of earthquakes along the Northern Japan subduction zone in terms of post- and interseismic fault creep, at least to the extent it affects observed GNSS surface displacement timeseries of the last decades. We start from the framework presented by Kanda et al. (2013), who performed endmember simulations of multi-cycle SEAS on the Japanese megathrust, assuming varying rheological properties and fully-locked patches on the interface ("asperities") in order to match GNSS-derived interseismic plate velocities. We extend this model by (1) taking into account the entire GNSS network timeseries available in Northern Japan (containing co-, post-, and interseismic periods) and (2) integrating it into a Markov chain Monte Carlo (MCMC) solver to solve for the most probable rheological properties of the plate interface.We present the results of selected forward simulations, including timeseries of slip rate on the fault and surface displacement, and compare them to GNSS observations. We furthermore present preliminary results of the inferred best-fitting rheological parameters (assuming rate-dependent friction).