Imaging the distribution of transient viscosity after the 2016 Mw 7.1 Kumamoto earthquake

See allHide authors and affiliations

Science  14 Apr 2017:
Vol. 356, Issue 6334, pp. 163-167
DOI: 10.1126/science.aal3422

Crustal rock strength from outer space

The response of crustal rock to stresses is challenging to estimate yet vital for determining risks from events such as earthquakes. Moore et al. take advantage of the recent Mw 7.1 Kumamoto earthquake in Japan to determine the rheology of crustal rocks in the region. The observed inversion of the crustal strain rates demonstrates that certain areas have stiff rock and others (e.g., under the Aso volcanic complex) have much weaker rock. The results match up with expectations, which means that the method can successfully measure rock properties over a wide range of strength and large spatial and temporal scales.

Science, this issue p. 163


The deformation of mantle and crustal rocks in response to stress plays a crucial role in the distribution of seismic and volcanic hazards, controlling tectonic processes ranging from continental drift to earthquake triggering. However, the spatial variation of these dynamic properties is poorly understood as they are difficult to measure. We exploited the large stress perturbation incurred by the 2016 earthquake sequence in Kumamoto, Japan, to directly image localized and distributed deformation. The earthquakes illuminated distinct regions of low effective viscosity in the lower crust, notably beneath the Mount Aso and Mount Kuju volcanoes, surrounded by larger-scale variations of viscosity across the back-arc. This study demonstrates a new potential for geodesy to directly probe rock rheology in situ across many spatial and temporal scales.

Crustal dynamics involves the nonlinear interactions of faulting, ductile flow, and fluid migration, entwined with a complex thermal and metamorphic history. In addition to elastic deformation in the crust, there are a number of anelastic deformation mechanisms, including slip on faults in earthquakes and off-fault distributed ductile flow. The spatial distribution of the viscous properties in particular is poorly known, yet these properties play a key role in plate tectonics (1) and earthquake triggering (24). Geodetic observation and modeling illuminates a range of these anelastic mechanisms activated in the aftermath of earthquakes. In particular, faults continue creeping aseismically for several months after the mainshock (afterslip), with accelerated ductile flow at depth (5) and in the surrounding rocks. Many studies exploit space geodetic data to infer the viscous properties of the lower crust and mantle using sophisticated models of stress relaxation (69), revealing the nonlinear rheology of rocks and transient creep, whereas other studies employ proxies such as magnetotellurics (10). Despite these efforts, knowledge of the spatial distribution of anelastic properties and mechanical strength variations remains limited. This is due in part to the complexity of distinguishing between deformation mechanisms, which often produce similar surface deformation (11). Additionally, kinematic models of crustal deformation have thus far been limited to slip on faults, hindering inference about the kinematics and rheology of distributed deformation. We address this by treating the kinematics of fault slip and distributed deformation concurrently, using novel Green’s functions for distributed deformation (12).

We exploited the spatial and temporal resolution of interferometric synthetic aperture radar (InSAR) observations from the European Sentinel-1A and Japanese Advanced Land Observing Satellite 2 (ALOS-2) satellites combined with dense Global Positioning System (GPS) measurements from the Japanese GNSS Earth Observation Network System (GEONET) to directly image the localized and distributed deformation that followed the 2016 Kumamoto, Japan, earthquake sequence (Fig. 1). We assimilated InSAR line-of-sight (LOS) and GPS displacement time series spanning before, during, and after the Kumamoto earthquakes to describe the details of coseismic slip, afterslip, and lower-crustal ductile flow. We found large variations of viscosities across Kyushu, primarily due to thermal fluctuations in the lower crust related to the volcanic arc.

Fig. 1 Postseismic deformation and coseismic slip from the 2016 Kumamoto earthquake sequence.

Postseismic GPS velocities from the GEONET array (black triangles indicate stations, with white arrows for station velocities). The stations north of the rupture predominantly move rapidly northward, while those south of the rupture move southward with comparatively smaller velocities. Known faults are mapped in light brown, the subducting Ryukyu plate is mapped in dashed contours, and volcanoes are marked with red triangles. We have included the National Research Institute for Earth Science and Disaster Resilience (NIED) moment tensor solutions for the foreshock, mainshock, and large aftershocks with known tensor solutions. Pink dashed lines indicate the graben boundary. The inset illustrates the result of our geodetic inversion for the coseismic slip distribution on the Hinagu and Futagawa faults, including contributions from both the foreshock and mainshock.

