Earthquakes and tectonics I

On 14^th November 2016, the north eastern coast of New Zealand’s South Island was struck by a large Mw 7.8 earthquake. The event occurred on the southern edge of the Marlborough Fault Zone and was the largest event to hit the area in more than 100 years. Shaking was widely felt throughout the whole of New Zealand with widespread damage across the northern South Island. A multitude of datasets including InSAR, GPS, field observations and seismology have revealed this to be one of the most complex earthquakes ever recorded. The rupture extends for more than 170 km across both mapped and unmapped faults, on and offshore. Many of the faults have slipped at the surface generating new fault scarps up to 10 m in height. Analysis of high rate GPS data indicates two discrete events separated by ~30 seconds starting in the south and moving north. Slip modelling suggests up to ~15 m of dextral and reverse slip along three main branches of the Marlborough fault zone with slip extending down to ~20-30 km depth. In addition, a fault bounded block located along the Kaikoura coast appears to have been thrust out of the ground by up to 10 m during the earthquake.

The 24th September Mw 7.7 Balochistan earthquake ruptured over 200 km of the Hoshab Fault in SW Pakistan. The fault is curved and dips moderately (50-70°) to the north, but gave predominantly strike-slip displacements in the earthquake, with left-lateral displacements up to ~14 m.  The slip vector for the 2013 earthquake therefore rotates by almost 60° along the length of the rupture, prompting a variety of models. The postseismic deformation following an earthquake of this size presents an opportunity to probe the structure, rheology and stress state of the crust in the transition region between the E-W oriented faulting related to the Makran subduction zone and the N-S faulting of the Chaman Fault zone.

We construct a Sentinel-1 InSAR time series of deformation covering the full 200 km length of the fault from two look directions, spanning the 2 year period from Nov 2014. We observe up to ~10 cm of deformation in satellite line-of-sight, but in contrast to the coseismic slip which peaked on the shallowest part of the fault, the surface displacement suggests longer wavelength and deeper deformation. A localised region of deformation at a fault bend near the epicentre is attributed to stress concentration arising from the geometric complexity. We model the InSAR time series to invert for the spatio-temporal evolution of slip at depth.

The Sentinel-1 time series also images the creeping segment at the northern end of the Hoshab fault, previously noted by Fattahi and Amelung 2016. The end of mapped surface ruptures coincides exactly with beginning of the creeping section, confirming previous suggestions that the rupture terminated at the creeping segment. Following the 2013 earthquake, this section of the fault shows enhanced creep, at up to 1 cm/yr line of sight displacement, an order or magnitude faster than the pre-seismic creep rate. The western tip of the fault shows a similar period of creep at the surface, extending from the end of the coseismic rupture. These observations are consistent with the earthquake perturbation to the regional stress field suggested by coulomb modelling.

The 2013 Mw 7.7 Balochistan earthquake occurred in the complex tectonic setting of the triple junction between Arabia, Eurasia and India tectonic plates. The earthquake ruptured a 200-km-long curved section of the Hoshab fault, after nucleating on the Chaman fault. Coseismic motion was dominated by left-lateral slip with a minor reverse component. Since the Hoshab Fault is mainly a thrust fault on which the rupture would have been expected to be mainly dip-slip, the strike-slip motion induced by the 2013 earthquake raises important questions on the mechanics of earthquake faulting.

We process TerraSAR-X ScanSAR data and RADARSAT-2 data using both interferometry and amplitude correlation. We also apply optical image correlation, using two pairs of Landsat-8 images together with a set of five pairs of SPOT-5 images, to cover the rupture on its entire length. By combining our radar and optical datasets, we first derive the full 3D coseismic displacement field at the surface. Retrieval of the fault-parallel and fault-normal components of slip along the fault enables us to show that surface slip exceeds 6-7 m over a distance of more than 100 km, with a maximum of 12 m in the central part of the rupture. The vertical component increases away from the nucleation area, reaching up to 3 m towards the southern extremity of the fault. Our analysis shows evidence for a North dipping fault, as already suggested in previous studies. Comparison of the relative vertical displacement across the fault against fault-normal motion yields a fault dip ranging between 45° and 70°.

