Postseismic deformation following the 2016 Kumamoto earthquake detected by ALOS-2/PALSAR-2

I have been conducting a study of postseismic deformation following the 2016 Kumamoto earthquake using ALOS-2/PALSAR-2 acquired till 2018. I apply ionospheric correction to interferograms of ALOS-2/PALSAR-2. L-band SAR gives us high coherence enough to reveal surface deformation even in vegetated or mountainous area for pairs of images acquired more than 2 years. Postseismic deformation following the Kumamoto earthquake exceeds 10 cm during 2 years at some spots in and around Kumamoto city and Aso caldera. Westward motion of ~ 6 cm/year was dominant on the southeast side of the Hinagu fault, while westward shift was detected on both sides of the Futagawa fault. The area of latter deformation seems to have correlation with distribution of pyroclastic flow deposits. Significant uplift was found around the eastern Futagawa fault and on the southwestern frank of Aso caldera, whose rate reaches 4 cm/year. There are sharp changes across several coseismic surface ruptures such as Futagawa, Hinagu, and Idenokuchi faults. Rapid subsidence between Futagawa and Idenokuchi faults also found. It is confirmed that local subsidence continued along the Suizenji fault, which newly appeared during the mainshock in Kumamoto City. Subsidence with westward shift of up to 4 cm/year was also found in Aso caldera. Time constant of postseismic decay ranges from 1 month to 600 days at selected points, but that postseismic deformation during the first epochs or two is dominant at point in the Kumamoto Plain. This result suggests multiple source of deformation. Westward motion around the Hinagu fault may be explained with right lateral afterslip on the shallow part of this fault. Subsidence along the Suizenji fault can be attributed to normal faulting on dipping westward. Deformation around the Hinagu and Idenokuchi faults cannot be explained with right lateral afterslip of Futagawa fault, which requires other sources. Deformation in northern part of Aso caldera might be the result of right lateral afterslip on a possible buried fault.


Introduction
A sequence of large earthquakes struck the city of Kumamoto and its surroundings, the central part of Kyushu, in April 2016, which claimed more than 200 fatalities including disaster-related deaths. This earthquake sequence includes M w 7.0 (USGS 2020) event and several events of M w 6.0 or larger. These earthquakes occurred on and around the Futagawa and Hinagu faults, which are right lateral strike-slip faults with a slightly different strike and meet between Kumamoto City and the Aso caldera (Figs. 1 and 2) (e.g., Asano and Iwata, 2016). The Futagawa fault runs eastward with a strike of N60°E and reaches the Aso caldera. On the other hand, the Hinagu fault trends in the N30 ~ 40°E direction and extends further south of the Yatsushiro city (Geological Survey of Japan, National Institute of Advanced Industrial Science and Technology (hereafter AIST) 2016). The first shock of Mw6.5 is considered to have occurred on the part of the Hinagu fault (Shirahama et al. 2016). Aftershocks are distributed along these faults, but there is difference in pattern of aftershock distribution in western and eastern parts. In eastern part from the epicenter of April 16 shock, aftershocks are aligned tightly along the Futagawa fault, while they are distributed in the fan-shaped area in its west. It is remarkable that there are few aftershocks south of the Futagawa and Hinagu faults. It is also emphasized that northeastern edge of aftershock distribution exceeds the northeastern rim of the Aso caldera.
There are also many reports of surface ruptures off these coseismic faults in the city of Kumamoto and on the western frank of Aso caldera (Goto et al. 2017;Fujiwara et al. 2016;Fujiwara et al. 2017;Kumahara et al. 2016;Toda et al. 2016; Geological Survey of Japan, AIST 2017) (Fig. 2). Most of them are considered to be nontectonic origin. Tsuji et al. (2017) and Fujiwara et al. (2017) reported that surface ruptures in the northern part of Aso caldera were generated by horizontal sliding of blocks or lateral spreading due to strong shaking. Goto et al. (2017) showed detailed distribution of surface ruptures in the Kumamoto Plain. One is the westward extension of the Futagawa fault, which they named the Akitsugawa flexure zone, and other is NW trending multiple traces of surface rupture, Suizenji fault zone, in Kumamoto City. They discussed relationship of them to topography and distribution of pyroclastic flow deposits of Aso volcano. Deformation due to these surface ruptures was also detected with InSAR measurements (Fujiwara et al. 2016;Fujiwara et al. 2017) and it is important to examine their temporal evolution following the earthquake sequence.
Kumamoto City is famous for its abundant groundwater. A lake, which is located close to the western extension of the Futagawa fault, suddenly dried up, which may be associated with the movement of Suizenji fault zone that appeared during the Kumamoto earthquake (e.g., Hosono et al. 2018). Hosono and Masaki (2020) and  reported hydrochemical changes of groundwater during the postseismic period. Groundwater flow may affect movement on the surface. Therefore, observation of surface movement contributes to the understanding of evolution of groundwater flow system in this area.
The Geospatial Information Authority (hereafter GSI) has been monitoring crustal movements with a continuous GNSS network in Japan, called GSI's Earth Observation Network (hereafter GEONET), while the Japan Exploration Agency (hereafater JAXA) has been operating a satellite (the Advanced Land Observing Satellite 2, hereafter ALOS-2) equipped with L-band radar (Phased Array L-band SAR 2, hereafter PALSAR-2). The European Space Agency also operates C-band SAR satellites called Sentinel-1. These sensors detected remarkable coseismic deformation of this earthquake sequence. Many authors processed the data provided by these sensors and presented coseismic fault models. According to these studies, the first shock was a right lateral strike-slip event on the Hinagu fault (Fukahata and Hashimoto 2016;Ozawa et al. 2016;Himematsu et al. 2016;Kobayashi 2017). On the other hand, both Futagawa and Hinagu faults slipped during the Mw7.0 event, but moment release on the Futagawa fault was dominant.
Postseismic deformation usually follows large earthquakes. There are several studies of postseismic deformation following inland earthquakes in Japan mainly using continuous and campaign GNSS data and their origins (e.g., Nakano and Hirahara 1997;Sagiya et al. 2005;Hashimoto et al. 2008;Ohzono 2011;Ohzono et al. 2012;Meneses-Gutierrez et al. 2019). These preceding studies speculated afterslip, viscoelastic relaxation and poroelastic rebound for possible mechanism of postseismic deformation, but they did not incorporate complicated geometry of faults or heterogeneous structure of crust due to the limited spatial resolution. To discuss generation mechanism of postseismic deformation, especially in relation to crustal heterogeneities, spatial resolution is important, but the density of GNSS stations is not high enough to detect detailed spatial distribution of postseismic deformation. Therefore, I must exploit synthetic aperture radar (hereafter SAR) images. Peltzer et al. (1996) discussed postseismic deformation following the 1992 Landers, Calfornia, earthquake using ERS interferograms and clarified relationship between complicated geometry of coseismic faults and poroelastic response. Geology affects groundwater distribution and flow direction. I wonder if there is correlation between the distribution of pyroclastic flow deposit and postseismic deformation. Moore et al. (2017) already studied postseismic deformation following the Kumamoto earthquake based on GNSS and InSAR data till the end of 2016. They mainly discussed large-scale deformation with reference to the viscoelastic structure beneath Kyushu. In this paper, I discuss finer-scale deformation that appeared in the vicinity of coseismic surface ruptures, which may convey invaluable information of property of shallow crust and active faults.