The Kumamoto earthquake sequence began on 14 April 2016 at 21:26 JST (Japan Standard Time) when a moment magnitude (Mw) 6.1 foreshock struck the southern Japanese island of Kyushu (13). It was closely followed on 16 April 2016 at 01:25 JST by the Mw 7.0 mainshock, which resulted in the loss of 72 lives (14) and widespread structural damage due to the intense ground shaking. These earthquakes and their aftershocks occurred on the Futagawa and Hinagu fault zones (Fig. 1), with a measured right-lateral strike-slip surface offset of 2 m (15) and a coseismic rupture front that terminated beneath Mount Aso (16). The Futagawa and Hinagu fault zones make up the westward extension of the Median Tectonic Line and take up 25% of the margin parallel component of deformation at the Ryukyu trench (17). In addition, there is north-south extension within the Beppu-Shimabara graben (18) and arc volcanism associated with the Ryukyu trench coincident with these structures (Fig. 1).

We inferred the kinematics of the earthquake sequence by assimilation of GPS coseismic offsets and InSAR data, using a combination of ascending and descending orbits. We inverted ALOS-2 InSAR images and GEONET GPS data to obtain the total coseismic rupture of both the foreshock and the mainshock (19). The coseismic slip inversion shows predominantly right-lateral strike-slip motion, a peak slip of 5 m, and surface rupture of ~2 m, with each earthquake occurring on separate fault segments (Fig. 1). We also find a normal-slip component on the Futagawa fault. These results are in close agreement both with those found by means of seismic waveform inversion (20) and the measured surface rupture (15). The two earthquakes induced a cumulative coseismic stress change on the order of 1 MPa in the lower crust, large enough to trigger accelerated viscoelastic flow.

We extracted the 91-day postseismic deformation from the GPS time series for 321 stations (19), revealing a coherent velocity field at Earth’s surface (Fig. 1), presumably caused by a combination of transient distributed deformation and afterslip. The greatest-magnitude velocities are generally found in the northern half of Kyushu. In the near field, there are a number of large-amplitude, short-wavelength features, indicative of local aftershocks and shallow aseismic afterslip. A few GPS stations exhibit anomalous velocities, most notably the one situated next to Tsurumi volcano in the Beppu-Shimabara graben, moving in the opposite sense to its neighbors, with the largest magnitude of any station. A Mw 5.7 aftershock near Tsurumi volcano triggered 30 s after the mainshock was responsible for the observed anomalous velocity, as we clearly saw a stepwise eastward offset in the GPS time series. The GPS velocity field may also be affected by the aseismic motion of small cracks near the surface (21). Additionally, we unwrapped the Sentinel-1A InSAR frames and used the GPS stations to correct orbital aberrations, providing absolute LOS velocities.

The kinematics of the postseismic deformation after the Kumamoto earthquake sequence can be largely reduced to localized slip on faults and distributed strain in finite volumes (22, 23) (Fig. 2 and fig. S4), corresponding to afterslip on the Futagawa and Hinagu fault zones and ductile flow in the lower crust, respectively. We divided up the Futagawa and Hinagu faults into 2-km by 2-km square patches (24), and the lower crust into 15-km by 15-km by 7.5-km cuboid volumes (12). We aligned the strike of the strain volumes with the average fault strike and located them in the bottom 15 km of the plate beneath Kyushu using Litho 1.0 for the Mohorovičić discontinuity (Moho) depth (25). This resulted in an allowed anelastic deformable region with depth range between 20 and 40 km. In addition, we included six cuboid strain volumes (30 km by 125 km by 30 km) in the mantle wedge, located between the Moho and the downgoing Philippine Sea plate at a depth of 50 to 80 km.

Fig. 2 Distributed and localized deformation imaged as lower-crustal strain rates under Kyushu, with afterslip rates on the Hinagu and Futagawa faults, immediately after the 2016 Kumamoto earthquake sequence.