We carry out an elastic inversion of the geodetic dataset in order to determine the slip distribution at depth. To determine the first-order features of the slip distribution, we use a simple geometry for the fault. The northern part is modeled as a 70°N dipping fault down to 18 km, consistent with linkage with the strike-slip Chaman Fault further to the North. The fault dip is gradually flattened towards the South to comply with the thrust fault morphology in the Makran accretionary wedge and with the analysis of on-fault relative displacements. At depth, the fault is modeled with a listric geometry that flattens at 10 km, consistent with a previously inferred décollement level. We use Okada's equations to invert for the strike-slip and dip-slip components of the earthquake. The proposed model explains up to 90% of our dataset. We show that the strike-slip component is segmented in two parts, with a slip of 10 m on a 75 km long section to the North, and a second section with a slip of 5-6 m on a 75 km long section to the South. The dip-slip component is mostly restricted to the southern section, with an average 1.5 m reverse slip. Transition between the two segments occurs around a major geometric complexity visible along the fault trace. This suggests that the fault geometry exerts a control on the coseismic slip distribution.

Synthetics are however almost systematically lower than the original data at short distance from the fault (less than 2 km). This suggests that the inversion cannot reconstruct inelastic deformation or apparent slip overshoot due to shallow flattening of the fault.

Furthermore, the predicted asymmetry of off-fault displacements in the mid-field (between 4 and 10 km) is systematically underestimated by the model. This suggests that the fault dip at depth could be steeper than previously proposed from the analysis of on-fault relative displacements. This result would be consistent with a steepening of the fault dip with depth, again suggesting that complexities within the geometry of the Hoshab Fault strongly influenced the coseismic slip distribution of this large earthquake. 

The extent of seismic versus aseismic slip along major faults is controlled by the rheological properties of the fault interface and the state of stress. Mapping subsurface fault slip during the different phases of the seismic cycle provides a probe of the mechanical properties these faults. Here, we focus on the megathrust in northern Chile, where the Pacific plate subducts beneath South America at an average rate of approximately 6.7 cm/yr. Using GPS-derived displacement rates, first order estimates of the distribution of fault locking in the interseismic period suggests little to no overlap in regions slipping seismically versus those that are dominantly aseismic. While the spatial distribution of slip associated with the 1868 and 1877 Mw 8+ earthquakes are relatively unknown, recent earthquakes, including the 2007 Mw 7.7 Tocopilla and the 2014 Mw8.1 Iquique earthquake ruptured portions of the megathrust that were inferred to be locked beforehand, confirming the first-order frictional description of active faults commonly assumed. However, most published distributions of slip, be they during seismic or aseismic phases, rely on unphysical regularization of the inverse problem (smoothing, damping), thereby cluttering attempts to quantify the degree of overlap between seismic and aseismic slip. Considering all the implications of aseismic slip on our understanding of the nucleation, propagation and arrest of seismic ruptures, it is of utmost importance to quantify the maximum spatial and temporal overlap of seismic and aseismic slip with corresponding uncertainties. Here, we take advantage of 20 years of InSAR observations and more than a decade of GPS measurements to derive probabilistic maps of inter-seismic coupling, as well as co-seismic and post-seismic slip along the northern Chile subduction megathrust.

We use InSAR observations from the ERS, Envisat and ALOS satellites to extract maps of ground displacements rates between 1993 and 2010. Because of the significant computational burden, most time series analysis methods developed in the past rely on a pixel-by-pixel approach. Such methods require a common reference for all interferograms by removing an empirically determined, long spatial wavelength field in all images. Furthermore, pixel-by-pixel methods ignore spatial covariances in interferograms. Residual tropospheric perturbations may be considered isotropic and their spatial distribution can be statistically described by an exponential decay as a function of distance between pixels. Accounting for such covariance pattern allows one to explicitly account for turbulent tropospheric perturbations in InSAR time series analysis. We have developed a method that simultaneously considers time series of all the pixels, accounting for the full spatial covariance between pixels of each interferogram. We consider the interferometric phase as the sum of the phase difference between two SAR acquisitions and of a long spatial wavelength field, here a linear function of range and azimuth describing what is considered as residual orbital contribution. We also include a parameterized evolution of the deformation through time, in order to tie disconnected observations from the ERS and Envisat satellites. We regularize our solution in space using an exponential covariance model. We therefore solve for approximately one million model parameters (i.e. the phase evolution for each pixel plus the orbital parameters) with approximately ten million data points (i.e. the pixels of all interferograms). Using full covariances allows to reduce the number of free parameters through a physical description of the residual noise in our data. Given the size of the estimation problem, we have adopted a conjugate gradient solver in a MPI-based, sparse formalism using the Python librairies petsc4py and mpi4py based on PETSc. We avoid direct matrix multiplications imposed by the use of full covariances by computing convolutions in the Fourier domain. Without inputs from additional GNSS data, we recover a continuous map of the rates of deformation from the Mejillones peninsula to the Arica bend in Northern Chile.

