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.@en

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.@en

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.@en

Geodetic observations of vertical land motion following a megathrust earthquake are key to a better understanding of processes and parameters controlling the dynamics at subduction margins. The relative contributions of dominant drivers during the postseismic phase, such as viscoelastic relaxation, afterslip and relocking, remain difficult to estimate individually and are often derived at the end of an observation period. Data assimilation can provide the means to estimate model parameters by combining physical models with observations whilst taking their uncertainties into account. We use Bayesian inference in the form of an ensemble smoother to estimate geodynamic parameters during the postseismic phase of the megathrust earthquake cycle. The ensemble smoother uses a Monte Carlo approach to represent the probability density distribution (pdf) of model states with a finite number of realizations. Prior estimates of the imperfect physical model are combined with the likelihood of noisy observations to estimate the posterior pdf of model parameters. We discuss a synthetic data experiment where observations are sampled from a 3D finite element model with noise added to represent errors in the data. With a smoother, observations at all time steps are assimilated in one go, to ensure consistency between estimated parameters and model outputs. We assimilate vertical and horizontal surface displacements into a 2D finite element viscoelastic earthquake cycle model with a power-law rheology. We incorporate heterogeneity into the viscoelastic structure and estimate the viscosity field based on a heterogeneous temperature field and experimental flow laws. Preliminary results show that model parameters, such as the extent of the cold nose, maximum depth of afterslip and flow law parameters can be recovered remarkably well by assimilating synthetic on- and offshore surface observations. We will apply the smoother to postseismic surface displacements of the Tohoku 2011 earthquake to estimate model parameters and driving processes across spatio-temporal domains.@en

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.@en

Large parts of the Netherlands have experienced subsidence due to anthropogenic soil deformation (peat oxidation and clay shrinkage and swell) since the Middle Ages. As of the start of gas production from the Groningen reservoir in 1968, the northeastern Netherlands has seen additional subsidence (>30 cm) originating from deep reservoir compaction. Due the land being situated close to sea-level this poses a significant societal problem. Geodetic surface observations (InSAR) contain components of both shallow and deep subsidence. This study is part of a larger project on subsidence forecasting, which aims to disentangle deep and shallow (soil) drivers of subsidence by assimilating geodetic time series in geomechanical models of the subsurface. This requires the models to be highly efficient. In this study, we focus on modelling the deep subsurface (reservoir and overburden). We perform a sensitivity analysis to investigate the level of geometrical and mechanical complexity of the reservoir and overburden that can be resolved from the data, to reach optimal model efficiency. We employ a semi-analytical mechanical model for the Groningen subsurface using the PSGRN/PSCMP code by Wang et al. (2006) of laterally uniform viscoelastic layers. We use the Visvalingam–Whyatt algorithm to simplify the existing geological model to versions with decreasing levels of complexity. The surface deformation resulting from the different model versions are used in a convergence test to define the required model complexity. We find that we can achieve orders of magnitude improvement in model run time efficiency, depending on the data uncertainty. Wang, R., Lorenzo Martín, F., Roth, F. (2006): PSGRN/PSCMP - a new code for calculating co- and post-seismic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation theory. Computers and Geosciences, 32, 4, 527-541.@en

Plate boundary deformation zones represent a challenge in terms of understanding their underlying geodynamic drivers. Active deformation is well constrained by GNSS observations in the SW Balkans, Greece and W Turkey, and is characterized by variable extension and strike slip in an overall context of slow convergence of the Nubia plate relative to stable Eurasia. Diverse, and all potentially viable, forces and models have been proposed as the cause of the observed surface deformation, e.g., asthenospheric flow, horizontal gravitational stresses (HGSs) from lateral variations in gravitational potential energy, and rollback of the Hellenic slab. We use Bayesian inference to constrain the relative contribution of the proposed driving and resistive regional forces. Our models are spherical 2D finite element models representing vertical lithospheric averages. In addition to regional plate boundaries, the models include well-constrained fault zones like north and south branches of the North Anatolian Fault, Gulf of Corinth and faults bounding the Menderes Massif. Boundary conditions represent geodynamic processes: (1) far-field relative plate motions; (2) resistive fault tractions; (3) HGSs from lateral density variations; (4) slab pull and trench suction at subduction zones; and (5) active asthenospheric convection. The magnitude of each of these is a parameter in a Bayesian analysis of ~100,000 models and horizontal GNSS velocities. The search yields a probability distribution of all parameter values including model error, allowing us to determine mean/median parameter values, robustly estimate parameter uncertainties, and identify tradeoffs (i.e., parameter covariances). The average viscosity of the overriding plate is well resolved 3-4 10^22 Pa.s, which is higher than from published models without faults. Significant trench suction forces from the Hellenic slab act on the overriding Aegean Sea, including along the Pliny-Strabo STEP Fault. Slab pull and convective tractions have a small imprint on the observed deformation of the overriding plate. HGSs are necessary to explain local features in the velocity field, particularly in the Aegean Sea, but are less important for fitting the regional pattern of velocities. Resistive tractions on most plate boundaries and faults are low.@en