(A) Schematic geometry of the lower-crustal strain volumes and upper-crustal faults. Right-hand side illustrates the six degrees of freedom for a deformable strain zone (where ε represents the individual components of the strain tensor) and two degrees of freedom for a fixed fault patch. (B) Strain rates in the lower-crustal strain volumes (depth 20 to 40 km) are shown by using the second invariant of the strain tensor in map view with afterslip rates on the fault and contours of coseismic slip (lower-right inset). Distribution of strain rates track the orange contours of coseismic stress change (assumed rigidity of 30 GPa) closely, with two notable exceptions: higher strain rates are observed beneath Mount Aso and Mount Kuju. GPS stations are indicated with black triangles, known faults are mapped with red lines, and volcanoes are indicated with red triangles.

We looked for the linear combination of strain rate in finite volumes and slip velocity on the Futagawa and Hinagu faults that best explained the instantaneous postseismic velocities in the GPS and InSAR data. Ductile flow in each of the deformable volumes comprises six independent directions of strain (fig. S4) and faulting requires two directions of slip, giving rise to a large number of free parameters. To avoid overfitting the data, we regularized the problem by penalizing isotropic strain rates and directions that are orthogonal to the induced coseismic stress change, and by minimizing the postseismic stress change rate, a form of smoothing. Slip is also penalized in the coseismic region and in directions orthogonal to the forces induced by the coseismic rupture (19).

Even if the excess of free parameters from the deformable strain volumes and fault patches are reasonably well balanced by a wealth of observational data and additional constraints, inferring the spatial distribution of localized and distributed deformation at crustal depths represents a substantial challenge. We constructed an outlier-insensitive hierarchical Bayesian model (19) by incorporating all constraints into the priors. We then used a Markov chain Monte Carlo algorithm to sample from the posterior distribution of the slip- and strain-rate parameters, and the constraint weights given for the GPS and InSAR data, resulting in estimates of these parameters as well as their uncertainties. This allowed us to carefully balance overfitting the data with possibly significant outliers and drawing inferences as to the most likely cause of deformation beneath Kyushu.

We illustrated the robustness of our approach by carrying out a series of synthetic tests (19), including checkerboards (figs. S6 and S7) and the recovery of forward models where we assumed uniform and spatially varying viscosity in the lower crust (figs. S8 to S10). To attempt to replicate the perils of using real data, we randomly introduced noise and large outliers into the displacement fields. The recovered models indicate that we can clearly differentiate between viscous flow and afterslip, though the viscoelastic flow is smoothed from the input. Perhaps unsurprisingly, the best resolution is obtained where we have both GPS and InSAR coverage, and the worst is located in the submarine regions. The smoothing of the recovered model is imposed to prevent overfitting the noise in the data or velocities due to sources not incorporated in the model. The synthetic tests also allowed us to quantify the error associated with the recovered viscosities of 30%, or half an order of magnitude.

Our approach allowed us to intrinsically disentangle the complex interactions between afterslip and ductile flow and to describe their relative contributions to the surface velocities. We found shallow afterslip on the Hinagu fault, with a broader region of afterslip on the eastern end of the Futagawa fault (Fig. 2), with velocities in the range predicted by typical frictional parameters under this stress perturbation (26). The strain rates in the lower crust broadly follow the independently obtained contours of coseismic stress change (Fig. 2), as expected for ductile flow, with distinct regions of high strain rate, most notably beneath the volcanic edifices of Mount Aso and Mount Kuju. Because of smoothing, the technique does not currently resolve differences between the top and bottom of the lower crust. The preferred model reproduces the observed InSAR data very closely, with 91% variance reduction (Fig. 3 and fig. S11), and matches the GPS data reasonably closely, with 84% variance reduction. The far-field velocities are dominated by viscoelastic flow in the lower crust, with only a handful of near-field GPS stations and InSAR picking up the afterslip (Fig. 3). For this data set, the statistical algorithm has preferentially matched the InSAR observations above the GPS when balancing constraints and slight inconsistencies between the data sources, a known issue with joint inversions (27). This is due to more significant errors associated with GPS, especially the vertical component, which has a large scatter [root mean square (RMS) of 10 mm or more] and is, in general, less suitable for use over this time period. In particular, we misfit some of the GPS stations and InSAR points located closest to the Futagawa fault, likely as a result of localized shallow deformation not included in our model such as surface cracks or poroelastic rebound.

Fig. 3 Observed and modeled postseismic velocities.

