High-resolution image on terminus of fault rupture: relationship with volcanic hydrothermal structure

SUMMARY

Earthquake-volcano interactions have been discussed to understand the underlying mechanisms of seismic ruptures or eruptions, yet the involvement of volcanic activity and the environment with fault slip termination remains unclear. Here, we present an unprecedented high-resolution image of fault motions and crustal structure at the rupture terminus in volcanic area from the 2016 Kumamoto earthquake by conducting synthetic aperture radar (SAR) data analysis and gravity inversion. We obtained a 3-D displacement field by applying multiple SAR analysis methods: standard SAR interferometry, split-bandwidth interferometry and pixel offset. We successfully mapped the ground displacements with a high-spatial resolution in the Aso caldera which was located on the eastern extension of the Futagawa fault that was the main source fault of this seismic event. We found that the rupture propagating on the Futagawa fault eastward penetrated into the Aso caldera and was divided into two major fault systems: a right-lateral fault system on the northern side and a left-lateral fault system on the southern side. However, they progressively converged immediately after penetrating into the caldera. A gravity-inferred 3-D density contrast structure revealed that a locally distributed low-density body existed in the shallow part (from the subsurface to a depth of ∼3 km) of the western edge of the caldera. The slip distribution model showed that the slips on the bifurcated faults penetrated into the low-density region and subsequently dissipated. A numerical simulation on 3-D dynamic rupture demonstrated that the low-stress state in the caldera played a role in suppressing the rupture evolution. A thermally activated hydrothermal field has developed in the area where the fault slips were attenuated. We interpret that the hydrothermal system may create conditions favourable for low-stress field, and plastic properties in the hydrothermal environment may facilitate a further decrease in rock brittleness owing to the high temperature, resulting in the terminus of fault rupture.

1 INTRODUCTION

Why an earthquake terminates remains a fundamental question in seismology. This has long been a contentious topic, during which some ideas have been proposed (King & Nabelek 1985; Kame & Yamashita 1999). A physical barrier which hampers rupture propagation is a potential underlying mechanism; a ‘hard’ barrier such as a seamount stops rupture propagation by strong frictional resistance (Kodaira et al. 2000; Yang, Liu & Lin 2013), and a ‘soft’ barrier which is a region that undergoes plastic deformation or a weak coupling area on fault, arrests rupture propagation by reducing the potential for further rupture (King 1986; Mochizuki et al. 2008; Kubo & Nishikawa 2020).

Earthquakes are often involved with volcanic activity and/or environment intimately, known as earthquake-volcano interactions. These interactions have been studied well using geophysical and geological approaches (Manga & Brodsky 2006; Eggert & Walter 2009; Gómez-Vasconcelos et al. 2020). In particular, phenomena on triggering have been numerously reported. The case of volcanic activities triggering earthquakes is one such example (Nostro et al. 1998; Nishimura et al. 2001; Toda, Stein & Sagiya 2002). Strain changes associated with a magma reservoir expansion and/or magma intrusion have been well discussed as a mechanism causing earthquakes based on observational studies and numerical calculations (Pulvirenti, Aloisi & Jin 2017; Kundu et al. 2020). Many more studies have been conducted on earthquakes that induce eruptions, and a wide range of triggering mechanisms has been discussed and proposed for many decades: static dilatation strain, static compressional strain and dynamic strain (Nakamura 1977; Rikitake & Sato 1989; Linde & Sacks 1998; Hill, Pollitz & Newhall 2002; Walter & Amelung 2007; Cannata et al. 2010; Nishimura 2021; Seropian et al. 2021). In contrast, studies that have focused on the interaction of suppressing either seismic or volcanic activity are few in number. In particular, the effects of volcanoes on the termination of fault rupture have been poorly documented.

Large inland earthquakes occur occasionally near active volcanoes. Volcanic regions may act as a ‘soft’ barrier for the rupture process owing to high temperatures, which should facilitate not brittle rupture but plastic deformation. However, it is unclear whether a volcano is involved in fault rupture propagation and, if so, how a volcano can affect fault ruptures. This is because opportunities to study the earthquake-volcano interactions on the rupture terminus are few. There is little direct evidence demonstrating an intimate relationship between the terminus of the fault rupture and the volcanic environment.

The 2016 Kumamoto earthquake was an inland earthquake with a moment magnitude (M_w) of 7.0 in western Japan (Kato, Nakamura & Hiyama 2016). It serves as an important case study of fault rupture termination. The extensive N–S graben running through central Kyushu Island (Loveless & Meade 2010), developed in the western extension of the Japan Median Tectonic Line (Fig. 1), includes active volcanoes and active faults. A few of these were involved in the 2016 Kumamoto earthquake. The earthquake began at the junction of the Futagawa and Hinagu faults, and the ruptures propagated eastward along the Futagawa fault (Kubo et al. 2016). The Aso caldera which holds one of the most active volcanoes in Japan exists at the eastern end of the Futagawa fault. The rupture stopped in the caldera without further propagating eastward; thus, some researchers have indicated that the Aso volcano is related to the terminus of the earthquake (Miyakawa et al. 2016; Yagi et al. 2016; Yue et al. 2017; Yamada et al. 2019). However, the spatial details of the terminus of fault ruptures within the Aso caldera have not yet been investigated well; hence, their relationship with volcanoes remains unclear.

Synthetic aperture radar interferometry (InSAR) is capable of creating detailed spatial maps of crustal deformation which can directly provide evidence of the extent to which a fault slip propagates. In particular, InSAR with an L-band SAR satellite, called Advanced Land Observing Satellite 2 (ALOS-2), is more suitable than C- or X-band microwaves for measuring ground displacements even in heavy vegetation area (Wei & Sandwell 2010) and/or with large displacement gradients (Hanssen 2001). However, for the 2016 Kumamoto event, even L-band InSAR could not fully map the deformation in the proximity of faults. This prevented us from obtaining the details of fault motions in the Aso caldera (Fujiwara et al. 2016). Although some previous studies have presented InSAR data even near faults, the spatial resolution is substantially degraded by applying smoothing filter and/or multilook processing, which prevented from mapping local fault motions in the proximity of faults and in the Aso caldera. Insufficient near-field data makes us difficult to obtain fault geometries involved with rupture terminus, and it is problematic that fault slips in the Aso caldera have not been investigated effectively yet. Hence, in this study, we apply not only standard InSAR but also sophisticated SAR techniques which enable to measure displacements in the caldera, to overcome the drawback, and show new findings on the terminus of fault motions. According to previous studies on crustal structure (Miyakawa et al. 2016), it has been pointed out that the presence of anomalous density structure is involved with the terminus of fault rupture in the Aso caldera. However, both the fault planes and the crustal structure discussed are large scale comparable to the caldera. For a more elaborate discussion, we also aim to obtain spatially detailed density structure by inversion with more gravity data. Against the background, the primary purpose of this study is to highlight the factors controlling the terminus of the fault rupture by identifying with high resolution where and how the fault rupture halts, and by exploring physical properties of the crust at and around the rupture terminus.