Tectonic setting
Central Kyushu is unique in Japan, because there is a large graben structure across the island. Aso and Unzen volcanoes sit right in its middle (Figs. 1 and 2). Centurylong geodetic surveys revealed N-S extension which is considered to tear the island. This idea seemed partly consistent with the existence of E-W-trending normal faults (Tada, 1984). Recent continuous GNSS observation, however, does not confirm the dominance of N-S extension (e.g., Nishimura and Hashimoto 2006). Now, dextral motion is considered to be appropriate across the Futagawa and Hinagu fault system. Aso volcano is one of the most active volcanoes in Japan and repeated large eruptions many times including at least 4 caldera forming eruptions. The last caldera forming eruption was the largest so far, whose pyroclastic flow deposits, ASO-4 (~ 90 ka BP) covers northern and central Kyushu (Ono and Watanabe, 1985) (Fig. 2). There are thick pyroclastic flow deposits of Pleistocene to Holocene in the surrounding area of the source faults of the 2016 Kumamoto earthquake sequence (#83, 95, 96, 99, 166 in Fig. 2). On the other hand, sedimentary rocks of Holocene are found in the Kumamoto Plain (#1 in Fig. 2). Goto et al. (2017) pointed out that the Suizenji fault zone that appeared during the 2016 earthquake sequence in Kumamoto City is located near the foot of terrace deposit of early-middle-late Late Pleistocence.

SAR images and processing procedure
I utilized ALOS-2/PALSAR-2 images acquired after the largest earthquake on April 16 in the Kumamoto sequence. JAXA made observations with PALSAR-2 for several different directions and modes, but there are not so many images that were acquired from the same orbits and with high frequencies. Among them, I collected strip-map mode images of high spatial resolution of path 23 (P23) of descending orbit and 130 (P130) and 131 (P131) of ascending orbits. Table 1 lists information of observed images with their parameters of observations. Figures 1, 2 and 3 illustrates footprints of images used in this study and temporal changes in perpendicular baselines, respectively. P23 covers the surrounding area of the Futagawa and Hinagu faults and Aso caldera and are frequently observed. It is because this path covers active volcanoes such as Aso, Kirishima, Sakurajima and Kuchinoerabujima. On the other hand, P131 and P130 cover the Kumamoto plain and Aso caldera, respectively, and there is no overlap between P130 and P131. There are 28, 13 and 7 images for P23, P131 and P130, respectively, during the period from April 18, 2016 to December 10, 2018. Perpendicular baselines are shorter than 400 m, which is good enough for interferometry. I did not use ScanSAR images because of their less frequent observations and lower spatial resolution. I did not use other SAR images acquired by other platforms than Sentinel-1, because their shorter wavelength of microwave causes decorrelation in vegetated and mountainous areas and with long temporal separations. I compared result with that of time series analyses of Sentinel-1 images later.
I performed 2-pass interferometry for pairs of collected SAR images with Gamma ® software (Wegmüller and ig. 3 Temporal changes in perpendicular baselines of ALOS-2/PALSAR-2 images used in this study. Blue, orange, and black lines with symbols are of P23, P131, and P131, respectively Werner 1997). For descending images (P23), the boundary between northern and southern images runs across the seismogenic zone of the Kumamoto earthquakes. I concatenated them to retain continuity of phase according to the procedure by Gamma ® . ASTER-GDEM ver.2 is used for the correction of topographic phase and geocoding (Tachikawa et al. 2011). I fixed the first images acquired after the April 16 earthquake as the reference and made interferograms for the pair of this reference and following images. Owing to L-band, coherence is high enough even for the pair with two-year-long separation. L-band SAR used to suffer from ionospheric disturbances and so does the present case. I exploited the technique developed by Gomba et al. (2016), Furuya et al. (2017), Wegmüller et al. (2018) to reduce ionospheric disturbances. We found ionospheric disturbances both in ascending and descending interferograms and sometimes large ramp in corrected interferograms. Therefore, I flattened ionospheric-corrected interferograms and then filtered them before unwrapping. I used the branch-cut technique for unwrapping of filtered interferograms. I stacked unwrapped interferograms for both ascending and descending orbits and converted them to E-W and U-D components.

Correction of ionospheric disturbances
Before discussing the detected surface deformation, it is worth mentioning the correction of ionospheric disturbances. Ionospheric disturbances that appear in interferograms of L-band SAR are considered to be related to medium-scale travelling ionospheric disturbances (MSTID) (e.g., Saito et al. 1998). There may be seasonality of MSTID and dependence on local time (Chen et al. 2019). Observation from ascending and descending orbits is made around mid-night and noon, respectively. Empirically, disturbances appear in ascending interferograms in summer, while those in descending interferograms are recognizable in winter. Figure 4 shows an example of correction of ascending interferograms of April 26 and June 13, 2016. I observed a large disturbance in the middle of the original interferogram (Fig. 4a). Similar disturbances appear Fig. 4 Interferograms produced during the procedure of correction of ionospheric disturbances for the pair of images April 26 and June 21, 2016 of P131-F640. a Original differential interferogram in radar coordinate. b Interferogram for higher band. c Interferogram for lower band. d Double-differenced interferogram; i.e., difference between higher and lower band interferograms. e Estimated ionospheric components. f Non-dispersive components in original interferogram; i.e., subtracted phase of ionospheric components from original one. g Flattened and filtered non dispersive components of interferogram. h Geocoded flattened and filtered nondispersive components of interferogram in interferograms of higher and lower frequencies, but there is slight difference between them (Fig. 4b, c). Double-differenced interferogram shows spatial variation (Fig. 4d), which leads to an estimate of effect of ionosphere (Fig. 4e). Taking difference between original and ionospheric component, I finally obtained ionosphere-corrected interferogram (Fig. 4f ). However, I still recognized significant trend in the azimuth direction. Therefore, I detrended it by fitting polynomial function in two dimensions and filtered it (Fig. 4g). I geocoded filtered ionosphere-corrected interferogram to detect surface deformation (Fig. 4h). In Fig. 4h, there is still a large disturbance that may be attributed to tropospheric disturbance in image of June 21, 2016, because I did not find similar signal in other ionosphere-corrected interferograms ( Fig. 6).
An example of descending interferogram is shown in Fig. 5. Empirically, ionospheric disturbances are considered to be not so serious as those in ascending interferograms. However, it is not always right. In original, higher-and lower-frequency interferograms, I recognized more than three cycles of fringes. Double-differenced interferogram is almost flat, but I have ionospheric disturbance with fairly long wavelength. Furthermore, ionosphere-corrected interferogram still has large trend of three cycles in the azimuth direction, to which I must apply flattening and filtering techniques.

Observed line-of-sight displacements
Owing to repeated acquisitions of ALOS-2/PALSAR-2, I obtained spatio-temporal variation in Line-of-Sight (hereafter LOS) displacements after the occurrence of the Kumamoto earthquake sequence. In this chapter, I discuss characteristics of observed LOS displacements from three different viewpoints; i.e., a) spatial distribution of averaged LOS displacements (Figs. 6,7,8,9 and 10b) profiles of displacements along selected sections (Fig. 11c) time series of LOS displacement at selected points (Fig. 12).
Additional file 1: Figures S1-S3 show all flattened filtered non-dispersive components of interferograms for P131-F640, P130-F650, and P23-F2950-2960, respectively. Close-ups of unwrapped interferograms around the source region and Aso caldera are shown in Figs. 6, 7 and 8, where displacement of GEONET stations during the corresponding period projected onto the LOS directions is also shown. All LOC displacements are referred to GEONET 960700 for the paths 131 and 23, 970833 for the path 130, respectively, considering distance from source faults and coherence around them. Coseismic interferograms are also shown in the top left panels of each figure. Comparison of LOS changes with those at GEONET sites with the same reference is given in Additional file 1: Figure S4. Average LOS change rates with GEONET average velocities are in Fig. 9. Time series of InSAR displacement roughly follow GNSS data at most sites with fluctuations. InSAR displacements in summer tend to depart from that of GNSS, which may be attributed to tropospheric disturbances related to heavy precipitation (Precipitation at Kumamoto is shown in Fig. 12f ). Because no correction of tropospheric disturbances nor temporal smoothing is applied, jumps appear in InSAR time series when there are storms or torrential rains. Furthermore, all interferograms refer to one  specific GEONET site in a scene and local disturbance around it affects entire image. GNSS data is daily averaged coordinate, while InSAR image is an instantaneous one. Therefore, local tropospheric disturbance on interferogram may affect more significantly than GNSS daily coordinates. Discrepancies are large at GEONET sites 950456 and 081169. I suspect soil condition or local topography around these sites affect the movement.
I also compare the present results with that of time series analysis of Sentinel-1 images. I processed Sentinel-1 images during the period from April 20, 2016 to April in 2018 using LiCSBAS developed by Morishita et al. (2020). Additional file 2: Figure S5 shows average LOS displacements of both ascending and descending images. Discrepancies are recognized, but this is attributable to the difference in strategies of analyses. The present result by stacking is the weighted average of changing rate between the first image and others. On the other hand, LiCSBAS calculates average of changing rates of LOS of pairs of consecutive images. Therefore, rapid movement in early stage, if any, may be emphasized in the present result, while LiCSBAS result gives us more slower rates in later stage. Despite this discrepancy, the same features of spatial distribution are recognizable. The most important issue is low coherence in mountainous area on the southeastern side of the Futagawa and Hinagu faults and on northern frank of Aso caldera in Sentinel-1 images. As already known, L-band SAR of ALOS-2/PALSAR-2 gives us higher coherence and can be utilized for the detection of movements. For conversion to E-W and U-D components, the same GEONET stations (960700 and 970833) were fixed in the overlapped area of ascending and descending images. In the following section of spatial variation of deformation, I mainly discuss E-W and U-D components in Fig. 10.

