Illuminating subduction zone rheological properties in the wake of a giant earthquake

We invert postseismic geodetic data to reveal subduction zone mechanical properties following a megathrust earthquake in Chile.


This PDF file includes:
Supplementary Materials and Methods Fig. S1. Continuous GPS network and interpolated interseismic velocity. Fig. S2. CGPS time series, empirical fits, and afterslip and viscoelastic flow modeling-based time series (also see subsequent pages). Fig. S3. Inversion parameters and sensitivity tests for select volumes. Fig. S4. Five checkerboard tests for slip on the megathrust and strain in the polyhedral volumes, which demonstrate our ability to resolve any up-and down-dip slip and viscous strain in the ductile region using our CGPS site distribution (white triangles).  Table S1. Dislocation creep rheological parameter estimates for the Maule region. Table S2. Temperature estimates for the polyhedra. References (48-67)

GPS data processing and preparation
For the Maule postseismic analysis we include continuous GPS (CGPS) data from a large regional geodetic network ( fig. S1). For each daily solution we process all available data using GAMIT (48), which estimates the 3-D relative position of ground stations and satellite orbits and GLOBK (49), which merges the daily regional with global solutions for International Global Navigation Satellite Systems Service (IGS) sites. Since we are interested primarily in regional deformation, we refrain from the global stacking common to large-scale scale plate motion studies and instead stack only stations on the South America and Nazca Plates and a few sites on adjacent plates and minimize motions of stations located on the stable core of the South American craton. We express our daily time series in a craton-fixed reference frame with horizontal and vertical velocities of ~1mm/yr. Specifically, we use 44 CGPS stations primarily located in Brazil to create the horizontal reference frame. These sites have a total RMS horizontal velocity of 0.75 mm/yr. The total RMS velocity of the 35 sites used to constrain the vertical reference frame is 1.13 mm/yr. Our GPS field, processing, and time series analysis methods in South America are described in detail elsewhere (50).
Prior to modeling, we isolate the postseismic portion of the CGPS time series by removing an estimate of the interseismic velocity. For sites where sufficient data exist prior to the Maule earthquake (~2 years) we do this by simultaneously estimating the best-fitting interseismic velocity, coseismic step, postseismic decay function, and seasonal component in a least squares sense (50). We subtract the interseismic and coseismic contributions to the time series but not the estimated seasonal term as our tests suggest that subtracting this component often results in an increase in the time series scatter.
For CGPS sites with little or no pre-Maule data, we require an estimate of the associated interseismic velocity, which we derive by combining the best-fit North, East, and Up velocities from sites with pre-quake data (described above) with published interseismic velocity estimates for the Maule region ( fig. S2) (24,(51)(52)(53)(54)(55). We use the most recent velocities for sites with numerous solutions and only include those expressed in a South America-fixed reference frame, which show relative uniformity and no noticeable reference frame bias. We calculate a spatially averaged dataset using the velocity uncertainties as weighting factors, create interpolated North, East, and Up velocity fields ( fig. S2), extract velocities for locations with little to no pre-quake data, and remove these estimates prior to computing the best-fit postseismic decay functions.
For daily position uncertainties we calculate the mean and standard deviation of the time series in a 7-day moving window, multiply the standard deviation by 1.5, and assign this value as the uncertainty for daily positions within the corresponding time window. We fill data gaps using the aforementioned best-fit postseismic decay functions and assign uncertainties to these daily positions by taking the mean of the weekly standard deviation estimates and multiplying this value by 10 so that the empirical fits across data gaps are down-weighted during the afterslip and viscoelastic strain inversion.
We use SLAB1.0 (56) to define the geometry of the subduction zone interface, which we discretize into a network of triangular boundary elements that extend along strike for ~400 km and from the trench to a depth of ~75 km (main text Fig. 3). We simultaneously invert the postseismic surface displacement time series at each CGPS site for localized afterslip on the megathrust and distributed deformation in the ductile regions using a modified version of the Kalman-based extended network inversion filter (7,57,58), which provides conditional probability estimates for the model parameters based on uncertainty in the data. The model parameters in this case correspond to the afterslip magnitude and direction for each of the triangular boundary elements (59) and six independent components of strain for each of the polyhedral volume elements (28) and their associated rates. The volumes vary in size, but their bases are typically 50-100 km wide/long and their heights range from 10-100 km. A linear summation of the six strain components determines the total deformation of each volume element and the 3D Green's functions map the strain at depth to the surface displacements. The amplitudes of the strain and fault slip Green's functions determine our ability to successfully explain the associated surface displacements and reveal tradeoffs between the two different deformation mechanisms. Due to the large number of model parameters we make simplifying assumptions to make the inversion more tractable and to force smooth parameter distributions.
These include a no back-slip, or positivity penalization, during the afterslip inversion and we assume that the viscoelastic flow is deviatoric, penalizing strain directions perpendicular to the induced coseismic stresses and promoting strain directions that align with the coseismic stress directions. We also penalize afterslip so that it preferentially occurs in the region surrounding the coseismic rupture patches and we apply spatial smoothing using stress kernels to both the afterslip and strain components to avoid large variations in slip or strain between neighboring Modeling approach elements that over-fit the postseismic time series. The temporal and spatial smoothing hyperparameters are estimated simultaneously with the model parameters (7,57,58). Following the previous work (7), these smoothing hyper-parameters are determined iteratively until the best fit to the GPS measurements is achieved.
It is impossible to directly measure the background stress and strain rate of the upper mantle.
Therefore, we again follow previous work (7) when estimating the total stress for evaluating the time evolution of effective viscosity for each deformable volume and consider three stress components; the initial stress from the coseismic rupture, the stress evolution for each strain volume, and the background stress corresponding to the steady-state viscosity multiplied by the background strain rate. We perform a grid search for the background strain rate and steady-state viscosity that best explain the time evolution of the effective viscosity curves compared to the Burgers rheology model (Figs. 5 and S6;7). The minimum misfit provides the most likely combination of both background viscosity and strain rate that matches the trend of stress relaxation within the polyhedral volumes and the corresponding effective viscosity time series.

Fault friction
To explore the frictional properties of the megathrust we assume that the temporal evolution of afterslip obeys the rate-and-state friction law under steady-state, with ss = * n + ( − ) n * where a and b are the rate-and-state frictional parameters (a-b describes the fault behavior), the driving shear stress, the effective normal stress, * the reference coefficient of friction, and * and the reference slip velocity and slip velocity, respectively. We estimate (a-b) (equivalent to the ratio of Coulomb shear stress change to slip rate) for each fault patch where CFS is the time-varying Coulomb stress on the creeping fault element (4, 10, 31, 60). If a-b > 0 the material is velocity strengthening and primarily stable sliding (i.e., creep) is possible.
In contrast, if a-b < 0, the material is velocity weakening and potentially seismic rupture (i.e. frictional instabilities may nucleate).
Steady-state rheology, temperature, and activation energy estimates The rheological properties (i.e. constitutive stress-strain rate relationship) of olivine aggregates in the mantle determined from deformation experiments are described using a flow law in the where ̇ is the strain rate, A is a pre-exponential factor, is deviatoric stress, is the stress exponent, is the activation energy, is the universal gas constant, is the temperature, is the activation volume, is the confining pressure, is grain size, is the grain size exponent, and and are the water concentration and associated exponent (16,41,61). Diffusion creep is associated with = 1 and = 2 − 3, and dislocation creep is associated with = 3 − 5 and = 0. To calculate the geothermal gradient of the continental lithosphere-asthenosphere system required for the depth-dependent viscosity we solve the steady-state, 1-D heat conduction equation with radiogenic heat production using Equation (3)  Based on these best-fit terms, we can further investigate the temperature and water content (i.e. activation energy) of the Maule lower crust and upper mantle with two end-member approaches.
The first involves using the best-fit rock properties from the 1-D profiles ( Fig. 5 and table S1) to estimate the 3-D temperature structure. We assume the rocks within the oceanic mantle, continental mantle and lower crust have the same properties (e.g. the mantle wedge has the same best-fit rock properties as the continental mantle). With this assumption and the inverted strain rate and stress over the steady-state stage (e.g., Figs. 5 and S5), we solve for the temperature of For our third forward calculation we recall that the inversion results suggest the mantle wedge is strong (i.e. relatively high viscosity) and thus is less deformed than the surrounding rocks during the postseismic period. Therefore, we penalize strain in the mantle wedge volumes and find that the model does a much better job of fitting both the horizontal and vertical displacements compared to the previous forward models. The penalized model predicts coastal subsidence and uplift across the forearc but subsidence across the volcanic arc. All of these features are consistent with both the data and the best model from the inversion pointing towards a powerlaw rheology.
To summarize, our UniCyclE 3D stress-driven forward models confirm that the inversion results provide a good first order estimate of the subduction zone rheological structure. A nonlinear power-law rheology is preferred, and the strong mantle wedge corner is a robust feature from both the inversion and forward modeling exercises. We reserve additional forward modeling exercises for a follow-up study.         *To obtain the optimized dislocation creep parameters listed above we start with values reported in Karato and Jung (2003). Specifically n =3.0±0.1, Q= 510±30 kJ/mol and 470±40 kJ/mol for dry and wet olivine, respectively, V =11x10 -6 m 3 /mol; C OH = 1000 H/10 6 Si, r =1.2±0.05, A =90 MPa -n s -1 , and strain rate=10 -14 s -1 .