(A) Unwrapped Sentinel-1A line-of-sight (LOS) InSAR velocity composite from three acquisitions between 20 April 2016 and 7 June 2016. GPS stations are indicated with black triangles, known faults are mapped with red lines, LOS vectors are indicated with orange arrows, and volcanoes are shown with red triangles. (B) LOS postseismic velocity model, with contributions from afterslip and ductile flow in lower-crustal strain volumes. The model captures the key features of the observations, excepting some high-frequency fluctuations, especially those following the fault traces of the Hinagu and Futagawa faults. (C) Postseismic observed and modeled GPS velocities. Observational data are shown with white arrows and recovered velocities are in black. Contributions due to afterslip on the Hinagu and Futagawa faults and ductile flow in the lower-crustal strain volumes are shown in red and blue, respectively. Subareal background color scale represents modeled uplift or subsidence velocity across Kyushu, where colored circles represent observed vertical velocity at GPS stations.

The viscosity of rocks depends on many physical parameters including local stress, strain rate, temperature, water content, composition, and degree of metamorphism, leading to tremendous uncertainties in the strength of the lower crust. We estimated the effective transient viscosity of the lower crust beneath Kyushu from the inferred kinematics usingEmbedded Imagewhere Embedded Image and Embedded Image are the effective stress and strain rate in the first hours after the earthquake (Fig. 4 and fig. S12), obtained from the second invariant of the respective tensors. We included the preearthquake stress assuming a steady-state uniform viscosity of 1019 Pa·s at a background strain rate of 10−15 s−1 (0.03 μstrains/year). We obtained effective viscosities in the range between 5 × 1016 and 1019 Pa·s, with the lowest viscosity beneath Mount Aso and Mount Kuju; this value is in reasonable agreement with that found in another volcanic arc beneath Santorini volcano, Greece, resulting from thermal activation (28, 29) and commensurate with predictions for transient creep of a thermally activated nonlinear rheology (6). A similar low-viscosity region is not found beneath Mount Unzen which, unlike the other volcanoes, is related to the opening of the Okinawa trench. The highest values are controlled by our assumptions for prestress. The low viscosities we recovered require an elevated geothermal gradient and also hint at the existence of melt pockets in the crust beneath these volcanoes. This low-viscosity region would also damp the rupture propagation and explain why the coseismic rupture terminated here. The large anelastic, predominantly normal, strain rate beneath Mount Aso may also have helped induce the eruption on 8 October 2016. The broader regions of low viscosity follow the stress contour, compatible with activation of dislocation creep or transient creep. We interpret the low transient viscosities as the result of a combination of steady-state and transient creep mechanisms. The low effective viscosities in the northern quadrant correspond to a large quaternary plutonic body associated with a notable heat flow anomaly (30). We also find little deformation within the mantle wedge (~10 μstrain/year), which would likely require a larger stress perturbation, such as the one created by the Mw 8.6 Indian Ocean (6) or Mw 9.0 Tōhoku earthquakes (7).

Fig. 4 Transient viscosity structure of the lower crust.

Transient viscosity of the lower crust with volcanoes marked in red triangles, the Hinagu and Futagawa faults in black, and coseismic stress contours in orange. The regions of low viscosity follow the pattern of coseismic stress change modulated by the distribution of arc volcanism and plutonic bodies in Kyushu, with noticeable low-viscosity anomalies beneath Mount Aso and Mount Kuju.

The postseismic deformation after the Kumamoto earthquake sequence yields insights into the distribution of brittle and ductile crustal processes beneath Kyushu, illuminating low transient effective viscosities in the lower crust, and thermally activated deformation beneath volcanic edifices. This provides a benchmark for the background variation in mechanical strength at a volcanic arc. The methodology introduced in this study reveals the potential for geodesy to illuminate the mechanics of distributed deformation and brings us one step closer to building a rheological model of the lithosphere-asthenosphere system.

Supplementary Materials

Materials and Methods

Figs. S1 to S12

References (3146)

References and Notes

  1. Materials and methods are available as supplementary materials.
  2. Acknowledgments: This research was supported by the National Research Foundation (NRF) of Singapore under the NRF Fellowship scheme (National Research Fellow Award NRF-NRFF2013-04), the Earth Observatory of Singapore under the Research Centres of Excellence initiative, and the Singapore Ministry of Education (Tier 2 project ARC5/14). We also thank Y. Ito for supporting access to the GEONET data. Data used in this study are stored at the Earth Observatory of Singapore and are freely available upon request.

Stay Connected to Science

Navigate This Article