Spatial distribution of average rate of postseismic deformation
Coseismic deformations are also shown at the top left in Figs. 6, 7 and 8. Comparing them with following postseismic interferograms, I confirmed that postseismic deformations are concentrated around the source area of the mainshock. However, spatial pattern is significantly different with each other, especially in ascending interferogram (Fig. 6). Fujiwara et al. (2016) already showed postseismic deformations in early stage, April-May in 2016, with ALOS-2/PALSAR-2 from both ascending and descending orbits. Interferogram from descending orbit is the same as that used in this study (P23; Second left panel of the top raw in Fig. 8). They used pairs of images from a different path with high elevation. There is a little difference in obtained spatial pattern of deformation in ascending interferogram, but the features of obtained postseismic deformations are basically the same. In this study, I put emphasis on their temporal evolution and deformation that arose afterward. Fujiwara et al. (2016) pointed out several spots of significant LOS changes; (1) deformation along the Futagawa fault, especially near the junction with the Hinagu fault (2) deformation around the Suizenji fault (they mentioned as the Suizenji Park), (3) deformation in the Ozu town. In Fig. 12 of Fujiwara et al. (2016), there are many signals in Aso caldera, but they did not mention in detail. I also recognized the same features and that they were amplified in the following 2.5 years (Figs. 6, 7 and 8). They pointed that there is no clear deformation around the outer rim of Aso caldera, where many surface ruptures were observed in coseismic interferograms. I did not observe clear deformation in later interferograms, neither.
The most prominent one is subsidence along the Futagawa fault and its western extension. Fujiwara et al. (2016) measured less than 10 cm displacement near the junction of the Futagawa and Hinagu faults during the first 2 weeks after the mainshock. Subsidence rate exceeding 6 cm/year in this zone is recognized during 2.5 years despite loss of coherence in most part (arrow a in Fig. 10). Another spot of subsidence is found between the junction and Aso caldera (arrow b). Westward shift is also prevailing in this area. There is a surface rupture along another fault, Idenokuchi fault . It is noteworthy that this area of subsidence is bounded by the Futagawa and Idenokuchi faults.
Rapid uplift is found on the south side of the Idenokuchi fault (arrow c). In Fig. 12 of Fujiwara et al. (2016), there is not notable signal in this area. Uplift is also recognized on the north side of the Futagawa fault (arrow d). A zone of slight subsidence and westward shift (arrow e) is surrounded by this uplift zone on the north side of the Futagawa fault. It is interesting that the boundary between these uplift and subsided zones nearly coincides with northern edge of a Pleistocene pyroclastic flow deposits (dark green line).
Westward shift is remarkable on the southeastern side of the Hinagu faults, reaching 6 cm/year (arrow f ). I also see eastward motion of < 2 cm/year around the epicenter. Further west, I observed subsidence in a fan-shaped zone near the coast (arrow g). It is interesting that its southern boundary roughly coincides with the western extension of the Futagawa fault.
I also found significant deformation off Futagawa and Hinagu faults, which is the same as that of Fujiwara et al. (2016). The most remarkable one is a NW-SE trending zone of subsidence of ~ 4 cm/year in the city of Kumamoto (arrow h). Large subsidence was also detected in coseismic interferograms (Upper left panel in Figs. 6 and 8) (e.g., Fujiwara et al. 2016). The zone of this subsidence coincides with the Suizenji fault zone found by Goto et al. (2017). The present results suggest that postseismic deformation also continued around this fault zone during 2.5 years. Several spots of subsidence can be observed in Aso caldera, as well. In the northernmost part of this caldera, coseismic surface ruptures were found (Tsuji et al. 2017;Fujiwara et al. 2017). I detected significant subsidence along these surface ruptures during the postseismic period (arrow i), implying continuing movement associated with these ruptures. Another remarkable motion was found on the northern frank of central cone of the Aso volcano (arrow j), where westward shift is also dominant here. Its northern boundary seems to be aligned along a line trending NE-SW.
Significant eastward motion was found at the central cone of Aso volcano (Fig. 10a). There were small explosions during February to May, 2016, and a significant explosion occurred on October 7-8, 2016 (JMA 2016). This eastward motion may be attributed to this activity.
I also found another small spot of westward shift of ~ 4 cm/year and slight subsidence north of Ozu Town, about 10 km north of the Futagawa fault (arrow k). This deformation was already pointed out by Fujiwara et al. (2016). This zone trends in the WNW-ESE direction, which corresponds to local trend of valley where Pleistocene sedimentary rocks are sandwiched by igneous rocks. I did not see any sign of such deformation in preseismic interferogram (Additional file 2: Figure S5). Therefore, this deformation may have been caused by strong shaking due to the Kumamoto earthquake sequence.