2 CRUSTAL DEFORMATION

2.1. SAR data and processing

We applied the InSAR method to ALOS-2 satellite data, acquired on 2016 April 15 and 2016 April 29, to map the spatially detailed ground deformation (Fig. 1b and Table S1, Supporting Information). We obtained InSAR images with two different viewing angles from ascending/left-looking and descending/left-looking orbit data. The difference in acquisition time is approximately 12 h, meaning that data combination is best for the ALOS-2 satellite to obtain a 3-D deformation map because effect of post-seismic deformation between the two orbit data can be negligible. The ALOS-2 data were processed using a software developed by the Geospatial Information Authority of Japan (GSI) (Fujiwara & Tobita 1999; Fujiwara et al. 1999; Tobita et al. 1999; Tobita 2003). To remove topographic fringes, we used a 10-m resolution digital elevation model (DEM) provided by the GSI. We conducted multilooking processings with 16 looks for both range and azimuth components, and then filtering operations using the modified Goldstein filter with a coefficient of 1.0 (Baran et al. 2003) for unwrapping processing. The grid spacing of geocoded images is approximately 50 m. We reduced troposphere-related noises using a numerical weather model; the Japan Meteorological Agency (JMA) mesoscale model (Kobayashi et al. 2014; Kobayashi 2016). In general, InSAR results produced using L-band SAR data often suffer from ionospheric noises. To reduce it, we applied a range split-spectrum method (Gomba et al. 2016).

In addition to the standard InSAR, we applied a split-bandwidth interferometry (SBI) method to the range and azimuth components (Jiang et al. 2017). Signals in an SAR image have an extent of bandwidth with a certain central frequency. For standard InSAR, the full-bandwidth image is used for interferometry, whereas SBI uses subband images. First, two images with lower and higher frequencies were produced by splitting the band from a full-bandwidth image, and then two interferograms were created from the subband images. Pixel shift can be measured by performing a phase-differential operation between the two interferograms. Since an SAR image has each bandwidth for the range and azimuth directions, the SBI is applicable in both directions. The range and azimuth SBI can measure displacement in the line-of-sight (LOS) and along-track directions, respectively. The application to the azimuth (AZ) component is often called multiple-aperture InSAR (MAI) (Bechor & Zebker 2006). To distinguish MAI, we here call the range-specified SBI a Range SBI (RSBI) method.

While the measurement accuracy of SBI is worse than that of the standard InSAR, the SBI can obtain better coherence even in areas with large displacement gradients. This implies that it can measure locally distributed large displacements of a few meters near faults. Owing to SBI techniques, we successfully mapped spatially detailed displacements; (1) not only in the far-field but also in the near-field, where the standard InSAR cannot effectively measure displacements, and (2) not only in the range component but also in the azimuth component which is sensitive to ground movement along the north–south orientation.

Furthermore, we applied a pixel offset method to incoherent areas where InSAR, including SBI, can hardly measure the displacement (Kobayashi et al. 2009). To obtain full displacement maps for the range and azimuth components, we performed the following steps. For the range, we appended the SBI data to the masked (incoherent) area of the InSAR, and then further appended the pixel offset data to the masked (incoherent) area of the combined image of InSAR and SBI. For the azimuth, we appended the pixel offset data to the masked (incoherent) area of the MAI. By combining these data, we obtained the full displacement field in both the range and the azimuth components.

Using the four independent components, we estimated the 3-D displacement field using the least-squares method (LSM). For a given pixel, we have an observed displacement data vector D = [d₁^LOS, d₂^LOS, d₃^AZ, d₄^AZ]^T, which comprises two LOS and two AZ components. The East–West (E–W), North–South (N–S) and Up–Down (U–D) components of displacement vector Û = [u_ew, u_ns, u_ud]^T can be solved with the following equation using LSM (Kobayashi, Morishita & Yarai 2017):

Here, G is a matrix that consists of direction cosines of the E–W, N–S and U–D components for the LOS and AZ components. For the LOS and AZ components, the direction cosines are expressed as

where θ and φ are the incidence angle of the microwave and the heading angle of the satellite flight direction measured clockwise from the north at each pixel, respectively. For the direction cosines of LOS (p_ew and p_ns), an upper sign should be used for the right-looking direction, whereas a lower sign should be used for the left-looking direction. C is the covariance matrix which is a diagonal matrix (Supporting Information). The covariance matrix consists of the errors of observation data. We calculated the errors based on coherence (Supporting Information), and then assigned the estimated errors to the covariance matrix.

2.2. SAR-unveiled terminus of fault slips in the Aso caldera

Here, we show the interferograms mapping deformation associated with the 2016 Kumamoto earthquake (Fig. 2 and Fig. S1, Supporting Information). Note that no filter is applied to InSAR data shown in Fig. 2(a) and (d). We can identify that the crustal deformation is distributed along the Futagawa fault; however, low-coherence areas spread in the proximity of the fault and in the eastern extension within the caldera (Fig. 2c and f), which prevents from measuring the displacements (Fig. 2b and e). Thus, in addition to the standard InSAR, we applied SBI and pixel offset methods, and have successfully measured the displacements even in the proximity of fault and in the caldera (Fig. 3 and Fig. S2, Supporting Information). Clear displacement discontinuities across which the orientation of displacement is opposite appear just along the geomorphologically recognized active fault traces of the Futagawa and Hinagu faults. Although the deformation is mainly produced by the two active faults, it is further noted that the displacement discontinuities can be also identified in the western edge of the caldera. The deformation pattern observed in the Aso caldera indicates partition of rupture into multiple faults, where no active fault is known to exist. This implies that the fault ruptures propagated eastward into the caldera without stopping.

