Next Article in Journal
Tree Species Classification and Health Status Assessment for a Mixed Broadleaf-Conifer Forest with UAS Multispectral Imaging
Next Article in Special Issue
Ground Deformation Modelling of the 2020 Mw6.9 Samos Earthquake (Greece) Based on InSAR and GNSS Data
Previous Article in Journal
A New Method for Automated Clearcut Disturbance Detection in Mediterranean Coppice Forests Using Landsat Time Series
Previous Article in Special Issue
Towards a PS-InSAR Based Prediction Model for Building Collapse: Spatiotemporal Patterns of Vertical Surface Motion in Collapsed Building Areas—Case Study of Alexandria, Egypt
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multi-Segment Rupture Model of the 2016 Kumamoto Earthquake Revealed by InSAR and GPS Data

1
School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China
2
Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, Wuhan 430079, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2020, 12(22), 3721; https://doi.org/10.3390/rs12223721
Submission received: 5 October 2020 / Revised: 6 November 2020 / Accepted: 9 November 2020 / Published: 12 November 2020
(This article belongs to the Special Issue Ground Deformation Patterns Detection by InSAR and GNSS Techniques)

Abstract

:
The 2016 Kumamoto earthquake, including two large (Mw ≥ 6.0) foreshocks and an Mw 7.0 mainshock, occurred in the Hinagu and Futagawa fault zones in the middle of Kyushu island, Japan. Here, we obtain the complex coseismic deformation field associated with this earthquake from Advanced Land Observation Satellite-2 (ALOS-2) and Sentinel-1A Interferometric Synthetic Aperture Radar (InSAR) data. These InSAR data, in combination with available Global Positioning System (GPS) data, are then used to determine an optimal four-segment fault geometry with the jRi method, which considers both data misfit and the perturbation error from data noise. Our preferred slip distribution model indicates that the rupture is dominated by right-lateral strike-slip, with a significant normal slip component. The largest asperity is located on the northern segment of the Futagawa fault, with a maximum slip of 5.6 m at a 5–6 km depth. The estimated shallow slips along the Futagawa fault and northern Hinagu fault are consistent with the displacements of surface ruptures from the field investigation, suggesting a shallow slip deficit. The total geodetic moment release is estimated to be 4.89 × 1019 Nm (Mw 7.09), which is slightly larger than seismological estimates. The calculated static Coulomb stress changes induced by the preferred slip distribution model cannot completely explain the spatial distribution of aftershocks. Sensitivity analysis of Coulomb stress change implies that aftershocks in the stress shadow area may be driven by aseismic creep or triggered by dynamic stress transfer, requiring further investigation.

Graphical Abstract

1. Introduction

In April 2016, a series of shallow earthquakes struck Kumamoto prefecture of Kyushu island, southwest Japan, causing more than 100 fatalities, more than 1000 injuries, and severe destruction to houses and infrastructure. The sequence started with an Mw 6.2 foreshock on 14 April 2016 (UTC), and the Mw 7.0 mainshock, which is thought to be the largest event in central Kyushu island recorded since 1990, occurred approximately 28 h later. The two events generated strong ground motions that reached the maximum seismic intensity of 7 (the largest value on the Japan Meteorological Agency (JMA) scale). Previous studies related to the kinematic rupture process have shown that the mainshock rupture propagated unilaterally to the northeast direction and terminated near the southwest edge of the Aso volcano [1,2,3]. The moment tensor solutions from seismic catalogs (e.g., Global Centroid Moment Tensor (GCMT) [4], United States Geological Survey (USGS) [5], and National Research Institute for Earth Science and Disaster Resilience (NIED) F-net [6]) indicate that the mainshock was dominated by right-lateral motion with significant non-double-couple components, suggesting that rather complex mechanisms were involved in this sequence (Figure 1).
This sequence produced approximately 34 km-long surface ruptures, mainly along the Futagawa fault zone and its eastern extension to the caldera of the Aso volcano, and the northernmost part of the Hinagu fault zone [7]. The Futagawa and Hinagu fault zones belong to part of the western extension of the Median Tectonic Line (MTL), known as the longest right-lateral strike-slip active fault in Japan. These fault zones are located within the Beppu–Shimabara graben, which is dominantly controlled by north–south extension and some northeast–southwest shear modifications from MTL [8]. The oblique subduction of the Philippine Sea Plate (PHS) and the northern extension of the opening Okinawa Trough (OT) are considered to be the sources that drove the formation of the graben structure [9].
The spatial distribution of seismicity and surface ruptures from the field investigation suggest that multi-segment faults ruptured during the earthquake. To better understand the complicated rupture and stress interaction, it is crucial to clarify the source characteristics of this earthquake. Several previous studies [1,2,3,10,11,12,13,14,15] have reported the source models from the inversion of teleseismic, strong-motion, and geodetic datasets, either independently or jointly. However, the fault geometry and maximum slip of various models differ significantly. Typically, Kubo et al. [10] estimated a maximum slip of 3.8 m from strong-motion waveforms with a curved fault model. Their results fitted the waveforms relatively well, yet the maximum slip seemed to be underestimated when compared with the peak slip of 5.7 m from the joint inversion of strong-motion, teleseismic body, and surface waveforms [11]. Fukahata and Hashimoto [12] determined the fault dip angles and slip distribution on two-segment fault planes by inverting the ALOS-2 InSAR data. Their model was well-constrained, but a large misfit (up to 30 cm) occurred near the caldera of Aso volcano, which suggested that their fault model needed to be improved. Additionally, Yue et al. [3] performed the joint inversion of GPS, strong-motion, Synthetic Aperture Radar (SAR) images, and surface offset data with a curved multi-segment fault model. Their results showed the maximum slip exceeding 10 m, which is rare for an Mw 7.0 intraplate event. These studies only assumed northwest-dipping faults, while Ozawa et al. [13], Yoshida et al. [14], and Zhang et al. [15] indicated that a southeast dipping fault within the Aso caldera region might also have ruptured. Therefore, a further study of the fault geometry and slip distribution is required to better understand the source mechanism of this complex earthquake.
In this study, the complex deformation field associated with the 2016 Kumamoto earthquake is firstly derived from both ALOS-2 and Sentinel-1A SAR satellite images. Subsequently, a combination of InSAR and GPS data is utilized to determine the optimal multi-segment fault geometry with the jRi method and to invert the slip distribution on the fault planes. Finally, the Coulomb stress change calculation and its sensitivity analysis are carried out to investigate the possible triggering mechanisms.

2. Coseismic Surface Displacements from InSAR and GPS Data

2.1. InSAR Data