LOS displacement profiles along selected sections
It is important to examine temporal variation in deformation for the discussion of mechanism of postseismic deformation. Because timing and frequency of observations are different between descending and ascending orbits, it is impossible to reduce E-W and vertical components at specific epochs. Therefore, I discuss LOS displacements in this section. For this purpose, I prepared two different views of time series of observed deformation. One is the temporal changes along selected profiles. I sampled LOS change from the area within 0.005° on the both sides of a profile and plotted them shifting according to the time of acquisitions of subsequent images. I chose 7 profiles, shown in Fig. 9b, that run through interesting spots of deformation discussed in the previous section, in which I can also grasp the characteristics of spatial distribution of deformation, especially discontinuities in deformation. 5 sections are along meridians, while 2 sections are in the E-W direction. I emphasize that correlation between LOS displacement and topography is not recognized though some sections runs in the areas of rough topography.
The Sect. 1 is the westernmost profile of LOS change, which runs off the main strand of Futagawa and Hinagu faults but crosses the area of local LOS increase around the Suizenji fault zone in Kumamoto City (Fig. 11a, b). I can see local LOS increase around 32.8°N in both interferograms (vertical line) and another local deformation a little bit north of 32.7°N in descending interferogram (red arrow in Fig. 11b). The former corresponds to local subsidence in Kumamoto City, while the latter is signal on the western extension of the Futagawa fault, i.e., the Akitsugawa flexure zone of Goto et al. (2017). These observations suggest that postseismic deformation occurred not only in the vicinity of coseismic faults but off the source. I notice two steps looking closely at the LOS change around 32.8°N in descending interferogram, implying at least two possible faults there (below SZ).
The Sect. 2 runs just west of the junction of the Futagawa and Hinagu faults (Fig. 11c, d). The LOS increase exceeds 30 cm in descending interferogram, the largest in the entire region under study. I observe sharp changes at the northern boundary of this zone of LOS increase (= subsidence) which corresponds to the Akitsugawa flexure zone (vertical line with AF). Southern half of subsidence zone has gradual change in both interferograms, but is limited by the Hinagu fault (vertical line with HF). Comparing the baseline of the last observation (orange lines), discrete shift of far-field displacement is noticeable on the both sides. The Sect. 3 is a profile running across a smaller local subsidence between the Futagawa and Idenokuchi faults (Fig. 11e, f ). There is a spike-like pattern of spatial distribution of LOS changes around 32.8°N (between vertical lines with HF and FF). Its width is much narrower than that found in the Sects. 2 and 4. There is also a shift in the far-field displacement, which is evident in Fig. 11f.
The Sect. 4 shows temporal evolution of LOS changes along the meridian passing the spot of large subsidence between the Futagawa and Idenokuchi faults (Fig. 11g, h). I recognize sharp changes across these two fault and large LOS increase (= subsidence) between them (vertical lines with IF and FF). This LOS change exceeded 10 cm about 1 year after. It is worth noting that the changes across the Idenokuchi fault are larger and sharper than that across the Futagawa fault especially in descending interferograms (Fig. 11h), which implies afterslip on the Idenokuchi fault is more active than on the Futagawa fault, if any. I also noted that there is another gradual step north of the Futagawa fault (red arrow next right of FF), suggesting a minor buried slip. There is another discontinuous change around 32.9°N (red arrow further right), corresponding to the area of westward shift north of Ozu Town in Fig. 10a. I should note convex pattern of the LOS change in ascending interferograms (double-headed arrow in Fig. 11g), while LOS change along the profile is almost flat in descending ones. This convex pattern of LOS change becomes noticeable about 200 days after.
The Sect. 5 runs across the Aso caldera. A sharp discontinuity is obvious around 33.0°N, just south of the northern caldera rim (RP). This point is located a little north of the surface rupture that was formed during the April 16 shock of Mw7.0 (Fujiwara et al. 2016;Fujiwara et al. 2017;Tsuji et al. 2017). I can notice the differential motion evolved according to elapsed time. There were several step-like pattern of deformation during the first 100 days, but most of them died out and the largest one continued for 2 years. LOS changes with relatively short wavelength of ~ 2 km can be seen in ascending interferogram in caldera floor and central cones, while long wavelength deformation is detected with local LOS increase centered around 32.9°N in descending interferogram (red arrow in Fig. 11j).
The Sects. 6 and 7 are LOS displacement profiles along two parallels. The Sect. 6 runs north of the Futagawa fault and northern part of Aso caldera (Fig. 11k, l). A spike-like change of LOS just east of the caldera rim (left red arrow) is related to coseismic surface rupture, the same signal in Sect. 5. Another notable deformation is rapid LOS increase around 131.2°E in the vicinity of central cone, which is as large as 10 cm (right arrow). This change obviously does not correlate with topography. I also recognize difference in level of LOS change between both sides of this zone in both ascending and descending interferograms.
The Sect. 7 crosses local LOS increase in Kumamoto City, junction of the Futagawa and Hinagu faults, and western frank of the Aso caldera. I can find a remarkable

