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. ... Read more

Before geodetically derived strain and rotation rates can be robustly compared to geological or seismological observations, we need reliable strain rate uncertainties. Various methods exist to compute strain rates from GNSS-derived interseismic velocities, but a realistic representation of interpolation uncertainties has remained a challenge. The main problem is that commonly used deterministic interpolation methods do not account for uncertainty resulting from the absence of information in between observation sites. We apply stochastic interpolation by means of ordinary kriging to propagate errors both from discontinuous data coverage as well as from observation uncertainties to our strain rate estimates. However, interseismic horizontal surface velocities in tectonically active regions are spatially highly non-stationary, with high spatial variability around active faults and lower velocity variability in tectonically more stable regions. This requires an extension of traditional ordinary kriging approaches. For interpolation uncertainties that reflect the local variability and spatial correlation of the observed surface velocities, we apply a novel method that incorporates the spatially variable statistics of the underlying data. We estimate realistic uncertainties and covariances of the interpolated velocity field. For regions with a high spatial velocity variability, we find a large increase in uncertainty with increasing distance from observation sites, while in areas with little spatial variability, we estimate a small increase in uncertainty with distance. Subsequently, we propagate interpolated velocity covariance to strain rate uncertainties, such that we can assess the statistical significance of the interpolated strain rate field. Applied to a number of actively deforming regions, including the Pacific coast of North America and Japan, we show to what degree we can robustly determine strain rates based on available GNSS-derived velocities. Realistic uncertainties assist the community to better discriminate continuous or localized deformation on active faults from the available geodetic data. ... Read more

We aim to better understand the overriding plate deformation during the megathrust earthquake cycle. We estimate the spatial patterns of interseismic GNSS velocities in South America, Southeast Asia and northern Japan and the associated uncertainties due to variations in network density and observation uncertainties. Interseismic velocities with respect to the overriding plate generally decrease with distance from the trench with a steep gradient up to a ‘hurdle’, beyond which the gradient is distinctly lower and velocities are small. The hurdle is located 500–1000 km away from the trench for the trench-perpendicular velocity component, and either at the same distance or closer for the trench-parallel component. Significant coseismic displacements were observed beyond these hurdles during the 2010 Maule, 2004 Sumatra–Andaman, and 2011 Tohoku earthquakes. We hypothesize that both the interseismic hurdle and the coseismic response result from a mechanical contrast in the overriding plate. We test our hypothesis using physically consistent, generic, 3-D finite element models of the earthquake cycle. Our models show a response similar to the interseismic and coseismic observations for a compliant near-trench overriding plate and an at least five times stiffer overriding plate beyond the contrast. The model results suggest that hurdles are more prominently expressed in observations near strongly locked megathrusts. Previous studies inferred major tectonic or geological boundaries and seismological contrasts located close to the observed hurdles in the studied overriding plates. The compliance contrast probably results from thermal, compositional and thickness contrasts and might cause the observed focusing of smaller-scale deformation like backthrusting. ... Read more

Greater landward velocities were recorded after 6 megathrust earthquakes in subduction zone regions adjacent to the ruptured portion. Previous explanations invoked either increased slip deficit accumulation or plate bending during postseismic relaxation, with different implications for seismic hazard. We investigate whether bending can be expected to reproduce this observed enhanced landward motion (ELM). We use 3D quasi-dynamic finite element models with periodic earthquakes. We find that afterslip downdip of the brittle megathrust exclusively produces enhanced trenchward surface motion in the overriding plate. Viscous relaxation produces ELM when a depth limit is imposed on afterslip. This landward motion results primarily from in-plane elastic bending of the overriding plate due to trenchward viscous flow in the mantle wedge near the rupture. Modeled ELM is, however, incompatible with the observations, which are an order of magnitude greater and last longer after the earthquake. Varying mantle viscosity, plate elasticity, maximum afterslip depth, earthquake size, and megathrust locking outside of the rupture does not significantly change this conclusion. The observed ELM consequently appears to reflect faster slip deficit accumulation, implying a greater seismic hazard in lateral segments of the subduction zone. ... Read more