During the 2016 Kumamoto earthquake, both ALOS-2 and Sentinel-1 SAR satellites captured the coseismic surface deformation field covering the entire epicenter area. The ALOS-2 satellite, equipped with the Phased Array type L-band Synthetic Aperture Radar-2 (PALSAR-2) and operated by the Japan Aerospace Exploration Agency (JAXA), provides excellent surface displacement observations. The L-band ALOS-2 satellite can observe the ground surface from the right or left direction and penetrate the heavy vegetation to obtain a high coherent interferogram. The Sentinel-1 satellite is composed of a constellation of two satellites (A and B) maintained by the European Space Agency (ESA) and provides C-band SAR images to detect crustal deformation signals. This constellation operates the novel Terrain Observation with Progressive Scans (TOPS) mode for its standard acquisition, which has a broad swath of up to 400 km wide, with a spatial resolution down to 5 m and a minimum 12-day repeat cycle [16]. In this study, we collected one pair of images from the ALOS-2/PALSAR-2 descending track (P023D) and two pairs of images from both the descending and ascending track (T163D and T156A) of Sentinel-1A, in order to map the coseismic surface deformation. Detailed information on the SAR data used in this study is given in Table 1.
GAMMA software [17] was utilized to process all SAR images. The interferograms were multi-looked with ratios between the range and azimuth direction of 10:2 and 8:8 for the Sentinel-1A and ALOS-2/PALSAR-2 images, respectively. The 90-m-resolution Digital Elevation Model (DEM) from the Shuttle Radar Topography Mission (SRTM) [18] was used to remove the topographic effects. The interferograms were then filtered with an adaptive-filter method [19] and unwrapped using statistical-cost network-flow approaches [20]. Finally, all interferograms were geocoded to the World Geodetic System (WGS-84) geographic coordinates with a 90 m resolution.
The interferograms obtained from SAR images show significant deformation signals along the mapped active faults associated with the earthquake (Figure 2). Compared with the L-band ALOS-2 interferogram, two C-band Sentinel-1A interferograms can hardly acquire the deformation near the faults due to the worse decorrelation caused by the large displacement gradient, surface ruptures, and landslides. Intensive fringes appear on both the Futagawa fault and the northern part of the Hinagu fault, along which displacement discontinuities can be clearly identified. The unwrapped interferograms of two descending tracks (P023D and T163D) both show three major asymmetric lobes of deformation and detect the line-of-sight subsidence signals at the Aso caldera. Two lobes are located along the mapped Futagawa and Hinagu fault zones, while another is located within the Aso caldera. The observed line-of-sight (LOS) displacements vary from –46 to 18 cm and from –44 to 67 cm for the P023D and T163D tracks, respectively. The slight discrepancy between the two descending tracks arises from the different incidence angles of ALOS-2 (~36°) and Sentinel-1A (~38°) satellite SAR sensors [22]. In the ascending track (T156A), the LOS displacement varies from –97 to 106 cm and the range is much larger than that of the two descending tracks. The LOS displacements on both sides of the Futagawa-Hinagu fault zone display opposite directions, consistent with a strike-slip event. Large displacements of up to 50 cm in the LOS direction are also detected within the Aso caldera.

2.2. GPS Data

In addition to InSAR data, coseismic GPS displacements near the epicentral region from the GPS Earth Observation Network system (GEONET) were also collected. The distribution of GPS stations is shown in Figure 3. The horizontal displacements are consistent with the predominant right-lateral rupture, while the vertical displacements are much more complicated. The largest ~1 m horizontal displacement to the southwest and 0.25 m uplift can be observed at the 960701 station, which is located within kilometers of the mapped fault trace near the Aso caldera. However, this station was not included in the inversion because it is located too close to the fault trace, and the displacement could not be properly modeled by Green’s functions [3]. Moreover, the small displacements with large errors were also not used in the inversion. Consequently, 45 horizontal and six vertical displacements were used to constrain the fault geometry and slip distribution. It is widely known that InSAR measures large-scale line-of-sight displacements with a high spatial resolution, whereas GPS can provide 3-D continuous displacements with a sub-centimeter accuracy. Therefore, a combination of InSAR and GPS datasets can better constrain the fault geometry and slip distribution associated with earthquakes.

3. Model

The surface displacements derived from InSAR and GPS data were modeled using the rectangular dislocations embedded in a homogeneous, isotropic, and elastic half-space [23] and assuming a Poisson’s ratio of 0.25. Millions of highly correlated data points in the unwrapped interferograms were down-sampled with a quadtree sampling method [24] to improve the computational feasibility and efficiency of the model. The InSAR and GPS datasets were weighted equally in the joint inversion, which could fit both datasets quite well, as shown by the relative weight ratio tests (Figure S1).

3.1. Determination of the Fault Geometry with the jRi Method

The surface ruptures, active fault traces, InSAR deformation maps, and seismicity indicate that multi-segment fault planes ruptured during the earthquake. Therefore, four-segment fault planes (hereafter termed F1, F2, F3, and FH) are considered to represent the rupture planes, among which segment F1 is set as a southeast dipping fault plane in the eastward extension of the Futagawa fault zone under the Aso caldera, and the other three northwest dipping segments are located along the Futagawa and northern Hinagu fault zones (Figure 1).
However, determining the fault geometry parameters of multiple segments and the slip distribution on the fault planes simultaneously is computationally expensive [25,26]. In this study, the location and extent of the four-segment fault planes could be identified from geologically-mapped fault traces, surface ruptures, and the seismicity distribution, but the fault dip angles of those segments were difficult to determine, especially for the F1 segment located at a seismic gap of the aftershock. Therefore, we considered fixing the location and extent of the four-segment faults and then searching for the optimal fault dip angles that minimized the data misfit between observed and predicted surface displacements. Regularization errors, such as weight root mean square errors (WRMSE) with the form of C d 1 / 2 ( d G ( δ ) s ) 2 or root mean square errors (RMSE) with the form of ( d G ( δ ) s ) 2 , are the most common measures of misfit [25,27,28], and the jRi value [29], including both regularization error and perturbation error from data noise, has also been proposed to find the best-fitting model [30]. The jRi value with respect to fault dip angles can be expressed as:
jRi ( δ ) = C d 1 / 2 ( d G ( δ ) s ) 2 1 k d i a g ( M I 2 M T ) + 1 k d i a g ( M M T ) ,
where C d 1 / 2 is the weight matrix from Cholesky decomposition of the inverse of the data covariance matrix C d ; G ( δ ) is the assembled matrix of Green’s functions for a vector of dip angle δ (with δ = [ δ 1 ,   ,   δ n ] ,   n = 4 ), relating the slip s on a dislocation to the observed data d (i.e., InSAR and GPS datasets); k is the amount of data used in the inversion; I 2 is a matrix consisting of 2 × 2   I matrices, where I represents a k-by-k identity matrix; and M = [ I , N ] , with N = G ( δ ) G g . G g is the generalized inverse with the following form:
G g = ( [ C d 1 / 2 G ( δ ) κ L ]   T [ C d 1 / 2 G ( δ ) κ L ] ) 1 ( C d 1 / 2 G ( δ ) ) T ,
where L is the Laplacian smoothing matrix and κ is a regularization factor that determines the slip roughness L s 2 . The optimal κ = 2 can be chosen from the trade-off curve [25] between the data misfit and the average slip roughness or the jRi curve [29], as shown in Figure S2.
It has been demonstrated that the determination of the fault dip angle via the jRi method with geodetic data is more robust than the performance of WRMSE or RMSE in simple one-segment fault scenarios, such as the 2015 Mw 7.9 Gorkha Nepal earthquake [30], but further verification was required for the complex multi-segment fault scenario in this study. Therefore, we extended the jRi method to determine the multi-segment fault dip angles and compared the results with the RMSE method to illustrate the effectiveness of the jRi method. More specifically, the optimal fault dip angles of four segments were iteratively searched, while the slip distributions on the four-segment faults were inverted simultaneously, until the jRi or RMSE value exhibited an insignificant difference. In order to avoid the effect of the imprecise estimation of the slip distribution on the fault planes during searching for optimal dip angles, each fault segment was further discretized into 2 km × 2 km sub-faults.
Figure 4 shows the result of estimating the fault dip angles with the jRi and RMSE methods from InSAR and GPS data. The optimal fault dip angles determined with the jRi method were estimated to be 77°, 57°, 63°, and 74° for the F1, F2, F3, and FH segments, respectively. The fault dip angles of F2, F3, and FH segments were also well-determined with the RMSE method, and only slight differences can be observed when compared with the result from the jRi method. However, the fault dip angle of F1 was poorly determined with the RMSE method, resulting in a shallower dip of 68° compared with the result obtained from the jRi method. The detailed fault parameters of the four-segment fault model used in this study are summarized in Table 2.

3.2. Distributed Slip Model