Time series of LOS displacement at selected points
The other is the time series of LOS changes at selected points, which is easier to understand the decaying history of deformation. We chose 5 points shown in Fig. 9. Because acquisitions were made frequently from descending orbit (P23) and were less from ascending orbits, I examine only time series of descending interferograms. I sampled LOS change rates in an area of 0.005° × 0.005° centered at the selected points and took average. To estimate characteristic time, I fit an exponential decaying function to observed time series; where u is LOS displacement, a and b are constants, t is elapsed time in day from April 16, 2016, τ is characteristic time. Red curves in each panel are estimated decaying time series. It is important to note that the LOS changes till the end of May 2016 are dominant during 2 years at most points, implying much faster motion during this period than this approximation. This fast motion may contribute to the difference between average velocities from stacking of ALOS-2/PALSAR-2 and time series analysis of Sentinel-1. Point A is located in the middle of local LOS increase in Kumamoto City. LOS changes rapidly decayed till the fall of 2016, though there is a fluctuation in 2017-2018 (Fig. 12a). If I fit exponential decaying function, I obtain characteristic time of only 29 days. Total LOS change amounts to ~ 5 cm. Point B is located south of the junction of the Futagawa and Hinagu faults, where westward horizontal motion is dominant around this point (Fig. 10a). This point also shows rapid decay with time constant of ~ 50 days and may have reached ~ 6 cm till the winter in 2016, though scatter is a little bit large (Fig. 12b).
On the other hand, points C-E have longer time constant than the previous points. Point C, located in the large subsidence between the Futagawa and Idenokuchi faults, gradually decayed till the beginning of 2017 with time constant of ~ 230 days (Fig. 12c). In 2017, it is stable at the level of 8 cm increase of LOS, and fluctuated in 2018. Point D is in the middle of uplift area on the western frank of the Aso caldera. During the first 2 weeks, this point moved rapidly, but suddenly was decelerated (Fig. 12d). Then, it continues to move in the same direction (= uplift) with slow decay rate of characteristic time of ~ 980 days.
Point E in Aso caldera shows a similar pattern of temporal change to Point C. Characteristic time is almost the same (~ 210 days) (Fig. 12e). Because these two points are located ~ 20 km away from each other, it may be hard to expect the possible mechanical link.
I add daily precipitation at the Japan Meteorological Agency's (JMA) Kumamoto station in Fig. 12f. Kumamoto area suffered from heavy rain mainly in summer during these 3 years, but the correlation with temporal change in LOS change is not clear at all points.

