We process the SuGAr raw data using the GPS-Inferred Positioning System and Orbit Analysis Simulation Software (GIPSY-OASIS) software (hosted at https://gipsy-oasis.jpl.nasa.gov/). To extract the postseismic transient due to the 2012 M 8.6 Indian Ocean earthquake, we identify and remove the signals from other earthquakes in the entire available time series for each component of each station. In this procedure, we simultaneously estimate the linear long-term rates, coseismic, postseismic, and seasonal signals by nonlinear least-squares fitting. A complete description of these steps can be found in refs 6 and 31. We isolate the postseismic deformation from the 2012 Indian Ocean earthquake sequence by removing the contribution of all other identified sources. The resulting time series is shown in Extended Data Fig. 1.
The effect of the transient creep of rocks, whether for diffusion creep or dislocation creep, is well known in the laboratory and some studies suggest that it is important to incorporate it in models of postseismic deformation18, 19, 20 or postglacial rebound32. Several attempts have been made to formulate a rheological law for transient creep33, 34, 35, 36, 37, but they have important shortcomings. For example, laws such as33, 34
where A, q and r are positive constitutive parameters, or35, 36
where B, ε , r and s are positive constants, can create only a single transient creep episode, as plastic strain systematically increases, and therefore they do not describe transient creep for repeated earthquake stress perturbations. Some other laws, such as the so-called Andrade creep37
where is the steady-state strain rate and A and B are positive parameters, are singular38 at t = 0, for n > 1. More fundamentally, laws such as equation (3) fail to satisfy fundamental principles such as time-frame invariance (the same flow should be predicted for any clock).
There are several ways to formulate a rheology for the transient creep of olivine by invoking work hardening. For simplicity, we propose a generalization of the Burgers rheology. The Burgers rheology may appropriately represent transient creep in the diffusion creep regime, and we adapt it to be compatible with the power-law steady state of dislocation creep.
In a Burgers material, where the Kelvin element and the Maxwell element are in series, the inelastic strain rate can be written
where is the inelastic strain rate in the Kelvin element and is the inelastic strain rate in the Maxwell element. The strain rate of the Maxwell element is given by
where the parameters are defined in the main text. In general, we can then formulate a generalized rheology for the Kelvin element of the form
where σ is the deviatoric stress (the same as in the Maxwell element in series) and is the stress in the Kelvin dashpot. Because the functional form f is unknown, we assume for simplicity that it is the same as for the steady-state creep
where G is a work hardening coefficient. In the absence of detailed laboratory data on wet olivine, we assume that the thermodynamic parameters for transient and steady-state creep are the same (for wet olivine, ref. 16 shows that flow law parameters are similar between transient and steady-state creep except for the pre-exponential factor). However, laboratory experiments indicate that if transient creep is due to the transition of slip from the soft slip system to the hard slip system of olivine, transient creep should be less sensitive to the water content22.
At the background stress σ , if the flow is at steady state (that is, = 0), we have where is the cumulative strain at the background stress. After a stress perturbation from an earthquake, the stress changes from and the rate of transient creep instantaneously becomes . So transient creep, unlike steady-state creep, is not sensitive to the background stress during the postseismic transient, unless it was not at steady state. Multiple stress perturbations also lead to multiple transient creep episodes.
We explore the predictions of the proposed rheology in a spring-slider system. Together with the constitutive relationship, conservation of momentum leads to the system of coupled ordinary differential equations
We solve these equations numerically with unit stress perturbations and typical laboratory-derived constitutive properties (Extended Data Fig. 2). The transient creep accelerates the immediate relaxation that follows the stress perturbation and the subsequent relaxation is slower, compared to when operating at steady state only. The hardening coefficient G controls the amount of strain that is relaxed by transient creep.
We generalize the constitutive relationship for transient creep for three-dimensional deformation with isotropic rheology. The constitutive relationship becomes
where the subscripts i and j are the tensor indices
is the internal deviatoric stress tensor and q2 = Q Q is the norm of the internal deviatoric stress tensor (we use Einstein’s summation convention). We implement these constitutive relationships in the community code Relax (www.geodynamics.org/cig/software/relax), such that we can simulate three-dimensional models of postseismic deformation that includes afterslip, viscoelastic flow and the transient creep of olivine in a self-consistent manner. Extended Data Figure 2b and c shows the predicted time series at a few SuGAr stations for various values of G .
The dynamics of postseismic deformation is sensitive to the stress preceding the coseismic perturbation in the asthenosphere because of the power-law dependence in the stress–strain rate relationship. Our mechanical model for the power-law flow of olivine includes (1) the background strain rate due to the shear arising from the vertical gradient of horizontal flow, which is associated with the long-term oblique subduction of the oceanic lithosphere below Sumatra, and (2) the pure shear associated with internal deformation by conjugate strike–slip faulting of the Wharton basin along the diffuse boundary between the Indian and Australian plates. We convert the background strain rate to pre-stress, accounting for the water content, temperature and the other physical parameters of olivine rheology. We assume a homogeneous velocity gradient from the surface (plate velocity of 5.6 cm yr−1 oriented 10° N; Fig. 1) to a relative fixed transition zone. The orientation of the conjugate faults in the Wharton basin suggests that the deformation in the lithosphere is horizontal pure shear39 with principal stress orientation between −30° N and −25° N for the compressive component and between 60° N and 65° N for the extensive component (Fig. 3).
The long-term strain rate is the sum of the linear (diffusion creep) and nonlinear (dislocation creep) strain rate contributions
where η and η are the effective viscosities for diffusion and dislocation creep, respectively. While the linear viscosity due to diffusion creep has little effect on the initial postseismic deformation and cannot be directly estimated (Fig. 2b), it does affect the background stress
So a low-viscosity diffusion creep reduces the background stress and the associated strain rate for dislocation creep.
We consider two end-member models for the background strain rate . In a first model, plate motion is accommodated in the mantle across a 100-km-thick region, leading to a background strain rate of 2 × 10−14 s−1. In a second model, the region is 400 km thick, leading to a background strain rate of 5 × 10−15 s−1. Considering these models one at a time, we investigate a range of strain rates for dislocation creep ranging from full background strain rate (corresponding to vanishing diffusion creep in the model) to 1 × 10−17 s−1 (corresponding to dominant diffusion creep at steady state). We find that our geodetic data are best explained with a dislocation strain rate of the order of 10−17 s−1 to 10−16 s−1, corresponding to grain size between 6 mm and 10 mm depending on the background strain rate. This indicates that diffusion creep and dislocation creep have similar strengths in the asthenosphere at steady state. The strength of diffusion creep is lower than the one for dislocation creep for grain sizes lower than 6 mm (Fig. 2b).
The recent great and giant earthquakes along the Sunda megathrust affected the background stress preceding the 2012 earthquake sequence40. To estimate this effect, we evaluate the deviatoric stress change in the asthenosphere at the time of the 2012 event caused by the nearby 2004 M 9.2 Aceh–Andaman and the 2005 M 8.6 Nias earthquakes (Extended Data Fig. 3). The deviatoric stress near the rupture area of the 2012 event at 100 km depth due to these earthquakes is about two orders of magnitude smaller than the coseismic stress change. For simplicity we ignore the stress caused by the Sunda megathrust earthquakes in our simulations.
The strength of the brittle lithosphere can be evaluated from the orientation of the conjugate faults using Byerlee’s law (τ = μ′σ ). The effective coefficient of friction is given by
where θ is the angle between the fault and the principal stress. The direction of the principle stress should bisect the conjugate faults. In the Indian Ocean earthquake rupture, the main conjugate faults are nearly orthogonal (Figs 1 and 3), providing an estimate of the effective coefficient of friction in the range 0.1–0.2. The exact orientation of the rupture fault is still subject to debate6, affecting our estimate of the effective friction coefficient. The inferred small friction coefficient in the Wharton basin reduces the strength of the brittle lithosphere (Fig. 3a).
The thickness of the lithosphere can be independently estimated from the depth of the coseismic rupture of the 2012 M 8.6 Indian Ocean earthquake as substantial coseismic slip occurred down to 60 km depth4. The reason for the great depth extent of the rupture is unclear, but strong dynamic weakening from frictional melting may have allowed the rupture to penetrate below the seismogenic zone. The depth range from 45 km to 60 km where coseismic slip tapers from its maximum value may correspond to the lithosphere–asthenosphere boundary (Fig. 3a).
We build a three-dimensional model of the lithosphere and asthenosphere in which the flow parameters depend on depth, except in the subducted slab where viscosity is infinite, resulting in a three-dimensional rheological model. The background depth-dependent viscosity is thermally activated following a half-space cooling model
with T = 0 °C, where T is the mantle temperature (a free parameter), the plate age t varies spatially following the model of ref. 24, the thermal diffusivity is Κ = kρ −1C −1 with conductivity k = 3.138 W m−1 K−1, the specific heat is C = 1.171 kJ kg−1 K−1, the density is ρ = 3.330 kg m−3, α = 0.4 °C km−1 is the adiabatic temperature gradient, H(z) is the Heaviside function, and z = 100 km is the depth of the thermal boundary layer for a 60-million-old oceanic plate (Fig. 2b). In our models, most of the postseismic deformation occurs around 100 km depth, where the coseismic stress change is high, so the adiabatic temperature gradient has almost no effect on our water content estimates. The elastic slab is constructed using the Slab 1.0 model41 with a uniform thickness of 80 km. Because the accelerated flow is concentrated in areas of high coseismic stress change, immediately below the mainshock, the effect of the slab does not greatly affect the fit to the geodetic data.
We create models of postseismic relaxation where the coseismic stresses from the mainshock and the largest aftershock are potentially relaxed by both afterslip in the lithosphere and viscoelastic flow in the asthenosphere, depending on the rheological parameters. We use the coseismic slip models of ref. 4 to produce the stress perturbation and to define the geometry of faults for afterslip. In our relaxation models, both afterslip and viscoelastic deformation are driven by stress and the coupling between the two mechanisms is taken into account. Afterslip is allowed only on the faults that ruptured coseismically, around the areas of negative coseismic stress change. We assume uniform friction properties outside the rupture area, as only near-field data would motivate finer-tuned models.
We model afterslip with rate-strengthening friction, a simplification of the rate-and-state friction law that is adequate to represent triggered aseismic slip at steady state8, 9, 10. In the rate-strengthening approximation the afterslip velocity V is given by25
where Δτ is the stress change after the mainshock, (a − b) is the steady-state friction parameter and σ is the effective normal stress on the fault. We assume a uniform (a − b)σ = 6 MPa in our simulations and we explore values of V that best explain the geodetic observations in combination with other mechanisms of deformation (Extended Data Fig. 6).
The calculations are performed using the Relax software (www.geodynamics.org), which employs a spectral method to evaluate the quasi-static deformation caused by stress perturbations. The method incorporates an equivalent body-force representation of dislocations that allows us to include detailed finite slip distributions. The equivalent body-force method is advantageous compared to other approaches such as finite-element methods because meshing of the domain around three-dimensional faults is automatic. The numerical approach has been validated using comparisons with analytic solutions for fault slip and comparisons with finite-element solutions in simple geometric settings23, 25, 42, 43. Another advantage of the Relax software and the equivalent body-force method is the ability to simulate several physical mechanisms of deformation simultaneously and to include nonlinear rheologies for plastic flow8, 9, 10.
Afterslip as a single mechanism of postseismic deformation predicts horizontal deformation compatible with the observation. However, the vertical displacement predicted by afterslip in the nearest stations is opposite to observations, suggesting another active process of postseismic deformation, such as viscoelastic relaxation in the asthenosphere (Extended Data Fig. 4).
Viscoelastic flow in the upper mantle can occur by diffusion creep, dislocation creep, or a combination of both7. The physical properties of dislocation creep in olivine are well known from a wealth of laboratory experiments documenting the effect of temperature, pressure7 and water44. The properties of diffusion creep of olivine are less well known because of its great sensitivity to grain size. Several field samples show that both diffusion creep and dislocation creep act in tandem to deform mylonite shear zones below the brittle-to-ductile transition45. Experimental work suggests that dislocation creep of olivine is the dominant mechanism of deformation in the upper-mantle conditions2, 46. But competition between grain growth during diffusion creep and grain size reduction by dislocation creep can promote comparable strain rate at these depths at equilibrium. In this case, dislocation creep should dominate postseismic deformation following a large stress perturbation because of the power-law stress–strain rate dependence. As the stress is reduced by postseismic relaxation the contributions of both mechanisms should become similar again.
Following these considerations we first tested the potential of dislocation creep to explain the postseismic transient. Initially, we ignored the transient creep of olivine. We find that nonlinear viscoelastic deformation produces vertical displacements compatible with observations and horizontal displacements in the correct azimuth. However, the amplitude of horizontal displacements is lower than observed. We conclude that steady-state dislocation creep cannot explain the entire set of observations. Instead, we find that only a combination of afterslip on the Indian Ocean coseismic faults and viscoelastic flow can satisfactorily explain the GPS observations.
We use a Bayesian approach to estimate the water content in olivine, assimilating prior information from geochemical estimates and additional constraints from geodetic observations. The posterior probability density is47
where ν is a constant, d is the GPS data vector, g(m) is the forward model based on parameters m, C is the data covariance matrix and ρ is the prior information. Geochemistry provides limits on the water content in the mid-ocean ridge basalt source. We include these constraints on water content using a log-normal distribution (Fig. 2a). For the afterslip parameter V , we assume the uniform prior 1/V , which corresponds to the limit of a log-normal distribution for infinite variance (Extended Data Fig. 6a). We assume that the geodetic data are independently and normally distributed. In principle, the joint probability density of all physical parameters can be estimated, but we limit our exploration to two physical parameters at a time to reduce the computational burden. In a first step, we jointly estimate the water content in olivine and the rate-strengthening parameters for afterslip, and assume all other parameters to be fixed. In a second step, we change the in situ parameters manually to explore a range of conditions.
We sample the posterior probability density using the Neighbourhood algorithm48, a derivative-free Markov chain Monte Carlo method. The inversion tool, dubbed Relax-Miracle, uses Relax for the forward models. We benchmark the approach using a synthetic data set using the GPS network configuration of the Sumatra GPS Array. We create a forward model with a known water content and afterslip parameter and we use this data set as an input data to our inversion scheme. The posterior probability density function not only recovers the target parameters but also informs us of the inherent tradeoffs between the two relaxation mechanisms (Extended Data Fig. 5).
We explore uniform water content C in olivine in the range 0.0003 to 0.04 wt% (50 to 6,000 H atoms per million Si atoms), mantle temperatures T from 1,350 °C to 1,400 °C, and reference velocities for afterslip in the range V = 0–2.75 μm s−1. We explore different values of the transient creep parameters A from A = 0 to A = 3A, and G from G/2 to 3G. We assume all other physical and in situ parameters the same for transient creep and steady-state creep. We assume that the laboratory-derived values for the constitutive parameters, which were carried out at much higher strain rate than at typical geological conditions, scale to natural conditions. Geodesy cannot independently constrain the water sensitivity of transient creep and the parameter A , so the product A (C )r for transient creep should be considered as a lump parameter. We also investigate a range of strain rates for dislocation creep from 10−17 s−1 to 2 × 10−14 s−1. Our best-fitting model in Fig. 2 has V = 1.75 × 10−6 m s−1, A = A/2, and G = G. Incorporating transient creep affects the estimate of water content in the mantle, as previously inferred in other studies49, 50, by lowering the water content required to fit the data (Fig. 2a). The inferred afterslip parameter V is similar to what is found in other tectonic settings25, 51.
We obtain the probability density function of the water content in the asthenosphere by either marginalizing out the afterslip parameter
where m is the water content in olivine and m is the afterslip parameter; or taking the conditional probability density function around the most likely value of the bivariate distribution
where m is the most likely value of the afterslip parameter. The posterior probability density function for water content and afterslip in the Wharton basin is shown in Fig. 2a and Extended Data Fig. 6.
The best-fit forward model is shown in Fig. 2 and Extended Data Fig. 7. The small misfit in the GPS time series may be due to our simplifying modelling assumptions, which ignore reactivation of the megathrust or the Sumatran fault and internal deformation of the accretionary prism.
The numerical software used in this study is hosted at https://bitbucket.org and is available from the corresponding author on request.