After the fault geometry parameters of the four segments were determined, a bounded variable least squares method [31] was adopted to invert the distributed slip on each segment. Regarding the limitation on slip, all segments were allowed to include both right-lateral strike-slip and free dip-slip components with a magnitude of less than 10 m, which is a reasonable constraint based on the focal mechanism of the mainshock and previous studies [1,2,6]. Moreover, the Laplacian smoothing constraint was also adopted to avoid non-physical oscillation of the slip distribution.
The preferred slip distribution model reveals a single asperity on each of the four segments, with a maximum slip of 5.6 m on the F2 segment at a 5–6 km depth (Figure 5). The predominant slip occurs at a depth of 0–13 km, which confirms the shallow hypocenters of this earthquake sequence. The relocated aftershocks are mostly located at a depth of between 6 and 15 km, where the slip decays quickly. The slip on the F2 segment is mainly right-lateral strike-slip, but with a significant normal slip of ~2.7 m. The shallowest slip on the F2 segment is characterized by 1–2 m strike-slip, together with 0–1.5 m normal slip, which is consistent with the lateral and vertical displacement components of surface ruptures from the field investigation [7]. On the F3 and FH segments, the slip is almost pure right-lateral, but the maximum slip on the F3 segment (3.3 m) is slightly larger than that on the FH segment (3.0 m). The slip on the southeast dipping F1 segment exhibits a dominant normal slip, but smaller magnitude than the other three segments, and the maximum slip is less than 2.5 m. The total geodetic moment release was estimated to be 4.89 × 1019 Nm, which is equivalent to a Mw 7.09 event, assuming a shear modulus of 32 GPa. The estimated moment release is in general agreement with that of Zhang et al. [15], but slightly larger than the seismic moment derived from teleseismic data [4,5,6] or regional strong-motion data [2]. Note that InSAR data cover the largest foreshock (Mw 6.0) and four days of aftershocks, which may explain the slightly higher moment release.
Figure 6 shows a comparison of observed and predicted InSAR data from our best-fitting slip model, as well as the corresponding residuals. The general deformation patterns from both ALOS-2 and Sentinel-1A observations can be interpreted well by the preferred model. However, there are some near-fault residuals for ALOS-2 descending track P023D with maximum misfit of up to 10 cm, which may be attributed to a combination of complicated near-field deformation and possible nontectonic influences of liquefaction, landslides, and inelastic deformation [32]. The histograms of the residual distribution of each InSAR track are illustrated in Figure 7. The RMSEs were estimated to be 5.64, 2,75, and 2.74 cm for P023D, T163D, and T156A tracks, respectively. The preferred model provides a good fit to the GPS horizontal and vertical displacements shown in Figure 3. A detailed comparison of three components of the observed and predicted GPS displacements is shown in Figure 7. The fits to horizontal displacements are relatively better than those to vertical displacements, with RMSEs of 1.70, 2.35, and 4.38 cm for east, north, and vertical components, respectively.

4. Coulomb Stress Change Analysis

Numerous previous studies have shown that the Coulomb stress changes caused by an earthquake may trigger subsequent events, change the seismicity rate, and affect the state of nearby faults through stress transfer [33,34,35]. Based on the Coulomb failure criterion, the change in Coulomb failure stress (ΔCFS) due to an earthquake on the receiver faults can be calculated as follows:
Δ CFS = Δ τ s + μ Δ σ n ,
where Δ τ s is the shear stress change parallel to the slip direction; μ is the effective friction coefficient, varying between 0 and 0.8, with a moderate value of 0.4 used as default; and Δ σ n is the normal stress change (positive for tension). Positive values of Δ CFS on the receiver faults with a typical triggering threshold of ~0.1 bar (0.01 MPa) are more likely to promote failure and increase regional seismicity, while negative ones may inhabit failure and decrease seismicity [36].
In order to investigate the possible triggering mechanism of the sequence, the static stress changes were calculated based on the distributed slip model in this study. The calculations were performed using Coulomb 3.4 software [37], with an effective friction coefficient of 0.4 and a shear modulus of 32 GPa. Receiver fault parameters (i.e., strike, dip, and rake) were specified as 220°, 68°, and −165°, respectively, and were derived from the equivalent focal mechanism of our preferred model, which shows a slight discrepancy with focal mechanisms of GCMT [4], USGS [5], and NIED [6] F-net (Table S1). Static stress changes at various depths between 0 and 20 km were compared with the locations of Mw > 3.0 aftershocks that occurred before 30 April 2016. Only large and early aftershocks were selected in our analysis to reduce possible location errors and multiple interactions between aftershocks.
Figure 8 shows the static Coulomb stress, shear stress, and normal stress changes at different depths with an interval of 5 km, compared with the spatial distribution of Mw > 3.0 aftershocks. The localized pattern of Coulomb stress changes is significantly observed in the vicinity of fault zones and dominated by shear stress changes. At a depth of 2.5 km, the calculated Coulomb stress increases larger than 10 bars are mainly located near the most northeastern and southwestern ends of the fault zone, where aftershocks primarily occurred. The pattern of stress changes at a 7.5 km depth is similar to that at a 2.5 km depth, but is more complicated at a 12.5 km depth. The area with positive Coulomb stress changes at a 12.5 km depth is larger than that at a shallower depth, but the majority of aftershocks do not concentrate on the Coulomb stress increase regions. As the depth increases to 17.5 km, the positive ΔCFS area has only a few aftershocks and may have a potential seismic risk in the future. Additionally, if focal mechanisms are available, a more stringent calculation that resolves Coulomb stress changes on the nodal planes associated with the aftershocks can be conducted to assess whether or not the nodal planes are brought closer to failure. The focal mechanisms of aftershocks with Mw >5.0 from NIED F-net were selected to calculate Coulomb stress changes on the two conjugate planes at each hypocenter. The result highlights that six out of the seven aftershocks are positively affected by a Coulomb stress increase on at least one conjugate plane (Table 3). Intriguingly, the largest aftershock (event 1) occurred just 20 min after the mainshock, but the Coulomb stress changes resolved on both conjugate planes show negative values. In general, static Coulomb stress triggering can explain the majority of Mw >5.0 aftershocks, but there are also some Mw >3.0 aftershocks within the stress shadow zone that may be triggered by other factors, which will be discussed in the following section.

5. Discussion

5.1. Seismogenic Fault Geometry

The fault model proposed in this study is composed of four-segment fault planes, including three northwest dipping planes (F2, F3, and FH) along the Futagawa and northern Hinagu faults and one southeast dipping plane (F1) along the northeastward extension of the Futagawa fault under the Aso caldera. The FH segment along the Hinagu fault zone has a strike of N205° and an estimated dip angle of 75°, in good agreement with other previous models [2,12]. The F3 segment with a strike of N226° connects the Futagawa and Hinagu faults, while the F2 segment has a strike of N236° along the Futagawa fault. The dip angles were estimated to be 63° and 57° for the F3 segment and F2 segment, respectively. The F1 segment with a strike of N40° and an optimal dip angle of 77° is located at the northeastward extension of the Futagawa fault, where dextral displacements exceeding 100 cm were observed in the field investigation [7]. Our preferred fault model is generally consistent with previous studies based on surface displacements and/or strong-motion waveforms [12,15]. For example, Fukahata and Hashimoto [12] obtained the dip angles of two faults, attaining values of 61° ± 6° for the Futagawa fault and 74° ± 12° for the Hinagu fault by inverting InSAR data. Zhang et al. [15] combined InSAR, GPS, and strong-motion data to construct a three-segment model, in which the Futagawa fault segment dips 60° to the northwest, the Hinagu fault segment dips 70° to the northwest, and a southeast dipping fault segment near the Aso volcano has a dip angle of 80°.
However, there are several studies proposing northwest dipping fault models in the Aso caldera region based on teleseismic and/or strong-motion data [10,38]. Is a southeast dipping fault model more reasonable than a northwest dipping one? We carried out additional tests to compare the RMSE and jRi values for the southeast dipping fault F1 and the northwest dipping one. The results (Figure S3) indicate that a southeast dipping fault F1 can better fit the observed data than a northwest dipping one. Ozawa et al. [13] suggested that a southeast dipping fault plane near the Aso caldera region should be considered to explain the complex deformation pattern obtained from SAR images. Yoshida et al. [14] also confirmed that the model with a southeast dipping fault predicted the observed subsidence around the central cones of Aso volcano better than a northwest dipping one. The vertical displacements calculated from our preferred model were utilized to confirm the predictive performance by comparing the data with those of previous models. The result (Figure S4) indicates that our four-segment model with a southeast dipping F1 segment can better explain the subsidence around the Aso volcano from the decomposition of InSAR data [3,39].