Trial of afterslip model
There are wide varieties of spatial and temporal characteristics in observed postseismic deformation and it may be difficult to explain them with one mechanism. Because I detected several sharp changes across some coseismic surface ruptures, it is reasonable to examine first to what extent afterslip model can explain observed deformation. For this purpose, I down-sampled average rates of LOS ( Fig. 9) using the quadtree algorithm (Additional file 2: Figure S7), and estimated slip on possible faults by inverting them.
It is obvious that there are at least four or five distinctive deformations in the vicinity of the Futagawa, Hinagu, Idenokuchi, and Suizenji faults and in the Aso caldera. Because there are too many parameters to simultaneously estimate, it is reasonable to separate areas into their surrounding zones as the first step. In this study, I divided dataset into four, considering distance from possible sources (Additional file 2: Figure  S7). Region (1) is the surrounding area of the Futagawa and Hinagu faults. L-band SAR gives us highly coherent phase data in mountainous regions, but I excluded data south of Midorikawa fault and north of 33°N, considering distance from Futagawa and Hinagu faults. I excluded data from the coast of Ariake and Yatsushiro Seas, because this region might have suffered from subsidence due to compaction of artificial land (Fig. 2). I also excluded data in region (2). The region (2) is the vicinity of Suzenji fault. Judging from spatial distribution of LOS displacements, data in about 10 × 10 km 2 wide area were extracted. These areas are covered with P130 and P23 images. The region (3) is Aso caldera, where images of P130 and P23 cover. I excluded data in the area of vicinity of surface ruptures. Tsuji et al. (2017) pointed out that deformation in the vicinity of surface ruptures in Aso caldera may be generated by a source as shallow as 50 m. It is reasonable to exclude them as noise in the following inversion of afterslip.
I applied methods of Fukahata and Wright (2008) and its extension to dual faults (Fukahata and Hashimoto 2016) to down-sampled LOS data.
According to Fukahata and Wright (2008), observed displacement d (N x 1 vector) can be expressed by the function of parameters m (M x 1 vector) and observation error e as below; where f is a vector function including Green's function. m consists of model parameters of faults p (location, length, width, strike, dip) and slip on them a. Thus, (2) can be written as where H is N x M matrix consisting of fault parameters and direction cosine of LOS. Thus, contribution of misfit to the system is where E is covariance matrix of observation data.
Then smoothness condition is added to this system; Finally, solution is obtained by minimizing ABIC in Eq. (20) in Fukahata and Wright (2008). Important parameter is α 2 , which is a hyperparametercontrolled trade-off between data and a priori information (assumption of smoothness). The larger α 2 gives smoother distribution of slip, but residuals between observed data and theoretical displacement become larger. The minimum ABIC can give us an optimal solution with appropriate α 2 .
For regions (2) and (3), I applied Fukahata and Wright's (2008) method, because a single fault is considered to be enough to explain the observed displacements. To reduce contribution of Futagawa and Hinagu faults, I carefully excluded data close to these faults as much as possible. For the region (1), I used the inversion procedure with dual faults by Fukahata and Hashimoto (2016). They modeled the Futagawa and Hinagu faults to explain coseismic deformation. Even with two faults, there are many degrees of freedom. Therefore, I fixed dip angles of two faults as their estimates; 61° and 74° for Futagawa and Hinagu faults, respectively, but length and width were changed considering spatial distribution of deformation. For the Suizenji fault, I assumed as the same strike as the surface ruptures and tried to estimate dip angle and location. In Aso caldera, there is no clear surface expression of faults, but I relied on spatial pattern of observed deformation. I put the modeled fault between zones of eastward and westward motions in Fig. 10a. In these models, slip on the edges except on the surface is fixed. By changing the location and dip angle, I tried to find its optimal model. List of model parameters are given in Table 2. Then slightly changing strike and location of these two faults, I tried to find optimal models that minimize ABIC.
During the course of inversion, covariance matrix is required. Its components are represented as follows assuming gaussian error with zero mean and covariance σ 2 E;  where x i and y i are easting and northing of site i, D is characteristic correlation distance of errors. D = 10 km is often used in many studies. Using covariance matrix with longer correlation length, deformation with short wavelength might be smoothed out. In this study, deformation with shorter wavelength than 10 km is dominant, especially around Suizenji, Futagawa and Hinagu faults. Therefore, I adopted 5 km in this study.
Distribution of ABIC is shown in Additional file 3: Figures S8-S10. Red circle in Additional file 3: Figure S8 and black dots in Additional file 3: Figures S9-S10 indicate optimal models. Overall, optimal models are located close to global minimum. Slight correlation between Xoff and strike for the Hinagu fault is recognized. I selected a model with minimum ABIC for Futagawa and Hinagu fault model, but chose a model with smoother slip distribution than that of minimum ABIC for Suizenji and Aso provisional faults. For model with smaller hyperparameter and minimum ABIC, constraints on slip distribution is weak, which sometimes arises physically unacceptable distribution. Therefore, I selected the second optimal with much smoother distribution of slip. Figure 13 is the compilation of 4 modeled faults with their estimated slip distribution projected onto the surface. Figure 14 shows distribution of estimated slip and its error projected onto a vertical plane along the strike of faults for optimal models. Motion of hanging wall side is shown relative to footwall side. Their residuals and theoretical LOS velocities are shown in Fig. 15 and Additional file 3: Figure S11, respectively. Larger residuals than 2 cm/year are found around the eastern tip of the Futagawa fault. Negative residuals are also seen along the central and western part of the Futagawa fault. These large residuals suggest complexity of deformation and possible other sources than afterslip.
Slips are concentrated in the depth shallower than 10 km for all models. Estimated errors are not larger than 8 cm/year. Optimal model for the Futagwa and Hinagu faults is very closely located to the surface ruptures. On the Futagawa fault, there are three main areas of large slip with a couple minor patches (Fig. 14a). Easternmost patch has left lateral slip of ~ 20 cm/year, which is against coseismic slip; e.g., Figure 4 in Fukahata and Hashimoto (2016). Normal faulting of up to 12 cm/year is dominant in central patch. This patch is located about 5 km east of the junction. The Fukahata and Hashimoto's (2016) model shows normal fault component in its eastern part. These left lateral and normal slip arises from westward motion and local subsidence on the north side of the Futagawa fault. Obviously, these motions cannot be created with right lateral slip on this fault. Therefore, it is not considered that westward motion around the eastern tip of Futagawa fault was caused by its afterslip. The westernmost patch with the largest slip is located west of the junction of the Futagawa and Hinagu faults. Right lateral slip exceeds 30 cm/year. As there is no significance slip in the coseismic model of Fukahata and Hashimoto (2016), this slip may be generated by stress concentration at the edge of coseismic slip. Hinagu fault has two patches of large slip (Fig. 14b). Northern patch is closely located to the westernmost patch of the Futagawa fault. Its normal faulting may be related to subsidence near the junction of these two faults, which also suggests interaction between two faults. Furthermore, considering geological condition there, this subsidence might be caused by the compaction of soil. The southern patch on the Hinagu fault has right lateral slip of ~ 20 cm/year. Its peak is estimated at the depth of ~ 3 km and slip almost reaches the surface. Observed displacements show clear discontinuity (e.g., Figure 9a) and creeping of surface ruptures were confirmed in this region (e.g., Shirahama et al. 2016). Therefore, right lateral afterslip is highly possible on this patch of the Hinagu fault. This model fails to explain subsidence between the Futagawa and Idenokuchi faults (Additional file 3: Figure S11a, b). Incorporation of Idenokuchi fault adds more complexities in inversion, which is beyond the present capability of inversion scheme. There might be contribution of compaction of soil in this area. Future work that incorporates these complexities is desirable. Figure 14(c) shows slip distribution of the Suizenji fault, where normal faulting of less than 10 cm/year was detected. Dip angle was estimated 64°, which is consistent with that used in stress calculation by Goto et al. (2017). Upper margin of this fault corresponds to one of the strands of surface rupture. Slip is concentrated in the depth range of 2-8 km. However, slip in very shallow part is negligible, which causes underestimate of observed displacements (Additional file 3: Figure S11c, d). Figure 14d is slip distribution of a provisional fault in Aso caldera. Dip angle was estimated 55° southward. I also made similar calculation with northward dipping fault model, but obtained ABIC is larger. Right lateral slip is dominated with its peak at a depth of ~ 4 km and maximum slip reaches 20 cm/year. This motion may cause subsidence in northern frank of central cone of Aso and uplift on the southwestern rim of caldera. Subsidence around the central cone cannot be explained by this model, which may be related to volcanic activity of Aso (Additional file 3: Figure S11e, f ).