Most of the seismic moment release of the complex earthquake sequence beneath the South Sandwich Islands occurred on the central part of the SS megathrust. Significant aftershock activity indicates that the central and southern megathrust was subsequently activated, i.e., where young South America lithosphere is subducted. Seismic activity thus seems to have been restricted by the lateral termination in the south of the SS Trench. Relatively little energy release occurred on the northern part of the megathrust. It was hypothesized by Govers and Wortel (2005) that here the South America slab breaks away from the surface part of the plate at the active STEP. Geochemical observations and earthquake P-axes orientations do not seem to agree with the hypothesis and we investigate the cause. We show results of new physical analog lab models that aim to elucidate what controls the geometry of the lithospheric STEP Fault. We study lithospheric tearing in the process of STEP evolution, which is dynamically driven by the buoyancy of the subducting slab. In our experiments, the lithosphere as well as asthenosphere are viscoelastic media in a free subduction setup. A stress-dependent rheology plays a major role in localization of strain in tearing processes of lithosphere such as slab break-off. The results show that the highly curved northern plate boundary is a STEP Fault following from lithospheric tearing at a depth of ~100km. This is a modification of the original STEP model of Govers and Wortel (2005). This is consistent with available observations along the northern Sandwich plate boundary, and likely exists in other STEP regions. The region’s largest recorded event, the 1929 Mw 8.3 earthquake, may reflect horizontal extension perpendicular to the STEP fault, which is also expected based on our experiments.@en

Stress-dependent nonlinear upper mantle rheology has a firm base in rock mechanical tests, where this nonlinearity results from dislocation creep of minerals. In the last few decades there has been some attention to nonlinear, power-law, materials for application in scaled analogue experiments for tectonic processes. However, studies describing the rheology of analogue materials with the same nonlinear dependency on stress as observed for lithospheric mantle materials at relevant stress levels, are still lacking. In this study we have developed and rheologically tested materials based on combinations of silicone polymers and plasticine, with the aim of obtaining a material that can serve as a laboratory analogue to the power-law rheology of olivine aggregates at lithospheric mantle conditions. From our steady-state creep tests we find that it is possible to obtain such a power-law material, with effective viscosities over relevant model stress ranges [5–4000 Pa] that allow for nonlinear deformation at laboratory time scales. We apply the developed material to a process where localized deformation of the lithosphere can be expected: slab break-off. We study this process using analogue models, where we apply the new nonlinear material to the lithospheric mantle domains, while we use Newtonian glucose to represent the low viscous asthenosphere. Now that we properly manage power-law behavior in our analogue lithosphere materials, we are able to model localized lithospheric tearing.

@en

The Gibraltar Arc (GA; Western Mediterranean) is located in between the slowly converging African and Eurasian continents and runs from the Rif in northern Morocco to the Betics in southern Spain. Geodetic velocities in this wide plate boundary region are different from the motions of both stable Eurasia and Africa. Various lithospheric and asthenospheric drivers have been called upon to explain the observed surface motions. We constrain the magnitude of these driving forces and the effective mechanical properties of the lithosphere and fault zones using regional geodynamic Finite Element models (FEMs). All of our FEMs have lithospheric forces that arise from lateral variations of gravitational potential energy, mainly from lateral variations in topography and Moho topology. They also all include far-field Africa-Eurasia convergence. The well-imaged Gibraltar slab suggests that slab pull and/or slab suction may contribute to the surface motions. Also, tractions that arise from lateral transport of the slab through the upper mantle are potentially relevant. We adopt parameterized versions of these slab forces in the FEMs. In our FEMs, we impose shear tractions on observationally constrained fault zones. We quantify the model parameter resolution and tradeoffs in light of the available kinematic surface observations with a Markov Chain Monte Carlo (MCMC) approach, implementing the Metropolis-Hastings algorithm. We compare the model results to the geodetic velocity observations, which have the smallest error margins of all the kinematic observations, and to the sense of shear along (potentially) active fault zones. Preliminary results based on a search of 60,000 models indicate that the observed GNSS velocity field and sense of slip on regional faults in the Gibraltar Arc appear to result mainly from Africa-Europe plate convergence and lateral GPE variations. Slab pull from the Gibraltar slab is very likely transmitted poorly into the overriding plate and probability distributions for the trench suction force do not display any favorable value to affect the kinematics in the Gibraltar Arc region. The best models have rms-misfits of 0.23 mm/yr, which is mostly due to a systematic SW motion of up to 4 mm/yr in the SW Rif.@en