5.2. Comparison with Previous Slip Distribution Models

The checkerboard and jackknife tests (Appendix A) show that the slip distribution model on each segment is well-constrained by the combination of InSAR and GPS data. Our preferred model is in general agreement with previous models based on teleseismic waveforms [1], strong-motion data [2,10], geodetic data [39,40], or combinations of them [3,11,15]. Consensus has been reached on at least three things, which are: (1) multi-segment faults with varying fault geometries ruptured during this earthquake; (2) the ruptures initiated at the Hinagu segment and propagated unilaterally to the northeastern extension of the Futagawa fault under the Aso caldera; and (3) the slip on the Hinagu fault (FH) was almost purely right-lateral strike-slip, and the Futagawa fault (F2 and F3) was dominated by right-lateral strike-slip, along with a significant normal slip component. However, there are some inconsistencies in the maximum slip and slip pattern. Our preferred slip distribution model shows that each segment has a concentrated slip area (asperity) and the largest asperity is located on the F2 segment, with a peak slip of 5.6 m at a 5–6 km depth. The maximum slip (5.6 m) is similar to the results of [1,2,11,12], but smaller than those of [3] (~10 m), [15] (7.4 m), and [41] (6.9 m). The slip pattern discrepancy may be attributed to the different data sets, model regularizations, and parameterizations employed. For example, Yue et al. [3] combined GPS, strong-motion, InSAR data from both ALOS-2 and Sentinel-1A images, and surface offsets to invert the kinematic rupture process on the curved multi-segment fault planes. They selected a smaller smoothing factor based on the consistency between the calculated Coulomb stress and aftershock distribution, which might overestimate the peak slip. Using Sentinel-1A imagery, GPS, and strong motion data, Zhang et al. [15] also inverted the slip distribution on the three-segment fault planes. Compared with our model, Zhang et al. [15] constructed the fault model with a similar southeast-dipping fault near the Aso volcano and generated a rougher slip pattern, choosing a smaller smoothing factor with the L-curve method. In our model, the Futagawa fault is parameterized into two segments (F2 and F3) and has two different asperities, which are consistent with the third and fourth asperity locations in Yue et al.’s model [3]. In summary, based on adequate constraints of InSAR (both ALOS-2 and Sentinel-1A images) and GPS datasets, our preferred model selects an optimal smoothing factor with the jRi method and presents a more plausible slip pattern.

5.3. Shallow Slip Deficit

The estimated slips on the shallowest sub-faults along the Futagawa and northern Hinagu faults are consistent with the displacements of surface ruptures from the field investigation [7] (Figure 5). On the F2 segment, the estimated slip has a peak of 5.6 m at a 5–6 km depth and diminishes to an average of 2 m near the surface, which implies a shallow slip deficit (SSD) of ~65% (Figure 5). A similar SSD also occurred on the other three segments. Previous studies have demonstrated that SSD might be caused by velocity-strengthening friction at a shallow depth that can be recovered through shallow afterslip [42,43] or distributed coseismic brittle damage in the uppermost crust layer [44,45]. Moreover, some studies have shown that SSD can also arise from artifacts due to a lack of near-fault data [46] or from model uncertainties [47]. However, the SSD in this study is real because the slip near the surface was estimated well in the inversion, and a similar SSD was also obtained using near-fault data in several previous studies [3,41]. Interestingly, Yue at al. [3] and Moore et al. [48] both observed significant post-seismic deformation close to the faults, which might arise from afterslip at a shallow depth. Therefore, velocity-strengthening friction at shallow depths is more likely to be the mechanism that produces significant SSD, but the “distributed brittle failure” explanation cannot be ruled out from the currently available observations.

5.4. Sensitivity Analysis of Coulomb Stress Change

Although the majority of Mw > 5.0 aftershocks with available focal mechanisms support static Coulomb stress triggering, the entire aftershock sequence cannot be explained well by the calculated Coulomb stress changes on the receiver fault that is parallel to the seismogenic fault of the mainshock. A few aftershocks occurred in the stress shadow regions characterized by negative stress changes. Given that the calculation of Coulomb stress change involves three variables ( Δ τ s , Δ σ n , and μ in Equation (3)), the sensitivity of Coulomb stress calculations to the effective coefficient of friction, receiver fault, and source model also needs to be investigated. The effective coefficient of friction varies from 0.2 to 0.8 for different types of fault. The mainshock focal mechanisms from NIED, GCMT, USGS W-phase, and this study were adopted to define receiver faults in four different scenarios, respectively. Furthermore, six available finite fault slip models inverted from different datasets shown in Table S1 were collected to address the sensitivity of the source model. The result of the sensitivity tests (Figures S5–S7) indicates that the source model plays a more important role in the Coulomb stress calculations, especially in the case of a complex rupture event, while the receiver fault geometry and effective friction coefficient have a slight effect on the size of the Coulomb stress change area. A comparison of Coulomb stress changes calculated with different source models shows that the differences are apparent near the faults, and source models (i.e., Yue et al. [3], Zhang et al. [15], and our study) inverted from various datasets jointly can generate a more detailed stress change pattern which can better explain the spatial distribution of seismicity. In all scenarios, there are some common features in that the static Coulomb stress increase significantly triggers aftershocks located in the northern part of the Aso volcano, and a few aftershocks in a stress shadow area may be driven by aseismic creep [3,49] or triggered by dynamic stress transfer [50], which requires further investigation.

6. Conclusions

The complex coseismic deformation field associated with the 2016 Kumamoto earthquake sequence is derived from Sentinel-1A and ALOS-2 SAR images. Together with GPS data, InSAR data provide powerful constraints on the fault geometry and slip distribution of the seismogenic fault. Based on surface ruptures, active fault traces, InSAR deformation maps, and seismicity, we constructed a four-segment fault model with one southeast-dipping fault plane and three northwest-dipping fault planes. The fault dip angles of the four segments were determined well with the jRi method, which exhibited a better performance than the conventional RMSE method. Our preferred slip distribution model reveals a single asperity on each segment with a maximum slip of 5.6 m on the F2 segment and suggests a significant shallow slip deficit that may be caused by velocity-strengthening friction at shallow depths. Areas with large slips are spatially complementary to the distribution of aftershocks. The total geodetic moment release was estimated to be 4.89 × 1019 Nm, equivalent to an Mw 7.09 event. The Coulomb stress calculation and its sensitivity analysis indicate that static Coulomb stress triggering can only explain part of the aftershock sequence, while aftershocks that occurred in a stress shadow area may be triggered by other factors, such as aseismic creep or dynamic stress transfer.

Supplementary Materials

The following are available online at https://www.mdpi.com/2072-4292/12/22/3721/s1: Figure S1: The choice of the relative weight ratio of GPS relative to InSAR datasets; Figure S2: The smoothing factor κ determination; Figure S3. Comparison of the RMSE and jRi values between the southeast dipping fault F1 and the northwest dipping fault F1; Figure S4: Predicted vertical displacements from our preferred model and previous source models; Figure S5: Sensitivity test of Coulomb stress changes to effective coefficients of friction ( μ = 0.2 ,   0.4 ,   0.6 ,   and   0.8 ) at different depths; Figure S6: Sensitivity test of Coulomb stress changes to receiver faults at different depths, calculated from the preferred model with μ = 0.4 ; Figure S7: Sensitivity test of Coulomb stress changes to source models at different depths, resolved on the receiver fault of 220°/68°/–165° with μ = 0.4 ; Figure S8: Checkerboard test for resolution of the joint or individual inversion of InSAR and GPS dat;. Figure S9: Jackknife test results of a. the mean slip, b. the standard derivation (Std.), and c. the coefficient of variation (CV) of slip; Table S1: Source parameters of the 2016 Kumamoto earthquake obtained from the inversion of various datasets.