Discussions
I presented the results of analysis of ALOS-2/PALSAR-2 images acquired after the 2016 Kumamoto earthquake sequence. In this section, I point several pros and cons in the present study and problems to be resolved in the future.

Efficiency of L-band SAR
Thanks to long wavelength of PALSAR-2, coherence is high even for pairs with longer temporal separation than 2 years (Fig. 5). The longest separation is 2.7 years (April 18, 2016 andDecember 10, 2018), but high coherence is obtained enough to detect deformation even in mountainous regions. Recently Sentinel-1 images are being used to study crustal deformation because its recurrence is 6 or 12 days and large amount of image of the same area have been already accumulated. However, temporal decorrelation is strong especially in vegetated area (e.g., Morishita et al. 2020), and it is difficult to obtain deformation with a single pair of images with long temporal separation. This is one of the biggest advantages of L-band SAR. I expect continuous accumulation of PAL-SAR-2 images as long as possible.

Correction of ionospheric disturbances
Ionospheric disturbances were observed in both ascending and descending interferograms, and their correction with Split Beam interferometry was effective especially for ascending interferograms (e.g., Fig. 4). It is interesting that distribution of ionospheric disturbance is different between ascending and descending interferograms (Figs. 4 and 5). Local time of acquisition is around midnight for ascending orbit, while observations are made around noon from descending orbit. This difference may be the cause of different pattern of ionospheric disturbances that appear in L-band interferograms. Chen et al. (2019) discusses variation of characteristic parameters of MSTID such as period, wavelength and phase velocity observed over Hongkong, and mention that wavelength of MSTID is slightly longer in daytime of spring, autumn, and winter than that in night in spring and summer, though the difference seems marginal.
To verify the ionospheric correction, people consider use of GNSS TEC. Comparison of ionospheric disturbances by GNSS and InSAR, however, is not straightforward. First, the timing of observation is different, even though recent continuous GNSS observation is made at the interval of 1 s. Second, incidence angle and azimuth are not the same. Coincidence of LOS of SAR and GNSS satellites might be rare. Finally, distribution of GNSS sites is sparse for this purpose. As shown in Fig. 4, wavelength of ionospheric disturbance is much shorter than length of one scene (~ 70 km) in the azimuth direction. Average spacing of GEONET in Japan is 20 ~ 25 km. It is hard to reproduce detailed distribution of ionospheric disturbances in interferogram with GNSS data. Therefore, I followed the method by Wegmüller et al. (2018) to verify the results.