To clarify the fault motions, we estimated the 3-D displacements calculated using the InSAR-, SBI- and pixel offset-derived displacements with different view angles (Fig. S3, Supporting Information). For the estimate, we calculated the errors of observation data (Fig. S4, Supporting Information), and assigned them (Fig. S5, Supporting Information) to the covariance matrix in the eq. (1). The estimated displacements are shown in Fig. 4 (Fig. S6 for the wide area and Fig. S7 for the standard deviation in the Supporting Information). The 3-D displacement fields clearly show the fault motions along the active faults and in the Aso caldera. The displacements are concentrated along the two active faults: the horizontal displacement reaches up to ∼2 m along the Futagawa fault, and the subsidence is distributed on the northern side of the Futagawa fault with the amount of ∼2 m. The fault motion is of right-lateral slip for both the two active faults, while normal fault motion is substantially included on the Futagawa fault, which are consistent with the previous studies (Kubo et al. 2016; Yoshida et al. 2017; He et al. 2020). Two major displacement boundaries across which orientation of displacement is opposite and/or displacement has a steep gradient, which correspond to top edge of faults, were systematically identified in the Aso caldera, hereafter referred to as F1 and F5 (Fig. 4). The strike of F1 was almost parallel to the Futagawa fault and had a right-lateral motion with the amount of ∼50 cm. F5 branched from the right-lateral fault system; it first ran in an east–west orientation, hereafter called F4, and then in northwest–southeast orientation. The displacement reaches up to ∼50 cm in horizontal component and ∼1 m in vertical component. The horizontal motion is of a left-lateral motion, as identified from displacement profiles and field survey (Fig. 5), although F1 is of right-lateral motion. We measured the displacement at the point Q in the field survey (Fig. 5f). The SAR-derived displacements at the corresponding point are ∼11 cm with left-lateral motion and ∼15 cm with subsidence in the northern side, which are almost consistent with the results from the field survey. The subsidence area is located on the southern and northern sides of F1 and F5 (F4), and the ground sandwiched between F1 and F5 subsides like a graben. If normal fault motions are present, F1 and F5 (and F4) should be south- and north-dipping fault planes, respectively.

In addition to F1, F4 and F5, we found locally distributed linear displacement boundaries called F2 and F3. Their strike angles are almost parallel to F1. Clear discontinuities were identified in the horizontal components. These displayed right-lateral motions similar to those of as F1. The horizontal displacement reaches up to ∼50 cm. The ground slightly subsides with the amount of ∼10 cm on the southern side of F2, whereas for F3, little vertical displacement can be recognized.

The SAR-derived displacement map reveals not only the bifurcation of fault rupture but also the termini of the fault ruptures. All the displacement boundaries F1–F5 exibit clear displacement offsets corresponding to shallow fault slips. However, the displacement boundaries are locally distributed with a maximum length of a few kilometers, suggesting that the slips terminate soon after penetrating the caldera. For F1, the displacement offset can be clearly observed at the western edge of the F1 boundary, indicated by a solid white line (Fig. 4e). However, the displacement boundary becomes blurred, as indicated by the broken line. The clear displacement boundaries for F2 and F3 have the length of ∼2 and ∼1 km, respectively; however, at the eastern edge, the displacement offsets sharply vanish. The boundary of F4 runs in E–W orientation with a length of ∼2 km, but the slip does not proceed further, and the strike changes to a northwest–southeast orientation, connecting to F5. The displacement boundary of F5 can be identified with a length of a few kilometers, but it gradually becomes unclear, as indicated by the broken white line (Fig. 4a).

One may think that the SAR-derived deformation in the caldera is contaminated by post-seismic deformation because it is derived by SAR images acquired about two weeks after the earthquake (Table S1, Supporting Information). However, we think that it almost consists of coseismic deformation although there is no direct evidence that the deformation was created by dynamic rupture. First, we can say that no significant aseismic slip occurred since 23 hr after the earthquake. This is because any fault motions cannot be identified in an InSAR image just after the earthquake (Fig. S8, Supporting Information), which is created by images acquired on 2016 April 17 (00:04 in Japan Standard Time, JST) and 2016 May 1. Second, according to the seismic wave analyses, seismic energy is substantially released also in the caldera (Asano & Iwata 2016; Kubo et al. 2016; Yoshida et al. 2017): shallow slips with ∼2 m, which is consistent with the slip amounts inferred from the SAR data (shown in chapter 4), are estimated in the west of the caldera. Further, the Global Navigation Satellite System (GNSS)-observed post-seismic deformation at 960701 site (Fig. 4), which is closest to the SAR-derived deformation area, is estimated to be −1.2, −1.5 and −1.2 cm in E–W, N–S and U–D components, respectively, during 24 hr just after the earthquake. During the two weeks after the earthquake, the deformation at 960701 is −1.4, −1.5 and 1.3 cm in E–W, N–S and U–D components, respectively. The displacements are much smaller than the displacements (−12.8, −54.1 and 14.9 cm) predicted by the fault models in the caldera (Flt1–Flt5 proposed in chapter 4). Thus, the SAR-derived deformation was unlikely produced after the earthquake, and it would be reasonable to assume that the deformation observed in the caldera is not post-seismic but coseismic signals: it reflects the eastward extension of dynamic rupture on the Futagawa fault.

3 GRAVITY-INFERRED DENSITY CONTRAST STRUCTURE

3.1 Gravity data

Upon considering the spatial details of the termination of fault slips, we inferred a 3-D density contrast structure using gravity data observed in and around the Aso caldera to examine whether there were any specific features in the crustal structure at and around the end of the fault ruptures. First, we will explain the gravity data used. We used 5185 land gravity data points for the inversion of the 3-D density contrast structure of the Aso caldera, provided by the Geological Survey of Japan (Miyakawa et al. 2015), Gravity Research Group in Southwest Japan (Yamamoto, Shichi & Kubo 2011), Kanazawa University (Honda et al. 2012), and GSI (Supporting Information). The positions of the gravity observation sites are shown in Fig. 1(b). Our data used for the inversion are substantially more voluminous than those used in previous study (Miyakawa et al. 2016): in the latter case, the used gravity sites were 512, meaning that there is a 10-fold difference between the two. In particular, in the west of the caldera, higher dense observations are achieved by adding the data of the GSI. These gravity data were obtained using relative gravimeters and referenced to the Japan Gravity Standardization Net 1975 (Geographical Survey Institute 1976).

We derived the complete Bouguer gravity anomalies at the observation points from the gravity data using the following equation (Kuhn et al. 2009):