Author Contributions

T.C. and Z.H. conceived the work. Z.H. performed the main experiments and wrote the manuscript. T.C., M.W. and Y.L. analyzed the results. All authors contributed to the reviewing and editing of the paper. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Key Research and Development Program of China (Grant no. 2018YFC1503605 and 2016YFC1401506).

Acknowledgments

We would like to thank the Editor and three anonymous reviewers for their helpful suggestions and comments, which greatly improved the manuscript. We thank Yangmao Wen for helping with InSAR data processing and Yukitoshi Fukahata for providing GPS data. This work benefitted from discussions with Yuqing He and Chendong Zhang. The Sentinel-1A SAR data were provided by the European Space Agency (ESA) through Sentinels Scientific Data Hub, and the ALOS-2/ PALSAR-2 data were provided by JAXA through the RA6 project (ID: 3048). Most figures in this paper were prepared by the Generic Mapping Tools (GMT) software [51].

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Checkerboard and Jackknife Tests

A set of checkerboard tests were performed to investigate the spatial resolution of the slip distribution inverted from InSAR and GPS data, both individually and jointly. The synthetic slip model, whose sub-fault size increases from 2 to 4 km along the downdip direction, has an alternating slip of 0 or 2 m. We first generated simulated InSAR and GPS data, without adding any noise, as we focused on the contribution of each dataset in the inversion. Then, the simulated data were inverted to see how well the synthetic slip model was recovered. The result (Figure S8) shows that InSAR data have a great importance in recovering the slip distribution on each fault segment. GPS data have a poor model resolution due to the sparse distribution and a small number of stations, but still place valuable constraints on the peak slip. The joint inversion presents an improved spatial resolution of the slip distribution compared with individual InSAR or GPS inversion. The slip distribution from joint inversion is well-recovered for the shallow sub-faults (<10 km along the dip direction), but smeared at a deeper depth, especially on the FH segment.
In order to better quantify the robustness of the slip inversion, a jackknife test was conducted by randomly removing 20% of the data and running the inversion with the same model parameters [52]. This process was repeated 200 times, and the mean slip model, the standard deviation of slip, and the coefficient of variation (CV) were then calculated. The CV given by dividing the mean by the standard deviation is a good measure of the slip variability on each sub-fault. Sub-faults with a low CV have a more reliable slip than those with a high CV. The result (Figure S9) shows that the main slip area on each segment is stable due to their CV being lower than 0.2, while a high CV occurs at the edge of each segment.