These deformation rates are combined with seismological data and available GPS observations collected by permanent networks installed by the École Normale Supérieure and Institut de Physique du Globe (Paris, France), the GeoForschung Zentrum (Potsdam, Germany) and the California Institute of Technology (Pasadena, CA, United States of America) to derive a complete probabilistic description of slip during all phases of the earthquake cycle in the past 20 years. We use AlTar, a massively parallel Monte Carlo Markov Chain algorithm exploiting the acceleration capabilities of Graphic Processing Units, to derive the probability density functions (PDF) of slip given the seismic and geodetic data available. Our solution accounts for full covariances between observations, including uncertainties on the GPS-derived rates and covariances between InSAR pixels due to spatially coherent noises such as tropospheric residuals. We compute the Green's functions, the surface displacement for a unit slip on each point of the megathrust, in an elastic layered space and account for the uncertainties on the elastic structure of the earth using a perturbation approach. We derive a probabilistic map of slip for the inter-seismic period preceding the 2014, Iquique earthquake, a probabilistic map of co-seismic slip for that earthquake and the corresponding map of post-seismic slip.

The PDFs we derive allow quantification of slip over three phases of the earthquake cycle, the mean model of slip (i.e. corresponding to what is usually derived using classic least-squares approaches and gaussian statistics), the most probable model and, more importantly, the whole range of models allowed by the seismic and geodetic data. Beyond the confidence levels in our solutions, we are able to answer fundamental questions in a probabilistic, principled way. In northern Chile, we find high probabilities for a complete release of the elastic strain accumulated since the 1877 earthquake by the 2014, Iquique earthquake and for the presence of a large, independent, locked asperity left untapped by recent events, north of the Mejillones peninsula. We evaluate the probability of overlap between the co-, inter- and post-seismic slip and consider recent developments suggesting the occurrence of slow, aseismic slip events along this portion of the subduction zone. Our model confirms previous estimates of the degree of locking of the subduction interface and of seismic slip, providing additional necessary information, the uncertainties, given by a systematic exploration of all possible models of slip along the northern Chile subduction zone.

We mapped complex fault ruptures for a number of large earthquakes in 2016, including the February 2016 MeiNong earthquake in Taiwan, the April 2016 Kumamoto earthquake sequence in Japan and the central Italy sequence, using analysis of SAR data from the Copernicus Sentinel-1A (S1A) and Sentinel-1B (S1B) satellites operated by the European Space Agency and the Advanced Land Observation Satellite-2 (ALOS-2) satellite operated by the Japanese Aerospace Exploration Agency (JAXA). We find that triggered slip on faults near main ruptures occurred during or soon after many of these events. The MeiNong main rupture at lower crustal depth triggered slip on another fault at upper crustal depth and shallow slip on several faults in the upper few km of southern Taiwan. The Kumamoto earthquake sequence ruptured two major fault systems over two days and triggered shallow slip on a large number of shallow faults in Japan. We combine less precise analysis of large scale displacements from the SAR images of the two satellites by pixel offset tracking or sub-pixel correlation and by burst overlap double-difference interferograms on ALOS-2 ScanSAR pairs, including the along-track component of surface motion, with the more precise SAR interferometry (InSAR) measurements in the radar line-of-sight direction to estimate all three components of the surface displacement for the events. Data was processed with customized workflows based on modules in the InSAR Scientific Computing Environment (ISCE).

Joint inversion of S1A and ALOS-2 InSAR, GPS, and strong motion seismograms for the Mw6.4 MeiNong earthquake shows that the main thrust rupture with N61°W strike and 15° dip at 15-20 km depth explains nearly all of the seismic waveforms but leaves a substantial uplift residual in the InSAR and GPS offsets estimated 4 hours after the earthquake. We model this residual with slip on a N8°E-trending thrust fault dipping 30° at depths between 5-10 km. This fault strike is parallel to surface faults and we interpret it as fault slip within a mid-crustal duplex that was triggered by the main rupture within 4 hours of the mainshock. In addition, InSAR shows sharp discontinuities at many locations that are likely due to shallow triggered slip, but the timing of these is uncertain.

The Kumamoto earthquake sequence in Japan started with Mw 6.2 and 6.0 earthquakes on 14 April (UTC) followed on 15 April by the Mw 7.0 mainshock. JAXA acquired one ALOS-2 scene between the foreshocks and mainshock that enables some separation of the surface deformation. InSAR shows M6 foreshocks were deeper, while M7 mainshock ruptured surface.

Earthquakes and tectonics I