where g is the observed gravity, ${\gamma }_0$ is the normal gravity on the surface of the reference ellipsoid, $\delta g_{\mathrm{FC}}$ is the free-air correction, $\delta g_{\mathrm{AC}}$ and $\delta g_{\mathrm{TOP}}$ are the gravitational attraction of the atmospheric masses (atmospheric correction) and topographic masses (topographic reduction), respectively (see the details in Supporting Information). The density of topographic masses in the Aso Caldera was estimated to be 2200–2300 kg m⁻³ by Komazawa (1995). Here, we used a value of 2300 kg m⁻³ for the computation of the topographic reduction. To reduce the computation time while maintaining accuracy, we computed the topographic reduction by splitting the integration area into three zones: a central-zone (0–50 km), a near-zone (50–300 km), and a far-zone (the rest of the Earth). The contributions of the central-, near- and far-zones were evaluated using 0.4 $\times$ 0.4 arcsec grid DEM provided by GSI (http://www.gsi.go.jp/kiban/faq.html), 1.0 $\times$ 1.0 arcmin and 1.0 $\times$ 1.0 arc-deg grid DEM from Shuttle Radar Topography Mission V4.1(Jarvis et al. 2008), respectively.

3.2 Three-dimensional gravity inversion for density contrast structure

We estimated the 3-D density contrast structure of the Aso caldera from the distribution of the complete Bouguer gravity anomalies using a gravity inversion software package GROWTH 3.0 (Camacho et al. 2021). The software determines the density contrast structure of an aggregate of small rectangular prisms using an LSM with a mixed minimization condition of model ‘fitness’ (least-square minimization of residuals) and model ‘smoothness’ (l₂-minimization of total anomalous mass). The software requires an input of the model parameters of a random search coefficient, a balance factor, a flattening coefficient, a lateral smoothing value and an upward weighting value. We used the default values offered in the software. The linear regional trend of the gravity anomaly field and the general downward density increase were jointly estimated using the 3-D density contrast structure. The area of the inversion ranged from 32.5°–33.2°, from 130.7°–131.4°, and from a depth of 10 000 m below sea level to an altitude of 1500 m above sea level. The mean length (cell size) of the 3-D grid was set to be 750 m.

3.3 Anomalous density structure in the Aso caldera

Fig. 6(a) presents a gravity-estimated density contrast structure at a depth of 2 km, where the estimated values represent the difference from the averaged density at the same depth. Low-density areas were significantly distributed in the central part and in the northern part of the Aso caldera. The density is up to ∼300 kg m⁻³ lower than the surrounding density. In the central part, the low density spreads around the Nakadake crater where eruptions have repeatedly occurred. Of note, a low-density body exists at the western edge of the caldera from the subsurface to a depth of approximately 3 km (Fig. 6c and Fig. S9, Supporting Information), where the fault slips terminated. The low density continuously extends from the crater at a depth of 1– 3 km.

Fig. 6(b) presents the spatial relationship between the low-density area and SAR-detected displacement boundaries. The bifurcated faults overlap with the low-density body. F1 is located along the northern edge of the low-density area. F2 and F3 are inside the low-density aera and extend towards the centre. F4 is also inside the low-density, but changes the direction southeastward (F5) on the half way, and finally disappears after contact with the low-density area. It suggests that the fault slips penetrated the low-density body and stop their propagation. The eastern edges of the displacement boundaries (white solid lines) reached around −300 kg m⁻³. The good spatial correlation implies that the low-density body is intimately related to the propagation of fault slips.

4 FAULT MODEL

4.1. Processing for the slip distribution model

To scrutinize the spatial relationship between the fault slips and the low-density body, we constructed slip distribution models for faults corresponding to F1–F5, hereafter referred to as Flt1–Flt5, by inverting the SAR-derived crustal deformation data. There are some displacement discontinuities, and the spatial pattern of crustal deformation is rather complicated in the source region, meaning that multiple faults are intimately involved with the Kumamoto earthquake. To unveil the detailed fault system and spatial pattern of fault slip in the area of interest of the Aso caldera, we constructed a fault model with the following two steps. First, we roughly determined the fault geometries with parameter searches with a uniform slip on a rectangular plane. Here, position, size and orientation of fault, and amount and orientation of slip are estimated. Then, we estimated the slip distribution on the fault plane determined in the previous step by LSM; however, dip angles of faults in the caldera are numerically searched to get more optimal geometry.

4.1.1. Fault model with uniform slip on rectangular plane

We constructed a fault model by assuming the presence of a rectangular fault with uniform slip in an elastic half-space (Okada 1985). The model parameters were estimated using a simulated annealing method (Metropolis et al. 1953; Cervelli et al. 2001). We randomly assigned parameters for search and adopted parameters as the optimal solutions, which gave the minimum value of the weighted sum of squared residuals between the observation d^obs and the calculation d^cal.

where N, M and σ represent the number of data samples, number of model parameters and standard deviations of the observation data, respectively. We assigned the estimated errors (Fig. S5, Supporting Information) as σ. We employed the bootstrap method to estimate the individual confidence of the inferred parameters (Efron 1979).

To reduce the number of SAR data for the modelling analysis, we resampled the SAR-derived data beforehand, using a quadtree decomposition method (Jónsson et al. 2002). For a given quadrant, if, after removing the mean, the residual is greater than a prescribed threshold, then the quadrant is further divided into four new quadrants. We set the threshold to 2 and 5 cm for the LOS and the AZ component data, respectively. This process was iterated until either each block met the specified criterion or until the quadrant reached a minimum block size of 8 × 8 pixels. Furthermore, we also utilized GNSS data for our modelling as auxiliary data (Fig. 1b, Supporting Information).

We set nine fault planes for modelling by following the displacement boundaries identified in the SAR displacement maps. Five faults were set up in the Aso caldera: Flt1, Flt2, Flt3, Flt4 and Flt5 (Fig. 1b). Based on previous studies (He et al. 2020) and a reasonable fault motion as pointed out in the main text, we assumed that Flt1 and Flt2 were dipping to the south. For Flt3, little discontinuity was observed in the vertical component in SAR displacement field; thus, we assumed that the fault plane of Flt3 was vertical to prevent from an increase in the number of model parameters solved. In contrast, Flt4 and Flt5 were set to be north-dipping planes. This is because it is thought to be unnatural that reverse fault motions are required to account for the observed subsidence in the case of the south-dipping plane.

In addition to these faults in the Aso caldera, we set four other faults (Flt6 to Flt9) along the Futagawa, Hinagu and Denokuchi faults which were mainly responsible for the crustal deformation of the 2016 Kumamoto earthquake (Fig. 1). Here, we assumed that these faults were north- (northeast-) dipping, based on previous researches (He et al. 2020). For the parameter search, we first fixed the position and strike with knowledge of the displacement boundaries so as to reduce the number of estimates, and then estimated depth, length of fault, width of fault, dip angle, rake angle and slip amount by the simulated annealing method. Note that depth and length are fixed for some faults that show clear displacement boundaries. The estimated fault parameters and their estimate errors are listed in Table 1.

4.1.2. Slip distribution model

Next, we constructed a slip distribution model with knowledge of the inferred fault geometries in the previous step. To reduce uncertainties for the modelling with the nine fault planes, we basically fixed some fault parameters which were reliable from the observation, except for dip angle, to estimate the slip distribution. Table 2 lists the preset locations (longitude, latitude and depth) and geometries (length, width and strike angle). There is little seismicity in the Aso caldera (Shito et al. 2017; Mitsuoka et al. 2020); thus, it is infeasible to determine the dip angles using aftershocks distribution. Thus, we further searched for the dip angles for faults in the Aso caldera. The corresponding faults Flt1, Flt2, Flt4 and Flt5 were searched by a grid search within the ranges from 20°–80°, from 60°–90° from 30°–70°, and from 40°–80°, respectively, with every 10°.

The fault was divided into square patches with a size of 1 × 1 km for the faults in the Aso Caldera (Flt1, Flt2, Flt3, Flt4 and Flt5); while a size of 2 × 2 km was assigned to faults outside the Aso caldera (Flt6, Flt7, Flt8 and Flt9). The fault planes for Flt3 and Flt4 consist of one and two in the strike orientation, respectively (Fig. 7 and Table 2), and their number may not be sufficient to represent the slips. However, even if increasing the number of patches, we could not obtain more significant results statistically (Fig. S10, Supplementary Information), thus we took the 1 × 1 km patch size in this study. We used the dislocation equations derived by Okada (1985) to calculate the surface displacement. Although there are more elaborated tools that can handle more realistic system with non-planar faults (Meade 2007; Nikkhoo & Walter 2015), observational information on fault geometry at deep part whether or not the faults are intersecting and/or bending is missing. In this study, we adopted Okada's formulation (planar fault), which is more suited to the approach of numerically searching for the fault dipping, although it is classical, and aimed to extract the features of the fault system.

To suppress the instability of the solutions due to the increase of model parameters, we imposed a spatial smoothness constraint on the slip distribution using a Laplacian operator. Consequently, we solved the following equation using LSM.

where d, G and m represent the displacement data vector, corresponding Green's matrix and fault slip vector comprising the dip-slip and the strike-slip motions, respectively. L and α represent a Laplacian operator matrix and hyperparameter which controls the degree of smoothness, respectively. The relative weights α were determined using Akaike's Bayesian information criterion (ABIC) (Akaike 1980). We incorporated the standard deviations of observation data (Fig. S5, Supporting Information) into the least-squares approach as weight. The contribution of the cross-terms of the covariance matrix for the InSAR data cannot be neglected, in as much as they have a spatial correlation in general, which largely results from the variations of atmospheric water vapour (Lohman & Simons 2005; Fukahata & Wright 2008). We incorporated the cross-terms of the covariance matrix in the inversion scheme, following the equation presented by Lohman & Simons (2005). We here took the characteristic correlation distance of errors to be 20 km. Through the LSM, we estimated the standard deviations for the slip amounts.

When there is no constraint on the depth of the lower edge of fault plane, the slip distribution sometimes becomes unnatural and/or no optimal ABIC is obtained. Thus, we imposed a constraint on the lower limit of fault depth. Based on the knowledge of the seismogenic zone in and around the source region (Matsumoto, Nishimura & Ohkura 2016; Shito et al. 2017), we assigned a lower limit of ∼5 km in depth for the faults in the Aso caldera (Flt1 to Flt5), and for the outside of the Aso caldera (Flt6 to Flt9), the lower limit was assumed to be ∼15 km in depth.

4.2. Features on spatial distribution of fault slip

We constructed slip distribution models (Fig. 7 and Fig. S11, Supporting Information). The slips on Flt1–Flt3, with right-lateral motions including slight normal motions, are concentrated at a depth of a few kilometers for all faults (Fig. 7). The depth is consistent with the depth of the low-density body inferred from the gravity inversion. The maximum slips on Flt1 and Flt2 are located at the western edge on each fault plane with ∼2 and ∼1 m in slip amount, respectively, which significantly exceeds the error level (Fig. S12, Supporting Information), and the slips gradually dissipate with eastward propagation. The fault system comprising Flt4 and Flt5 have left-lateral slips with normal motion which becomes dominant on Flt4 and the western edge of Flt5. The slips are concentrated at depths of 2–3 km samely as Flt1–Flt3, and the slip on Flt5 gradually dissipates as it propagates eastward. Although there are large slips at some patches in the fault edges, which may raise a concern on physically unnatural displacement gradients, it seems that the shear stress does not necessarily exceed the shear strength of the surrounding rocks (Supporting Information). The total released M_w is estimated as 7.05 for all the faults, while the seismic moment released by the faults inside the caldera (Flt1–Flt5) reaches up to 9.9 per cent of the total seismic moment (Table 2).

In the inversion, we numerically searched for optimal dip angles for Flt1 to Flt5. Fig. S13 (Supporting Information) shows the results of the numerical search. Herein, the residual distribution is plotted as a function of the dip angle of Flt1. The vertical axis indicates the residuals. The horizontal axis indicates the dip angle of Flt1. In each panel, the dip angle of Flt2 is fixed at a certain angle, as indicated in the upper left. Individual lines on each panel are the residuals of when the dip angles of Flt4 and Flt5 are fixed at certain angles. The line types represent the dip angles of Flt4, whereas the colours represent the dip angles of Flt5. Dip angles that yield the minimum value of the residual are adopted as the optimal solutions. The final solutions for the dip angles are listed in Table 2. The numerically searched fault planes of Flt1, Flt4 and Flt5 have low- and moderated-dip angles, whereas those of Flt2 and Flt3 have high-dip angles (Tables 1 and 2), suggesting that the fault system in the caldera forms a kind of graben structure.

5 DISCUSSION

5.1 Numerical simulation on three-dimensional dynamic rupture: effect of low-stress field

The observed distribution of the surface faults is similar to that of the so-called ‘horse-tail structure/splay faults’. This type of fracture style is often observed at the ends of faults. If the elastic energy release of the main fault is distributed onto individual bifurcated faults, the driving force for rupturing each fault is reduced, suppressing further rupture growth. Here, we tested through a numerical simulation of 3-D dynamic rupture propagation whether the slip termination was attributed to the distributed energy release onto such multiple faults. We briefly examined the rupture evolution for Flt1–Flt5 for the two cases of modelled stress states: one is a homogeneous regional stress field over the area both in and off the caldera (Matsumoto et al. 2015), and the other is an inhomogeneous stress field where the low-stress state corresponding to the low-density body is incorporated into the medium in the caldera.

The dynamic fault ruptures were numerically simulated following the algorithm and modelling concept presented in Ando et al. (2017). We applied the fast domain partitioning boundary integral equation method developed by Ando (2016) for the fully dynamic simulation. Here, we assumed the frictional strength as the linear slip weakening friction law (Ida 1972), where the static and dynamic frictional coefficients and critical slip-weakening distance were set to be 0.3, 0.2 and 0.4 m, respectively, by following previous models (Ando et al. 2017).

The objective of the numerical test is to roughly estimate the sensitivity of the rupture evolutions for the Flt1, Flt2, Flt3 and Flt4–Flt5 faults with respect to the spatial change of the stress state. Fault geometry and stress state are parameters that can be constrained by observations to characterize individual faults and are primary control of the rupture process and the resulting slip heterogeneity (Ando et al. 2017). For the simulation, we set the fault geometries as shown in Fig. 8, where we connected Flt5 with Flt6 continuously at depth, considering both faults dipping north. The fault planes were discretized by triangular meshes using Gmsh software (Geuzaine & Remacle 2009). The rupture started on Flt6, corresponding to the Futagawa fault.

Here, we assumed a stress state constrained by the principal stress orientations and stress ratio, based on the stress inversion analysis by Matsumoto et al. (2015). The observed stress regime is characterised by the N–S extension and the E–W compression as shown in Fig. 8a. A ratio between the principal stresses (σ₂–σ₃)/(σ₁–σ₃) was constrained to be 1 leading to the trans-tensional stress regime, where σ₁, σ₂, and σ₃ represent the maximum, intermediate and minimum principal stresses, respectively. Here, σ₂ is a vertical principal stress. The vertical principal stress (σ₂) was determined by the hydrostatic overburden pressure as σ₂ (h) = (ρ_s − ρ_f)gh ∼ 16.7 (MPa km⁻¹) from the average densities of the crustal rocks ρ_s and fluid ρ_f, the gravitational acceleration g, and the depth h. The stress ratio σ₃/σ₁ was given so that the average stress drop at the epicentre was approximately 10 MPa.

In this study, we simply simulated two end-member cases: one is the case of the uniform stress distribution over the entire area without the stress heterogeneity in the Aso caldera, and the other is the case with a low-stress area in the caldera where the low density is observed. We numerically tested these two cases because the available observation has insufficient spatial resolution of the stress heterogeneity in the caldera scale. In the latter case, we modified the initial shear stress at a low level to prevent dropping the stress on the faults in the caldera (Fig. S14, Supporting Information). The low-stress state was implemented by making the initial shear stress the same as the dynamic frictional strength. Furthermore, as in the source modelling, with the knowledge of the lower limit of the seismogenic zone in the caldera (∼5 km in depth), the modelled fault area extends down to 5 km to inhibit the rupture propagation beyond the depth. We limited our study to qualitative explanations. More detailed sensitivity tests to frictional and fault geometry parameters are beyond the scope of this study.

For Flt1, the final slip distribution of the numerical simulation (Fig. 8 and Fig. S15, Supporting Information) produces a similar spatial pattern to that of the geodetic data inversion (Fig. 7): the slip is concentrated at the western edge of the fault plane in the shallow part, and then weakens and becomes stronger again as it moves eastward. There was no substantial difference in slip distribution between the two modelled stress states. The fault orientation of Flt1 with a low-dip angle 30^o is not preferred by the observed trans-tensional (strike-slip and normal faulting) stress regime in this area. Flt1 was poorly suited to hosting a slip even before imposing the low-stress condition. However, we could recognize the effect of low stress on Flt5 with a higher dip angle of 50^o. In the homogeneous stress field, the fault rupture on Flt5 easily propagates to the eastern end of the fault plane without stopping, while in the case of low stress, the fault slip gradually terminates at the deep part by propagating eastward without reaching the eastern end (Fig. S16, Supporting Information). The final spatial patterns of the slip distribution (Fig. 8f and Fig. S15f, Supporting Information) are consistent with those of the slip distribution model (Fig. 7). This result is robust and insensitive to the chosen values of frictional coefficients because the obliqueness of Flt1 and preference of Flt5 were not significantly altered by the variation in the internal frictional angle of ∼15° for the significant difference in the frictional coefficient of 0.3. Also for Flt2, low-stress state plays a role to produce a similar spatial pattern to that of the slip distribution model determined by the inversion. The slip simulated in the homogeneous stress field is almost uniformly distributed over the entire plane (Fig. 8 and Fig. S15c, Supporting Information), which is inconsistent with the slip distribution model (Fig. 7). In contrast, in the low-stress field, no significant slip is distributed at the depth, instead, the slip concentrates at a shallow part only, which is similar to the slip distribution model (Fig. 8 and Fig. S15d, Supporting Information). These results suggest that the idea of distributed energy release under the assumption of homogeneous stress field cannot fully explain why ruptures in the Aso caldera terminate within a length of a few kilometers.

Locally developed stress that is different from the regional stress orientation, which may be produced by a volcanic activity such as a volume change of a magma chamber, could be a factor that controls fault slips; however, this effect would be small. For the Aso volcano, it has been inferred that a magma chamber is possibly located at the caldera centre, near the Nakadake crater (Fig. 1), at approximately 5 km depth (Sudo et al. 2006; Nobile et al. 2017). This is a candidate source which could produce a local stress change within the caldera. The magma chamber is inferred to be located due east (131.07° E, 32.89° N) of the fault system of Flt1–Flt5. If the magma chamber deforms at the inferred position, the maximum principal axis of stress produced by the deformation should be oriented approximately east–west. It deviates only ∼10° from the orientation of the regional stress assigned in the simulation. In this case, the stress perturbation has little effect on the slip distribution (Fig. S17a–d, Supporting Information), and the simulation results are agreement with the observations.

The validity of the stress field assigned in the numerical simulation can be also confirmed from the simulated slip directions. Here, it should be noted that the numerical simulations have successfully reproduced not only the spatial patterns on slip amounts but also the slip orientations: our simulation result produced the left-lateral slip on Flt5 (Fig. 8g and h, and Fig. S18, Supporting Information), although the right-lateral slip is predominant for other faults as expected. Some may find it unnatural that the coexistence of right-lateral and left-lateral slips in a fault system, but it is reasonable, considering the orientation of maximum principal stress axis and the fault geometries. The principal stress axis in the Aso caldera is in a favourable orientation to promote right-lateral motion on Flt1–Flt3 and left-lateral motion on Flt5. If the stress state would largely change, the features obtained in the geodetically inferred slip distribution could not be reproduced. For example, when assuming that the maximum principal stress axis is rotated 45° clockwise, not a left-lateral but a right-lateral motion is produced on Flt5 although right-lateral slips are kept on other faults (Figs. S17g and h, Supporting Information).

The fault system in the caldera is complex, where faults are closely located one another, thus the slips could be affected by the proximity faults. In particular, slips in the shallow part can be passively driven by the surrounding fault motions because the confining pressure by lithostatic stress is weak. Flt2 for which only the shallow slips remain is the good example. The simulation results show that the slips on Flt2 tend to be suppressed by low-stress state, but the shallow slips still remain (Fig. S15d, Supporting Information). It suggests that the shallow slips were passively driven by not the deep slips on the same fault plane but the stress perturbation produced by the surrounding fault motions. By analogy with similar mechanism, the interaction between fault motions may have contributed to the shallow slips also for other faults. On the other hand, the main features of the slips on Flt1 and Flt5 are not significantly related to the presence of Flt2 and Flt3. Numerical simulations that exclude Flt2 and Flt3 do not change the main features of the slip distribution (Fig. S19, Supporting Information).

Theoretical studies suggest that a large number of subsidiary faults (Ando & Yamashita 2007) and off-fault inelasticity (Wollherr et al. 2018) may play roles in the dynamic rupture processes. Considering the low-density and low-velocity zone, the western part of the caldera may be highly fractured and there could be many unrecognized shallow faults. We here limited our model to include the least geometrical complexity visible in the observations because hidden faults are difficult to be validated with the currently available observations, while further extensive study incorporating the effects of hidden faults would be needed for better understanding on rupture behaviours in the caldera.

The high reproducibility of the slip distribution indicates that the numerical simulation obtains reasonable solutions, meaning that it would validate our interpretation that low stress corresponding to low density can effectively act as a barrier for rupture propagation.

5.2. Anomalous crustal structure in and around terminus of fault

Anomalous structures at the western edge of the caldera have also been detected in other geophysical surveys. According to the seismic wave velocity structure (Huang et al. 2018), a low S-wave velocity (Vs) is located at shallow section west of the caldera. The low Vs with an area of ∼2 km s⁻¹, indicated by the contours in Fig. 6(c), overlaps with our estimated low-density body. Furthermore, according to previous studies on electrical resistivity, a low resistivity was also inferred in the same region (Matsushima et al. 2020).

Hydrothermal system may cause such anomalous physical properties. A magma chamber could be one of candidates for low Vs and low resistivity, but the magmatic source involved with the volcanic activity of Aso volcano has been geodetically estimated to be below not the west of the caldera but the caldera centre, near the Nakadake crater (Fig. 1), at approximately 5 km in depth (Sudo et al. 2006; Nobile et al. 2017). Most of eruptions recorded in historical times occurred at the Nakadake crater. The only other eruption outside the Nakadake crater occurred in 1816 at Yunotani hot spring, which is located in our study area (Fig. 6a). It was not a magmatic but a phreatic eruption. Based on the knowledge, a hydrothermal field, to which heat is supplied from the magma chamber positioned at the central part of the caldera, may have developed in and around the anomalous structure. Fluid-filled media of hydrothermal system could be responsible for low Vs and low resistivity. Porous media filled with fluids, which is supposed to be developed in hydrothermal system, can account for the low density (Tiede et al. 2005). The anomalous structure possibly suggests hydrothermal environments where hydrothermal reservoirs and/or pathways containing volcanic fluids supplied from the magma chamber are predominant developed in the shallow part of the western area of the Aso caldera. The hydrothermal field may cause the low-stress field in the caldera which suppresses the rupture evolution as shown in the numerical simulation.

The area where the fault slips are attenuated is a thermally activated hydrothermal field. This anomalous area correlates with hydrothermal surface manifestations, such as fumaroles (Fig. 6a). The geothermal gradient is anomalously high in this area: the temperature observed at the drilled site reaches ∼200° C at a depth of 1 km (Matsushima et al. 2020). As rocks get more fluid with increase of temperature (above 300° C), the breakdown stress drop decreases, and the critical distance Dc increases (Kato et al. 2003). It may allow for fault slip to terminate as observed. To confirm it, we conducted an additional numerical simulation on dynamic rupture incorporating higher temperature condition. To incorporate the effect of high temperature, we changed the frictional parameters of linear slip weakening friction law: the value of Dc and the breakdown stress drop were set to twice and half of those outside the caldera, respectively. As a result, the fault slip for Flt5 gradually terminated at the deep part by propagating eastward without reaching the eastern end (Fig. S20, Supporting Information), similar to the result without the effect of high temperature (Fig. 8). Also for Flt1, the characteristic of slip distribution does not change substantially. Our simulation suggests that the low-density material is substantially responsible for the fault rupture termination among the three factors: complex bifurcated fault geometries, low-density material corresponding to low-stress/weak fault-strength, and enhancement of fault slip stability by temperature increase.

The thermally active state can hamper brittle rupture (Scholz 1998). The transition from brittle to plastic deformation typically occurs at ∼350° C (Scholz 2019). Hydrothermal system can also be responsible for suppression of fault rupture because hydrothermal alteration or high-pore-fluid pressure causes substantial reduction of rock strength (Wong & Baud 2012; Heidari, Khanlari & Torabi 2014; Ball et al. 2018; Heap et al. 2021). In addition to the low density that is responsible for low-stress field, the presence of locally distributed hydrothermal system can facilitate decrease of rock brittleness due to weakening of the crust, resulting in terminus of fault rupture. The few aftershocks that occurred on the western side of the caldera (Fig. 1a) may also support this idea that the crustal region has little potential for brittle fracture.

The complete picture obtained in this study is illustrated in Fig. 9. Considering Flt1 to Flt5 as one fault system, our fault model suggests that the outer faults (Flt1, Flt4 and Flt5) are of low- and moderated-dip angles, whereas the inner faults (Flt2 and Flt3) have high-dip angles (Tables 1 and 2), which is similar to the so-called ‘negative flower structure’ (Fossen 2016). A trench survey conducted before the 2016 event supports such an estimated fault-structure system (New Energy & Industrial Technology Development Organization 1995), which may indicate that the fault system has already been developed. The inward dipping faults have been often developed within a volcanic structure during caldera formation (Acocella 2007), which may be a cause of why the dipping in the caldera (Flt 1) is opposite to the Futagawa fault. The fault slips observed in the caldera might have occurred on pre-existing fault planes, similar to the Futagawa and the Hinagu faults. The elevation in the SAR-detected subsidence area sandwiched between F1 and F4–F5 is slightly lower than the surrounding area (Fig. 5), implying that the area has repeatedly subsided historically. Assuming that the topographic difference has been created by historical fault slips, the mean recurrence interval is estimated to be 3188 yr (Fig. S21, Supporting Information). It is consistent with that estimated by trench investigation near F3: 2000–3700 yr (Toda et al. 2019). The geomorphogically estimated mean slip rate and recurrence interval support the idea that the rupture on the Futagawa fault had repeatedly penetrated the caldera and stopped every time historic earthquakes had occurred.

6 CONCLUSIONS

We addressed the issue on terminus of fault rupture through cross-disciplinary research that consists of SAR-based crustal deformation analyses, gravity-based crustal structure analysis, and a numerical simulation of dynamic rupture. Our findings and interpretation are conceptually summarized in Fig. 9. The following conclusions were derived from the research;

1. We successfully mapped the high-resolution 3-D crustal deformation field within the Aso caldera for the 2016 Kumamoto earthquake, by applying standard InSAR, SBI and pixel offset methods to ALOS-2 satellite data.

2. The high-resolution deformation map unveiled that the rupture propagating on the Futagawa fault eastward penetrated into the Aso caldera without stopping, and was divided into two major fault systems: one is a right-lateral fault system on the northern side, the other is a left-lateral fault system on the southern side. The former runs parallel to the Futagawa fault although the north-dipping plane changes to the south-dipping one. Some locally distributed fault motions can be also identified. On the other hand, the latter initially runs in E–W orientation with a length of ∼2 km, but the strike changes to a northwest–southeast orientation.

3. The displacement map also revealed the details of the termini of the fault ruptures. The displacement boundaries suggesting fault motions within the caldera are locally distributed with a maximum length of a few kilometers, suggesting that the slips terminate soon after penetrating the caldera.

4. The estimated density contrast structure revealed that a low-density body existed at the western edge of the caldera from the subsurface to a depth of approximately 3 km. The bifurcated fault motions overlap with the low-density body.

5. The slip distribution model showed that the slips were mainly distributed at the western edge on each fault plane at depths of 2–3 km, and gradually dissipated with eastward propagation. The slip areas are consistent with the depth of the low-density body inferred from the gravity inversion. The spatial consistency suggests that the fault slips penetrated into the low-density body and stop their propagation.

6. Numerical simulation of 3-D dynamic rupture demonstrated that low-stress state corresponding to the low-density in the caldera played a role in suppressing the rupture evolution. The simulation reproduced well the spatial patterns of slips for each fault plane and both the right- and left-lateral fault motions estimated by geodetic inversion, suggesting that the fault ruptures and their terminations were controlled by the regional stress field and the low-stress state in and around the Aso caldera.

7. The area where the fault slips were attenuated is a thermally activated hydrothermal field. The presence of locally distributed hydrothermal system may be responsible for the low-stress. Further it might facilitate decrease of rock brittleness due to hydrothermal alteration and/or high-pore-fluid pressure. This, in turn, may result in the terminus of fault rupture. Our results demonstrate that a ‘soft’ barrier can potentially control rupture propagation.

ACKNOWLEDGMENTS

ALOS-2 data were provided from the Earthquake Working Group under a cooperative research contract with the Japan Aerospace Exploration Agency (JAXA). The ownership of ALOS-2 data belongs to JAXA. Numerical weather model data were provided from the JMA under a cooperative agreement with the JMA. Hypocentre data processed by the JMA were used. This work used computational resources of Earth Simulator provided by JAMSTEC through the HPCI System Research Project (Project ID: hp220105). We used Generic Mapping Tools (GMT) provided by Wessel & Smith (1998) to construct the figures. This study was partially supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI grant programme (Grant Number JP18K03810). Finally, we thank the editor (Dr Ana Ferreira) and the three anonymous reviewers whose constructive comments helped improve and clarify the paper.

AUTHOR CONTRIBUTIONS

TK conducted SAR data analysis, constructed fault models, and drafted the manuscript. KM inferred gravity-inferred density contrast structure. RA and GW conducted numerical simulation of 3-D dynamic rupture. TN conducted geomorphological researches. All the authors read and approved the final version of the manuscript.

CONFLICT OF INTEREST

The authors declare that they have no competing interests.

DATA AVAILABILITY

We used a 10-m resolution DEM provided by the GSI (https://fgd.gsi.go.jp/download/mapGis.php?tab=dem). ALOS-2 data were provided from the Earthquake Working Group under a cooperative research contract with the JAXA. Hypocentre data can be downloaded from https://www.data.jma.go.jp/svd/eqev/data/bulletin/hypo.html. Numerical weather model can be downloaded from http://database.rish.kyoto-u.ac.jp/arch/jmadata/. The gravity data are referred to Yamamoto et al. (2011), Miyakawa et al. (2015) and Honda et al. (2012). The gravity data provided by GSI are available on request. GNSS data can be downloaded at https://terras.gsi.go.jp/pos_main.php. The analysis data obtained in our paper are available on request. Mapping tools used in this study: The GMT, QGIS, and Paraview can be downloaded at https://www.generic-mapping-tools.org/download/, https://qgis.org/en/site/forusers/download.html and https://www.paraview.org/download/, respectively. Gmesh can be downloaded at https://gmsh.info/. A 3-D gravity inversion tool GROWTH3.0 can be downloaded at https://github.com/josefern/GROWTH_3.0_software. The other codes, excluding a part of SAR analysis code, are available upon request.

REFERENCES