References

  1. Yagi, Y.; Okuwaki, R.; Enescu, B.; Kasahara, A.; Miyakawa, A.; Otsubo, M. Rupture process of the 2016 Kumamoto earthquake in relation to the thermal structure around Aso volcano. Earth Planets Space 2016, 68, 118. [Google Scholar] [CrossRef] [Green Version]
  2. Asano, K.; Iwata, T. Source rupture processes of the foreshock and mainshock in the 2016 Kumamoto earthquake sequence estimated from the kinematic waveform inversion of strong motion data. Earth Planets Space 2016, 68, 147. [Google Scholar] [CrossRef] [Green Version]
  3. Yue, H.; Ross, Z.E.; Liang, C.; Michel, S.; Fattahi, H.; Fielding, E.; Moore, A.; Liu, Z.; Jia, B. The 2016 Kumamoto Mw = 7.0 Earthquake: A significant event in a fault–volcano system. J. Geophys. Res. Solid Earth 2017, 122, 9166–9183. [Google Scholar] [CrossRef]
  4. Global Centroid Moment Tensor Catalogue (GCMT). Available online: http://www.globalcmt.org/CMTsearch.html (accessed on 30 June 2020).
  5. United States Geological Survey (USGS). Available online: https://earthquake.usgs.gov/earthquakes/eventpage/us20005iis/moment-tensor (accessed on 30 June 2020).
  6. National Research Institute for Earth Science and Disaster Resilience (NIED). Available online: https://www.fnet.bosai.go.jp/event/tdmt.php?_id=20160415162400&LANG=en (accessed on 30 June 2020).
  7. Shirahama, Y.; Yoshimi, M.; Awata, Y.; Maruyama, T.; Azuma, T.; Miyashita, Y.; Mori, H.; Imanishi, K.; Takeda, N.; Ochi, T.; et al. Characteristics of the surface ruptures associated with the 2016 Kumamoto earthquake sequence, central Kyushu, Japan. Earth Planets Space 2016, 68, 191. [Google Scholar] [CrossRef] [Green Version]
  8. Matsumoto, S.; Nakao, S.; Ohkura, T.; Miyazaki, M.; Shimizu, H.; Abe, Y.; Inoue, H.; Nakamoto, M.; Yoshikawa, S.; Yamashita, Y.; et al. Spatial heterogeneities in tectonic stress in Kyushu, Japan and their relation to a major shear zone Seismology. Earth Planets Space 2015, 67, 172. [Google Scholar] [CrossRef] [Green Version]
  9. Kato, A.; Nakamura, K.; Hiyama, Y. The 2016 Kumamoto earthquake sequence. Proc. Jpn. Acad. Ser. B Phys. Biol. Sci. 2016, 92, 358–371. [Google Scholar] [CrossRef] [Green Version]
  10. Kubo, H.; Suzuki, W.; Aoi, S.; Sekiguchi, H. Source rupture processes of the 2016 Kumamoto, Japan, earthquakes estimated from strong-motion waveforms. Earth Planets Space 2016, 68, 161. [Google Scholar] [CrossRef] [Green Version]
  11. Hao, J.; Ji, C.; Yao, Z. Slip history of the 2016 Mw 7.0 Kumamoto earthquake: Intraplate rupture in complex tectonic environment. Geophys. Res. Lett. 2017, 44, 743–750. [Google Scholar] [CrossRef]
  12. Fukahata, Y.; Hashimoto, M. Simultaneous estimation of the dip angles and slip distribution on the faults of the 2016 Kumamoto earthquake through a weak nonlinear inversion of InSAR data. Earth Planets Space 2016, 68, 204. [Google Scholar] [CrossRef] [Green Version]
  13. Ozawa, T.; Fujita, E.; Ueda, H. Crustal deformation associated with the 2016 Kumamoto Earthquake and its effect on the magma system of Aso volcano. Earth Planets Space 2016, 68, 186. [Google Scholar] [CrossRef] [Green Version]
  14. Yoshida, K.; Miyakoshi, K.; Somei, K.; Irikura, K. Source process of the 2016 Kumamoto earthquake (Mj7.3) inferred from kinematic inversion of strong-motion records. Earth Planets Space 2017, 69, 64. [Google Scholar] [CrossRef]
  15. Zhang, Y.; Shan, X.; Zhang, G.; Gong, W.; Liu, X.; Yin, H.; Zhao, D.; Wen, S.; Qu, C. Source model of the 2016 Kumamoto, Japan, earthquake constrained by InSAR, GPS, and strong-motion data: Fault slip under extensional stress. Bull. Seismol. Soc. Am. 2018, 108, 2675–2686. [Google Scholar] [CrossRef]
  16. Yague-Martinez, N.; Prats-Iraola, P.; Gonzalez, F.R.; Brcic, R.; Shau, R.; Geudtner, D.; Eineder, M.; Bamler, R. Interferometric Processing of Sentinel-1 TOPS Data. IEEE Trans. Geosci. Remote Sens. 2016, 54, 2220–2234. [Google Scholar] [CrossRef] [Green Version]
  17. Werner, C.; Wegmüller, U.; Strozzi, T.; Wiesmann, A. GAMMA SAR and interferometric processing software. In Proceedings of the ERS ENVISAT Symposium, Gothenburg, Sweden, 16–20 October 2000. [Google Scholar]
  18. Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L.; et al. The shuttle radar topography mission. Rev. Geophys. 2007, 45. [Google Scholar] [CrossRef] [Green Version]
  19. Goldstein, R.M.; Werner, C.L. Radar interferogram filtering for geophysical applications. Geophys. Res. Lett. 1998, 25, 4035–4038. [Google Scholar] [CrossRef] [Green Version]
  20. Chen, C.W.; Zebker, H.A. Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization. J. Opt. Soc. Am. A 2001, 18, 338. [Google Scholar] [CrossRef] [Green Version]
  21. He, P.; Wen, Y.; Xu, C.; Chen, Y. Complete three-dimensional near-field surface displacements from imaging geodesy techniques applied to the 2016 Kumamoto earthquake. Remote Sens. Environ. 2019, 232, 111321. [Google Scholar] [CrossRef]
  22. Wen, Y.; Xu, C.; Liu, Y.; Jiang, G. Deformation and source parameters of the 2015 Mw 6.5 earthquake in Pishan, Western China, from Sentinel-1A and ALOS-2 data. Remote Sens. 2016, 8, 134. [Google Scholar] [CrossRef] [Green Version]
  23. Okada, Y. Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 1992, 82, 1018–1040. [Google Scholar]
  24. Jónsson, S.; Zebker, H.; Segall, P.; Amelung, F. Fault slip distribution of the 1999 Mw 7.1 Hector Mine, California, earthquake, estimated from satellite radar and GPS measurements. Bull. Seismol. Soc. Am. 2002, 92, 1377–1389. [Google Scholar] [CrossRef]
  25. Burgmann, R.; Ayhan, M.E.; Fielding, E.J.; Wright, T.J.; McClusky, S.; Aktug, B.; Demir, C.; Lenk, O.; Turkezer, A. Deformation during the 12 November 1999 Duzce, Turkey, earthquake, from GPS and InSAR data. Bull. Seismol. Soc. Am. 2002, 92, 161–171. [Google Scholar] [CrossRef]
  26. Fukahata, Y.; Wright, T.J. A non-linear geodetic data inversion using ABIC for slip distribution on a fault with an unknown dip angle. Geophys. J. Int. 2008, 173, 353–364. [Google Scholar] [CrossRef] [Green Version]
  27. Chen, T.; Newman, A.V.; Feng, L.; Fritz, H.M. Slip distribution from the 1 April 2007 Solomon Islands earthquake: A unique image of near-trench rupture. Geophys. Res. Lett. 2009, 36, 6–11. [Google Scholar] [CrossRef] [Green Version]
  28. Xu, C.; Liu, Y.; Wen, Y.; Wang, R. Coseismic slip distribution of the 2008 Mw 7.9 Wenchuan earthquake from joint inversion of GPS and InSAR data. Bull. Seismol. Soc. Am. 2010, 100, 2736–2749. [Google Scholar] [CrossRef]
  29. Barnhart, W.D.; Lohman, R.B. Automated fault model discretization for inversions for coseismic slip distributions. J. Geophys. Res. Solid Earth 2010, 115, 1–17. [Google Scholar] [CrossRef]
  30. Luo, H.; Chen, T.; Xu, C.; Sha, H. Fault dip angle determination with the jRi criterion and coulomb stress changes associated with the 2015 Mw 7.9 Gorkha Nepal earthquake revealed by InSAR and GPS data. Tectonophysics 2017, 714, 55–61. [Google Scholar] [CrossRef]
  31. Stark, P.; Parker, R. Bounded-variable least-squares: An algorithm and applications. Comput. Stat. 1995, 1–13. [Google Scholar]
  32. Fujiwara, S.; Morishita, Y.; Nakano, T.; Kobayashi, T.; Yarai, H. Non-tectonic liquefaction-induced large surface displacements in the Aso Valley, Japan, caused by the 2016 Kumamoto earthquake, revealed by ALOS-2 SAR. Earth Planet. Sci. Lett. 2017, 474, 457–465. [Google Scholar] [CrossRef]
  33. Reasenberg, P.A.; Simpson, R.W. Response of regional seismicity to the static stress change produced by the Loma Prieta earthquake. Science 1992, 255, 1687–1690. [Google Scholar] [CrossRef]
  34. King, G.C.P.; Stein, R.S.; Lin, J. Static stress changes and the triggering of earthquakes. Bull. Seismol. Soc. Am. 1994, 84, 935–953. [Google Scholar]
  35. Freed, A.M. Earthquake triggering by static, dynamic, and postseismic stress transfer. Annu. Rev. Earth Planet. Sci. 2005, 33, 335–367. [Google Scholar] [CrossRef] [Green Version]
  36. Toda, S.; Stein, R.S.; Reasenberg, P.A.; Dieterich, J.H.; Yoshida, A. Stress transferred by the 1995 Mw = 6.9 Kobe, Japan, shock: Effect on aftershocks and future earthquake probabilities. J. Geophys. Res. Solid Earth 1998, 103, 24543–24565. [Google Scholar] [CrossRef] [Green Version]
  37. Toda, S.; Stein, R.S.; Sevilgen, V.; Lin, J. Coulomb 3.3 Graphic-Rich Deformation and Stress-Change Software for Earthquake, Tectonic, and Volcano Research and Teaching—User Guide; U.S. Geological Survey Open-File Report 2011–1060; U.S. Geological Survey: Reston, VA, USA, 2011; 63p.
  38. Kobayashi, H.; Koketsu, K.; Miyake, H. Rupture processes of the 2016 Kumamoto earthquake sequence: Causes for extreme ground motions. Geophys. Res. Lett. 2017, 44, 6002–6010. [Google Scholar] [CrossRef]
  39. Himematsu, Y.; Furuya, M. Fault source model for the 2016 Kumamoto earthquake sequence based on ALOS-2/PALSAR-2 pixel-offset data: Evidence for dynamic slip partitioning. Earth Planets Space 2016, 68, 169. [Google Scholar] [CrossRef] [Green Version]
  40. Jiang, H.; Feng, G.; Wang, T.; Bürgmann, R. Toward full exploitation of coherent and incoherent information in Sentinel-1 TOPS data for retrieving surface displacement: Application to the 2016 Kumamoto (Japan) earthquake. Geophys. Res. Lett. 2017, 44, 1758–1767. [Google Scholar] [CrossRef] [Green Version]
  41. Scott, C.; Champenois, J.; Klinger, Y.; Nissen, E.; Maruyama, T.; Chiba, T.; Arrowsmith, R. The 2016 M7 Kumamoto, Japan, earthquake slip field derived from a joint inversion of differential lidar topography, optical correlation, and InSAR surface displacements. Geophys. Res. Lett. 2019, 46, 6341–6351. [Google Scholar] [CrossRef] [Green Version]
  42. Tse, S.T.; Rice, J.R. Crustal earthquake instability in relation to the depth variation of frictional slip properties. J. Geophys. Res. 1986, 91, 9452. [Google Scholar] [CrossRef] [Green Version]
  43. Marone, C. Laboratory-derived friction laws and their application to seismic faulting. Annu. Rev. Earth Planet. Sci. 1998, 26, 643–696. [Google Scholar] [CrossRef] [Green Version]
  44. Fialko, Y.; Sandwell, D.; Simons, M.; Rosen, P. Three-dimensional deformation caused by the Bam, Iran, earthquake and the origin of shallow slip deficit. Nature 2005, 435, 295–299. [Google Scholar] [CrossRef]
  45. Kaneko, Y.; Fialko, Y. Shallow slip deficit due to large strike-slip earthquakes in dynamic rupture simulations with elasto-plastic off-fault response. Geophys. J. Int. 2011, 186, 1389–1403. [Google Scholar] [CrossRef] [Green Version]
  46. Xu, X.; Tong, X.; Sandwell, D.T.; Milliner, C.W.D.; Dolan, J.F.; Hollingsworth, J.; Leprince, S.; Ayoub, F. Refining the shallow slip deficit. Geophys. J. Int. 2016, 204, 1867–1886. [Google Scholar] [CrossRef] [Green Version]
  47. Huang, M.H.; Fielding, E.J.; Dickinson, H.; Sun, J.; Gonzalez-Ortega, J.A.; Freed, A.M.; Bürgmann, R. Fault geometry inversion and slip distribution of the 2010 Mw 7.2 El Mayor-Cucapah earthquake from geodetic data. J. Geophys. Res. Solid Earth 2017, 122, 607–621. [Google Scholar] [CrossRef] [Green Version]
  48. Moore, J.D.P.; Yu, H.; Tang, C.H.; Wang, T.; Barbot, S.; Peng, D.; Masuti, S.; Dauwels, J.; Hsu, Y.J.; Lambert, V.; et al. Imaging the distribution of transient viscosity after the 2016 Mw 7.1 Kumamoto earthquake. Science 2017, 356, 163–167. [Google Scholar] [CrossRef] [PubMed]
  49. Perfettini, H.; Avouac, J.P. Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks, application to the Chi-Chi earthquake, Taiwan. J. Geophys. Res. Solid Earth 2004, 109. [Google Scholar] [CrossRef]
  50. Nissen, E.; Elliott, J.R.; Sloan, R.A.; Craig, T.J.; Funning, G.J.; Hutko, A.; Parsons, B.E.; Wright, T.J. Limitations of rupture forecasting exposed by instantaneously triggered earthquake doublet. Nat. Geosci. 2016, 9, 330–336. [Google Scholar] [CrossRef]
  51. Wessel, P.; Smith, W.H.F.; Scharroo, R.; Luis, J.; Wobbe, F. Generic mapping tools: Improved version released. Eos Trans. Am. Geophys. Union 2013, 94, 409–410. [Google Scholar] [CrossRef] [Green Version]
  52. Melgar, D.; Ganas, A.; Geng, J.; Liang, C.; Fielding, E.J.; Kassaras, I. Source characteristics of the 2015 Mw 6.5 Lefkada, Greece, strike-slip earthquake. J. Geophys. Res. Solid Earth 2017, 122, 2260–2273. [Google Scholar]