Comparison of postseismic deformation with preceding inland earthquakes in Japan
I detected postseismic deformation following the 2016 Kumamoto earthquake sequence. The maximum displacement exceeded 20 cm near the junction of the Futagawa and Hinagu faults (Fig. 11d). I observe several spots of larger LOS changes than 10 cm (Figs. 6, 7 and 8). Are these large postseismic displacements special for the Kumamoto earthquake? Observations of postseismic displacements were made for previous inland earthquakes in Japan as listed in Additional file 1: Table S1. Postseismic displacements are definitely dependent on size of and distance from the mainshock. Therefore, I should compare those with mainshock of similar size to the Kumamoto earthquake. Of course, it is not suitable to strictly compare results because of sparse distribution of GNSS sites around the epicenter, but it may give some insights into characteristics of postseismic deformation.
First, I compare with strike-slip events. The first example is the Kobe earthquake in 1995 (M JMA 7.3, M w 6.9; all following Mws are from USGS (2020)). Nakano and Hirahara (1997) reported postseismic displacements detected by campaign Global Positioning System (hereafter GPS) surveys and early GEONET. They detected about 2.5 cm displacement at Iwaya station, northern tip of Awaji Island, which is closely located to the epicenter (~ 2 km), till the end of 1995. Hashimoto (2017) detected subsidence between two active faults along the NE extension of the source fault of the 1995 Kobe earthquake with ERS-1/2, Envisat and ALOS/PALSAR. Its maximum was less than 1 cm/yr, which is one order smaller than that of Kumamoto case. Sagiya et al. (2002) detected only ~ 3 cm postseismic displacements at the station right above the aftershock area after the 2000 Western Tottori earthquake (M JMA 7.3, M w 6.7) during half year. In case of Kumamoto earthquake, a GEONET site 021071 west of the Hinagu fault recorded 8 cm displacement during 2 years. Considering moment magnitude, it is acceptable that postseismic deformation of the Kumamoto earthquake is larger than Kobe and Tottori events.
What about thrust events? Takahashi et al. (2005) observed postseismic displacement of 3 cm or larger during about 2 months after the 2004 Niigata Chuetsu earthquake (M JMA 6.8, M w 6.6). For the 2007 Noto Peninsula earthquake of M JMA 6.9 (M w 6.7), only 2 cm displacements were observed by campaign GPS surveys by Hashimoto et al. (2008). After the 2007 Chuestu Oki earthquake (M JMA 6.8, M w 6.6), Ohta et al. (2008) detected postseismic displacements less than 2 cm at a GEONET station during ~ 50 days. Although distance from the epicenter is larger than 15 km, distance from the edge of aftershock area is much shorter. Ohzono (2011) showed postseismic deformation of up to 13 cm at a GEONET station located ~ 11 km from the epicenter during 800 days after the 2008 Iwate-Miyagi Nairiku earthquake of M JMA 7.2 (M w 6.9). Ohzono (2011) also detected ~ 11 cm postseismic deformation at their original site 2.5 km from the epicenter. Moment magnitude of earthquakes other than Iwate-Myagi event is much small than the Kumamoto earthquake, though observation periods are short to compare. Iwate-Miyagi earthquake has as large displacement as the Kumamoto earthquake, implying correlation with magnitude of mainshock.
Postseismic deformation, however, may be controlled not only by magnitude of mainshock, but also by geometrical relationship between the source and observation points, local geological conditions, flow of groundwater, etc. These factors should be pursued in the future.

Possible correlation with geological structure
Considering these different features of postseismic deformations between the Kumamoto earthquake and other inland earthquakes in Japan, it is speculated that there may be different characteristics in the Kumamoto area. Spatial pattern of deformation and distribution of pyroclastic flow deposits seem to be correlated with each other (Fig. 10). For example, large LOS increase in the Aso caldera is located in the region covered with igneous rock of Cenozoic Quaternary Holocene. Uplift zone north of the Futagawa fault corresponds to the area of early Late and Late Pleistocene volcanic rocks. Local subsidence is distributed in a narrow zone about 10 km north of the Futagawa fault. This zone corresponds to the area of middle-late Late Pleistocene (Fig. 10). These observations imply that the age of igneous and sedimentary rocks might affect the response to coseismic loading. It is important to re-examine postseismic deformation following previous inland earthquakes from this viewpoint.

Temporal characteristics of postseismic deformation
Postseismic deformation following the 2016 Kumamoto earthquake sequence may have decayed during 2 years, though it may still continue in some areas (Fig. 11c). Although observation periods are short for other inland earthquakes discussed above, they may have decayed with short time constant as well. It is noteworthy that the LOS changes during the first epochs or two are dominant in the whole time series and cannot fully be explained with a simple exponential decay. A possible cause of deformation with short time constants is poroelastic rebound or movement of groundwater. As Hosono et al. (2018) reported, water level rapidly dropped in the lake near the Suizenji fault, suggesting fast flow of groundwater.
I also found deformations that arose with delay such as concave pattern in Fig. 11g, acceleration of motion on the southeastern side of the Futagawa fault in Fig. 11m. The former is westward motion on the north side of the Futagawa fault in Fig. 9a. These delayed onsets of deformation might not be related to afterslip. Recently,  proposed a model of flow of groundwater in this area. They performed hydrogeochemical study of groundwater and suggested that precipitated water came down from surface ruptures on the western frank of Aso caldera and flew toward the Kumamoto Plain. They also implied rise of water level on the north side of the Futagawa fault and in the Kumamoto Plain. The uplift detected in the present study might be related to such a phenomenon.

Deformation in and around Aso Caldera
There are other issues to be solved by the future works. For example, uplift and westward motion on the western frank of the Aso caldera cannot be explained by afterslip on the Futagawa or Idenokuchi faults (Fig. 10). At present, I would like to rule out the possibility of magma intrusion or large-scale landslide. This area is about 10 km away from central cones. I cannot accept the magmatic activity there. As shown in the preceding section, flow of groundwater may be one of candidates. Large-scale landslide may not be candidate, neither, because uplift is dominant. The InSAR technique, however, has little sensitivity to displacement in N-S direction. There might be possibility that movement dominantly occurred in N-S direction. It may be a good idea to incorporate image acquired with different incidence angles and directions, which help resolve three-dimensional displacements.

Conclusions
I processed ALOS-2/PALSAR-2 images acquired after the 2016 Kumamoto earthquake sequence with correction of ionospheric disturbances and revealed spatiotemporal variation in LOS changes during 2 years. I could draw conclusions below: 1. L-band SAR gives us high coherence enough to reveal surface deformation even in vegetated or mountainous area for pairs of images acquired more than 2 years. 2. Ionospheric disturbances are seen both in the ascending and descending images, but spatial characteristics may be different each other. 3. Notable features of postseismic deformations are as follows: a. Deformation earthquake exceeds 10 cm during 2 years at some spots in and around Kumamoto city and Aso caldera. b. Westward motion of ~ 6 cm/year was dominant on the southeast side of the Hinagu fault, while westward shift was detected on both side of the Futagawa fault. The area of this westward motion has spatial correlation with distribution of pyroclastic flow deposits. c. Significant uplift of 4 cm/year was found around the eastern Futagawa fault and on the southwestern frank of Aso caldera. d. Sharp changes were found across several coseismic surface ruptures. e. Rapid subsidence between Futagawa and Idenokuchi faults was also detected.