Figure 1. Tectonic setting of the 2016 Kumamoto earthquake, Kyushu island, in southwest Japan. The upper-left inset shows the seismogenic background of southwest Japan, while the blue rectangle delineates the study area. The lower-left inset illustrates the magnitude time evolution of the first three-day earthquake sequence, and blue and red stars denote the largest foreshock and mainshock, respectively. Blue and red beach balls indicate the focal mechanism solutions of the foreshock and mainshock determined by NIED F-net [6], respectively, while gray beach balls show aftershocks with Mw >5.0. Moreover, the mainshock focal mechanisms from other earthquake catalogs (e.g., GCMT [4], USGS [5] W-phase, and body wave) are also shown for comparison. Color-filled circles denote the relocated foreshocks (blue) and aftershocks (the filled color varies with hypocentral depth) until 30 April and the circles’ size is proportional to the magnitude [3]. The red ellipse illustrates a ~10 km-long seismic gap of aftershocks near the Aso volcano. Active faults are shown with red lines. The gray squares denote the GPS stations. The red triangles show the active volcanoes and the purple square is Kumamoto city. Surface projections of the assumed fault planes (marked as F1, F2, F3, and FH) are represented by black dashed rectangles (the top edge marked by thick lines). PHS, Philippine Sea Plate; NT, Nankai Trough; OT, Okinawa Trough; MTL, Median Tectonic Line; EU, Eurasian plate; FF, Futagawa fault; HF, Hinagu fault.
Figure 1. Tectonic setting of the 2016 Kumamoto earthquake, Kyushu island, in southwest Japan. The upper-left inset shows the seismogenic background of southwest Japan, while the blue rectangle delineates the study area. The lower-left inset illustrates the magnitude time evolution of the first three-day earthquake sequence, and blue and red stars denote the largest foreshock and mainshock, respectively. Blue and red beach balls indicate the focal mechanism solutions of the foreshock and mainshock determined by NIED F-net [6], respectively, while gray beach balls show aftershocks with Mw >5.0. Moreover, the mainshock focal mechanisms from other earthquake catalogs (e.g., GCMT [4], USGS [5] W-phase, and body wave) are also shown for comparison. Color-filled circles denote the relocated foreshocks (blue) and aftershocks (the filled color varies with hypocentral depth) until 30 April and the circles’ size is proportional to the magnitude [3]. The red ellipse illustrates a ~10 km-long seismic gap of aftershocks near the Aso volcano. Active faults are shown with red lines. The gray squares denote the GPS stations. The red triangles show the active volcanoes and the purple square is Kumamoto city. Surface projections of the assumed fault planes (marked as F1, F2, F3, and FH) are represented by black dashed rectangles (the top edge marked by thick lines). PHS, Philippine Sea Plate; NT, Nankai Trough; OT, Okinawa Trough; MTL, Median Tectonic Line; EU, Eurasian plate; FF, Futagawa fault; HF, Hinagu fault.
Remotesensing 12 03721 g001
Figure 2. Wrapped and unwrapped InSAR maps. (a–c) Coseismic interferograms derived from the ALOS-2 and Sentinel-1A images. Each color cycle represents 11.8 cm of line-of-sight (LOS) deformation. (d–f) Unwrapped LOS deformation corresponding to (a–c). The active faults are shown as a red line. Yellow star and red triangle denote the mainshock epicenter and Aso volcano, respectively. Positive value indicates LOS displacement moving toward the satellite.
Figure 2. Wrapped and unwrapped InSAR maps. (a–c) Coseismic interferograms derived from the ALOS-2 and Sentinel-1A images. Each color cycle represents 11.8 cm of line-of-sight (LOS) deformation. (d–f) Unwrapped LOS deformation corresponding to (a–c). The active faults are shown as a red line. Yellow star and red triangle denote the mainshock epicenter and Aso volcano, respectively. Positive value indicates LOS displacement moving toward the satellite.
Remotesensing 12 03721 g002
Figure 3. Observed co-seismic GPS displacements compared with modeled displacements from the preferred model. (a) Horizontal displacements with exaggerated 95% confidence ellipses and (b) vertical displacements. Red star and red triangle denote the mainshock epicenter and Aso volcano, respectively.
Figure 3. Observed co-seismic GPS displacements compared with modeled displacements from the preferred model. (a) Horizontal displacements with exaggerated 95% confidence ellipses and (b) vertical displacements. Red star and red triangle denote the mainshock epicenter and Aso volcano, respectively.
Remotesensing 12 03721 g003
Figure 4. Comparison of the estimations of fault dip angles with the jRi and root mean square error (RMSE) methods. (a) Plot of the jRi values as a function of dip angles for four segments. (b) Plot of the RMSE values. Optimal dip angles corresponding to minimal jRi and RMSE values are marked with cross symbols.
Figure 4. Comparison of the estimations of fault dip angles with the jRi and root mean square error (RMSE) methods. (a) Plot of the jRi values as a function of dip angles for four segments. (b) Plot of the RMSE values. Optimal dip angles corresponding to minimal jRi and RMSE values are marked with cross symbols.
Remotesensing 12 03721 g004
Figure 5. Slip distribution of the preferred model inverted from InSAR and GPS data, and comparison with surface ruptures and aftershocks. (a–b) Surface projection of the slip distribution. The color-filled circles in (a) show the lateral displacements of surface ruptures from the field investigation (positive values for right-lateral displacements), while those in (b) show the vertical displacements (positive values for south-side uplift displacements). The red star denotes the epicenter, and red lines show the active faults. The white contours highlight the asperities with a slip magnitude of over 2 m. Four segments are marked with F1, F2, F3, and FH from northeast to southwest. (c) 3-D slip distribution viewing from SSE. The black arrows show the slip directions of each sub-fault. (d) The slip averaged at different depths on each segment and number of aftershocks versus depth. Slips on four segments are plotted with curves. The gray histogram illustrates the depth distribution of aftershocks.
Figure 5. Slip distribution of the preferred model inverted from InSAR and GPS data, and comparison with surface ruptures and aftershocks. (a–b) Surface projection of the slip distribution. The color-filled circles in (a) show the lateral displacements of surface ruptures from the field investigation (positive values for right-lateral displacements), while those in (b) show the vertical displacements (positive values for south-side uplift displacements). The red star denotes the epicenter, and red lines show the active faults. The white contours highlight the asperities with a slip magnitude of over 2 m. Four segments are marked with F1, F2, F3, and FH from northeast to southwest. (c) 3-D slip distribution viewing from SSE. The black arrows show the slip directions of each sub-fault. (d) The slip averaged at different depths on each segment and number of aftershocks versus depth. Slips on four segments are plotted with curves. The gray histogram illustrates the depth distribution of aftershocks.
Remotesensing 12 03721 g005
Figure 6. Observed and model-predicted InSAR data and the residuals. (a–c) Subsampled InSAR line-of-sight displacement data used in the inversion for ALOS-2 descending track P023D, Sentinel-1A descending track T163D, and ascending track T156A, respectively. (d–f) Model-predicted InSAR line-of-sight displacements corresponding to (a–c) from our preferred model. (g–i) Residuals after subtracting the model predictions from the observations. Positive value indicates the LOS displacement moving toward the satellite.
Figure 6. Observed and model-predicted InSAR data and the residuals. (a–c) Subsampled InSAR line-of-sight displacement data used in the inversion for ALOS-2 descending track P023D, Sentinel-1A descending track T163D, and ascending track T156A, respectively. (d–f) Model-predicted InSAR line-of-sight displacements corresponding to (a–c) from our preferred model. (g–i) Residuals after subtracting the model predictions from the observations. Positive value indicates the LOS displacement moving toward the satellite.
Remotesensing 12 03721 g006
Figure 7. Results of GPS and InSAR data fittings. (a–c) Comparison of the observations of GPS data and the predictions from the preferred model. (d–f) Histograms of residuals after subtracting the model predictions from the observations of InSAR data.
Figure 7. Results of GPS and InSAR data fittings. (a–c) Comparison of the observations of GPS data and the predictions from the preferred model. (d–f) Histograms of residuals after subtracting the model predictions from the observations of InSAR data.
Remotesensing 12 03721 g007
Figure 8. Static stress changes caused by the mainshock at different depths, calculated from the preferred model with μ = 0.4 . (a1d1) show the static Coulomb stress changes (positive ΔCFS may promote failure); (a2d2) show the shear stress changes (positive in the direction of slip); (a3d3) show the normal stress changes (positive when the fault is unclamped). The receiver fault is approximately parallel to the seismogenic fault of the mainshock, with strike, dip, and rake of (220°, 68°, and –165°, respectively. Aftershocks with Mw >3.0 are denoted by green dots in each panel, with hypocenter depth ranges of d ± 2.5 km. The white contours highlight the asperities with a slip amplitude of over 2 m. The red lines show the active faults. The black dashed rectangles are surface projections of fault planes.
Figure 8. Static stress changes caused by the mainshock at different depths, calculated from the preferred model with μ = 0.4 . (a1d1) show the static Coulomb stress changes (positive ΔCFS may promote failure); (a2d2) show the shear stress changes (positive in the direction of slip); (a3d3) show the normal stress changes (positive when the fault is unclamped). The receiver fault is approximately parallel to the seismogenic fault of the mainshock, with strike, dip, and rake of (220°, 68°, and –165°, respectively. Aftershocks with Mw >3.0 are denoted by green dots in each panel, with hypocenter depth ranges of d ± 2.5 km. The white contours highlight the asperities with a slip amplitude of over 2 m. The red lines show the active faults. The black dashed rectangles are surface projections of fault planes.
Remotesensing 12 03721 g008
Table 1. Synthetic aperture radar (SAR) interferometric pairs used in the study.
Table 1. Synthetic aperture radar (SAR) interferometric pairs used in the study.
TrackSensorOrbitModePrimary Image (yyyy/mm/dd)Secondary Image (yyyy/mm/dd) B
(m)
B T
(days)
σ
(cm)
Np
P023DALOS-2DescendingStrip-map2016/03/072016/04/1881.5422.3657
T163DSentinel-1ADescendingIW2016/03/272016/04/201.9242.5549
T156ASentinel-1AAscendingIW2016/04/082016/04/2065.4123.2568
B represents the perpendicular baseline between the orbits; B T is the interval between the master and slave images acquisition dates; Interferometric Wide (IW) swath mode is the main acquisition mode over land; σ, denoting measured error from a combination of various error sources (e.g., ionospheric and atmospheric effects, topographic correction residuals, and orbital and unwrapping errors), was estimated from the far-field area [21]; Np is the number of down-sampling points.
Table 2. Fault parameters of the four-segment fault model used in this study.
Table 2. Fault parameters of the four-segment fault model used in this study.
SegmentLatitude
(°)
Longitude
(°)
Depth
(km)
Length (km)Width (km)Strike (°)Dip (°)Max Slip (m)Moment
(Nm)
Mw
F1131.017832.90880.2121840772.34.89 × 10197.09
F2130.922232.83860.21218236575.6
F3130.821632.77220.21218226633.3
FH130.767232.69180.21618205743.0
The position (latitude, longitude, and depth) indicates the top-center of each fault segment. The rigidity was assumed to be 32 GPa. The fault dip angles were determined with the jRi method.
Table 3. Coulomb stress changes at the hypocenter resolved on the two conjugate nodal planes of the focal mechanisms of aftershocks with Mw >5.0.
Table 3. Coulomb stress changes at the hypocenter resolved on the two conjugate nodal planes of the focal mechanisms of aftershocks with Mw >5.0.
IDOrigin Time (UTC)HypocenterMwNodal Plane 1Nodal Plane 2
yyyy/mm/dddd:mmLatitude (°)Longitude (°)Depth (km)Strike/Dip/Rake Δ CFS 1
( bar )
Strike/Dip/Rake Δ CFS 2
( bar )
12016/04/1516:4532.8632130.899010.555.7286/35/–70−18.1581/57/–104–20.11
22016/04/1518:0332.9638131.08686.895.5209/60/–1747.35116/85/–307.79
32016/04/1518:5533.0265131.191010.895.5220/72/–1671.09126/78/–192.02
42016/04/1600:4832.8470130.835015.915.2230/38/–112−4.4377/55/–736.95
52016/04/1607:0232.6992130.720012.305.1255/30/–883.2472/60/–91–2.78
62016/04/1811:4133.0020131.19988.645.4314/86/32.38224/87/1761.07
72016/04/1908:5232.5352130.63539.965.3221/60/–1692.45126/81/–301.51
Coulomb stress changes were evaluated using our preferred slip distribution model and an efficient coefficient of friction of 0.4. The focal mechanisms originated from the F-net catalog of the National Research Institute for Earth Science and Disaster Resilience (NIED), Japan.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

He, Z.; Chen, T.; Wang, M.; Li, Y. Multi-Segment Rupture Model of the 2016 Kumamoto Earthquake Revealed by InSAR and GPS Data. Remote Sens. 2020, 12, 3721. https://doi.org/10.3390/rs12223721

AMA Style

He Z, Chen T, Wang M, Li Y. Multi-Segment Rupture Model of the 2016 Kumamoto Earthquake Revealed by InSAR and GPS Data. Remote Sensing. 2020; 12(22):3721. https://doi.org/10.3390/rs12223721

Chicago/Turabian Style

He, Zhongqiu, Ting Chen, Mingce Wang, and Yanchong Li. 2020. "Multi-Segment Rupture Model of the 2016 Kumamoto Earthquake Revealed by InSAR and GPS Data" Remote Sensing 12, no. 22: 3721. https://doi.org/10.3390/rs12223721

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop