Seasonal and transient surface displacements in the Kumamoto area, Japan, associated with the 2016 Kumamoto earthquake: implications for seismic-induced groundwater level change

The 2016 Kumamoto earthquake sequence on April 14 (Mw 6.2) and April 16 (Mw 7.0) altered the regional groundwater level. To better understand the relationship between groundwater level change and surface displacement, we estimated surface displacement in the Kumamoto area (Japan) using persistent scatterer interferometry from 19 ALOS/PALSAR images acquired between January 7, 2007 and March 5, 2011, 28 ALOS-2/PALSAR-2 images acquired between April 17, 2016 and December 10, 2018, and 113 Sentinel-1 images acquired between May 26, 2016 and December 30, 2018. Our estimation shows that transient surface displacement occurred following the 2016 Kumamoto earthquake sequence, together with seasonal surface displacement that was not detected from the 2007–2011 images. We suggest that a portion of the transient displacement occurred via groundwater drawdown through new ruptures that formed owing to the 2016 Kumamoto earthquake sequence and sediment compaction. Seasonal surface displacements detected after the 2016 Kumamoto earthquake sequence are linked to groundwater level variations.


Introduction
Substantial changes in groundwater levels following large-scale earthquakes have been extensively reported (e.g., Montgomery and Manga 2003;Manga et al. 2012;Manga and Wang 2015). The groundwater level change mechanism is related to changes in pore pressure (e.g., Jonsson et al. 2003), permeability (e.g., Elkhoury et al. 2006), and seismogenic dilatant cracks (e.g., Sibson and Rowland, 2003). A notable example is the 2016 Kumamoto earthquake sequence, which began on April 14 (M w 6.2) followed by a mainshock (M w 7.0) on April 16 along the Hinagu and Futagawa fault systems in southwestern Japan (Fig. 1a, b). Coseismic groundwater drawdown was observed in the 30-45 days following the mainshock using high-resolution groundwater level network data Nakagawa et al. 2019;Kagabu et al. 2019) in both confined and unconfined aquifers in central Kumamoto where newly formed surface fissures had also been reported . In contrast, groundwater levels continuously increased in the confined aquifer of the eastern terraces in the first 2 weeks following the mainshock ) and immediately increased after the mainshock in the confined aquifer along the western coast. Nakagawa et al. (2019) showed that confined aquifer groundwater levels changed over a Ishitsuka et al. Earth, Planets and Space (2020) 72:144 period of several weeks to months; whereas, the groundwater in shallow unconfined aquifers rapidly recovered to initial levels. Kagabu et al. (2019) suggested that the postseismic increase in water levels observed in the eastern terraces of Kumamoto were owing to mountain water release.
Hydraulic head variations can lead to the compaction and expansion of granular aquifer skeletons, which can induce surface displacement. To quantify these features, synthetic aperture radar interferometry (InSAR) (Massonnet and Feigl 1998) and time-series InSAR analysis, such as persistent scatterer interferometry (PSI) (Ferretti et al. 2000(Ferretti et al. , 2001, have been used to determine the spatial pattern and temporal evolution of surface displacement. Valuable information has, thus, been obtained regarding surface displacement characteristics and groundwater level response (e.g., Bell et al. 2008;Osmanoglu et al. 2010;Erban et al. 2014;Khakim et al. 2014;Ishitsuka et al. 2014;Chaussard et al. 2014;Chen et al. 2017;Hu et al. 2018). InSAR and time-series InSAR analysis have also been proven useful for understanding how earthquakes influence the relationship between surface displacement and groundwater level (e.g., the 1992 Landers earthquake, Johnston et al. 1995; the 2000 Iceland earthquake, Jonsson et al. 2003; the 2011 Tohoku earthquakes, Ishitsuka et al. 2017). Previous studies in the Kumamoto area characterized the groundwater level change associated with the 2016 Kumamoto earthquake sequence using high-resolution groundwater level network data Nakagawa et al. 2019;Kagabu et al. 2019). The Kumamoto area is, thus, highly suitable for comparing surface displacement patterns and groundwater level changes associated with earthquakes.
In this study, we estimated surface displacement in the Kumamoto area before and after the 2016 Kumamoto earthquake sequence using PSI analysis of the Advanced Land Observing Satellite (ALOS)/phased array-type L-band synthetic aperture radar (PALSAR) images, ALOS-2/PALSAR-2 images, and C-band Sentinel-1 images. We further combined the satellite data analysis with groundwater level variations to constrain the mechanisms of surface displacement following the 2016 Kumamoto earthquake.

Hydrogeological setting
The Kumamoto area, located in central Kyushu in southwestern Japan, has abundant groundwater resources that are recharged along the western flank of the Aso volcano  (1995). Black square is the location of the city of Kumamoto, and black circles indicate the locations of GNSS stations (Kumamoto (North station) and Jonan (South station)). Dashed rectangle is the area of interest. G1, G2 an G3 indicate the locations of groundwater wells plotted in Fig. 8. Black lines indicate surface ruptures associated with the 2016 Kumamoto earthquake based on Toda et al. (2016). A black triangle is the reference location of PSI analysis. The focal mechanism for the foreshock (M w 6.2) and mainshock (M w 7.0) was derived by F-net managed by the National Research Institute for Earth Science and Disaster Resilience (Asano and Iwata, 2016) and flow toward the Ariake Sea under topographic constraints ( Fig. 1b) (AIST-GSJ 2014). Most of the water used in the city of Kumamoto is derived from groundwater (Kumamoto Prefecture, 1995). The surface geology of the western part of the Kumamoto area consists of Quaternary alluvium (Fig. 2a); whereas, the eastern part is characterized by pyroclastic deposits from the Aso volcano that erupted over four periods between 30 and 90 ka as well as terrace deposits (Fig. 2a). The aquifers consist mainly of pyroclastic deposits from the Aso volcano. The upper aquifers are either unconfined or confined and composed of the most recent pyroclastic deposits and marine beds; whereas, the lower aquifers are confined and consist mainly of older pyroclastic deposits (AIST-GSJ 2014). The altitude of the upper aquifer is approximately 0-200 m at the terrace and − 50 to 0 m in the lowland area; whereas, that of the lower aquifer is approximately 0-100 m at the terrace and − 250 to 0 m in the lowland area (National Institute of Advanced Industrial Science and Technology (AIST)-Geological Survey of Japan (GSJ) 2014) (Fig. 2b, c). The spatial groundwater flow patterns in the upper and lower aquifers are generally similar (AIST-GSJ 2014).

Synthetic aperture radar images
For the pre-earthquake period, we used 19 ALOS/PAL-SAR images taken on a descending orbit between January 7, 2007 and March 5, 2011 (Table 1, Additional file 1:  Table S1). For the post-earthquake period, we used 143 images taken over a period of 32 months, including (1) 28 ALOS-2/PALSAR-2 images obtained on a descending orbit between April 18, 2016 and December 10, 2018, (2) 65 Sentinel-1 images obtained on a descending orbit between July 1, 2016 and December 30, 2018, and (3)   The ALOS-2 images are expected to be less affected by temporal decorrelation owing to the longer wavelength (~ 23.8 cm); whereas, the higher temporal sampling of the Sentinel-1 images maintains coherence despite the smaller wavelength (~ 5.6 cm). The aim of using PAL-SAR-2 and Sentinel-1 descending images is to demonstrate the reliability of the estimated surface displacement by comparing the spatial and temporal patterns obtained from different satellite images. We used the Sentinel-1 ascending images to understand the quasi-vertical and quasi-EW (East-West) displacement by combining Sentinel-1 images from descending and ascending orbits.

Persistent scatterer interferometry and characterization of time-series displacement pattern
Persistent scatterers (PSs) are point-like scatterers that are coherent for all baseline conditions (Ferretti et al. 2000(Ferretti et al. , 2001 and, therefore, show little or no decorrelation noise that typically appears in pixels without PSs. Manmade objects or bared surfaces are physical examples of PSs that coherently reflect microwaves. Accurate surface displacements can, thus, be estimated from the temporal phase changes of PSs. PSI analysis provides the processing flow that identifies PSs and estimates surface displacement using SAR images (Ferretti et al. 2000(Ferretti et al. , 2001Werner et al. 2003;Hooper et al. 2004;Kampes 2006;Ferretti 2014). PSI analysis also includes strategies to mitigate atmospheric phase contributions based on their spatial and temporal characteristics. The PSI processing sequence used in this study was based on the approach of Kampes (2006). We used the radar interferometry calculation tool (Ozawa et al. 2016) to generate single-look differential interferograms for the PALSAR images and GAMMA software (Wegmüller and Werner 1997;Wegmüller et al. 2016) for the Sentinel-1 and PALSAR-2 images. Other PSI processing steps were implemented using Fortran and Matlab (version R2019a).
Further implementation details are provided in Ishitsuka (2015) and Ishitsuka and Matsuoka (2016), and other PSI processing applications are discussed in Ishitsuka et al. (2016a;2016b;. We selected all interferometric pairs with respect to a single reference image (Ferretti et al. 2001) (Table 1). The PS candidates (PSCs) were selected based on an amplitude dispersion index (Ferretti et al. 2001) of 0.30. We then created single-look differential interferograms of the interferometric pairs at the PSC pixels. We used a 10-m mesh digital elevation model (DEM) provided by the Geospatial Information Authority of Japan to remove topographic phase components. Phase stability was subsequently evaluated using phase coherence (Ferretti et al. 2001). We used a coherence threshold of 0.85 for the PALSAR images and 0.80 for the PALSAR-2 and Sentinel-1 images. Because the phase contribution of the DEM error was estimated during the phase coherence evaluation, we subtracted the phase contribution from the interferometric phases of each PS pixel. A minimum cost flow algorithm was used for phase unwrapping (Costantini and Rosen, 1999). The line-of-sight time-series displacement m est was estimated by combining all of the interferograms using the least squares method. To reduce the atmospheric phase contribution, we simultaneously smoothed the temporal phase variations based on the assumption that the atmospheric phase contribution in the time domain is within the high-frequency end of the spectrum (Eq. (1)) (Schmidt and Bürgmann 2003). Spatial filtering is also often used to mitigate atmospheric contributions; however, we did not use spatial filtering in this study because surface displacement can occur after a large earthquake with a variety of spatial characteristics.
where G is a matrix representing the temporal distribution of the SAR acquisitions with values of 1 or 0

Table 1 Basic information of SAR images used in this study
All of acquisition dates of SAR images used in this study was provided in a supplemental information (Additional file 1: Tables S1 and S2). (r, az) for the spatial resolution represents the spatial resolution of the slant-range and azimuth direction, respectively depending on the data availability at a given time, d is a displacement dataset obtained from the PSI analysis with vector dimensions of 18 × 1 in the case of the ALOS/ PALSAR images, m represents the time-series surface displacement of the 19 SAR acquisition time, and τ is a weighting factor that indicates the degree of temporal smoothing of the time-series surface displacement and determined based on Akaike's Bayesian information criterion (ABIC) (Akaike 1980). We assume that atmospheric contributions in the temporal domain appear as a Gaussian distribution; ABIC can, thus, be calculated as: with where H is the designed matrix of d/dt in Eq.
(1), m is the optimal value of the time-series surface displacement, and N , P , and M are the number of images, rank of matrix H , and rank of m , respectively. To determine τ , we calculated the ABIC of τ over a range of 0.001-1.0 and interval of 0.001 using the time-series displacement at a PS pixel to obtain the optimal τ that minimizes ABIC. More specifically, we first calculated the optimal τ values of randomly selected 10% of the total PS pixels from the range and interval. The τ value in Eq. (1) was then determined from the averaged values. As a reference point, we assumed that no surface displacement occurred in the northern section of the study area (Fig. 1b).
To validate the estimated surface displacement, we compared the line-of-sight surface displacements derived by the PSI analysis with those obtained from global navigation satellite system (GNSS) observations of GEONET, which is a nationwide GNSS network operated by the Geospatial Information Authority of Japan (Hatanaka et al. 2003) and provides daily coordinate values. Observations were made for the Kumamoto GEONET station with respect to the Jonan GEONET station within the study area (Fig. 1b). To reduce random noise, we applied a moving average filter to the GEONET displacement time series using a window size of 5 days before and after the date of interest (i.e., 11 days) except in the case of the GNSS observations on April 18, 2016 (a PALSAR-2 acquisition date), which were collected 2 days after the mainshock of the 2016 Kumamoto earthquake sequence. To obtain the reference surface displacement of the date, we averaged the daily coordinate values of the date in reference and the 5 subsequent days. To compare with the GNSS observations, we averaged the surface displacement of the PSI analysis using PS pixels within 300 m from a GNSS station. (2) To better understand the surface displacement characteristics, we modeled the time-series displacements as the sum of sinusoidal and exponential functions. The sinusoidal function was used to model the seasonal displacement over a period of 1 year. The exponential function was used to model the transient displacement. We estimated four unknown parameters ( α , θ , β , and k ) using the least squares method: where α represents the magnitude of the seasonal displacement and θ is the time shift at the beginning of the seasonal uplift in the sinusoidal function relative to the acquisition time of the first SAR image. The factor sin(θ) was introduced such that U (0) = 0 . Here, β is the magnitude of the exponential function and k is the temporal scaling time. The scaling time represents the day when 63.2% of the total transient displacement had occurred. We estimated θ and k using a grid search algorithm over a search range of − 180 to 180° and 5° intervals for θ and search range of 20-10,000 with an interval of 10 for k . During the grid search, we simultaneously estimated α and β using a least squares approach. We then determined the optimal values of α , θ , β , and k that minimize the square difference between the modeled and PSI timeseries displacement. The temporal model in Eq. (4) consistently satisfies U(0) = 0, where t = 0 is the acquisition date of the first SAR image.
Transient displacement after an earthquake can occur in the vertical and/or horizontal directions and reflects the displacement mechanism. We, therefore, converted the first year of transient surface displacement, modeled as an exponential function in Eq. (4), into quasi-vertical and quasi-EW directions by combining the Sentinel-1 descending and ascending line-of-sight displacement data. This processing is referred to as 2.5-D analysis (Fujiwara et al. 2000). The acquisition geometry of the SAR images is insensitive to the NS (North-South) direction; displacements in the quasi-NS direction can, therefore, not be retrieved using the analysis. Owing to data availability, the Sentinel-1 ascending images used for the 2.5-D analysis were from between November 16, 2016 and June 15, 2018. To combine with descending images in a corresponding time span, we used the Sentinel-1 descending images from November 22, 2016 to June 9, 2018 for the ascending images. We then converted the interpolated line-of-sight displacements in the quasi-vertical and quasi-EW directions. The dip angle in the quasi-vertical direction was 81.7° toward the south and the strike angle of the quasi-EW direction was 91.3° from the north.   Fig. 1b (Toda et al. 2016) first SAR image ( β exp(−t/k) − 1 and t = 365 days in Eq. (4)) and seasonal (sinusoidal) displacement magnitude ( α in Eq. (4)) using PSI analysis of the ALOS/PAL-SAR images acquired between January 2007 and March 2011. Figure 3c-e shows examples of time-series displacement and best-fit model in Fig. 3a, b. Although the best-fit models are not accurately fitted to the observed PSI, the results demonstrate that nearly no significant steady surface displacement occurred during the time period in most of the study area (Fig. 3a). The time-series displacement in the northern section of the study area fluctuated more than in other areas (Fig. 3e); however, the seasonal displacement magnitude was less than 0.25 cm during this period (Fig. 3b). Figure 3f shows a comparison of the PSI analysis and GNSS data at the Kumamoto station. PSI and GNSS line-of-sight displacement are compared at the GNSS station of Kumamoto using the GNSS station of Jonan as a reference point. The absolute differences of the PSI analysis and GNSS data were 0.15 cm for the seasonal displacement magnitude and 0.30 cm for the transient displacement during the first year from the first SAR image.

Surface displacement during the 32 months following the 2016 Kumamoto earthquake sequence
The spatial distribution of the transient displacement during the first year and seasonal displacement magnitude of the PALSAR-2 and Sentinel-1 descending images are shown in Figs. 4a, b, 5a, b. Examples of the time-series displacement and best-fit model are shown in Figs. 4c-e and 5c-e. A comparison of the line-of-sight time-series displacement between the PSI analysis and GNSS data at the Kumamoto station with respect to the Jonan station are shown in Figs. 4f and 5f. The absolute differences of the seasonal displacement magnitude between PSI analysis and GNSS data are 0.072 cm (PALSAR-2 descending images) and 0.027 cm (Sentinel-1 descending images), and the absolute differences of the transient displacement during the first year are 0.12 cm (PALSAR-2 descending images) and 0.19 cm (Sentinel-1 descending images).
The PSI analysis of the PALSAR-2 and Sentinel-1 descending images shows a similar surface displacement pattern. For the transient surface displacement, displacements toward the satellite were found mainly in the northern section of the Futagawa-Hinagu fault system. Surface displacements away from the satellite were identified particularly around the central part of the study area (indicated by I in Figs. 4 and 5) and around the western coastal area (indicated by II in Figs. 4 and 5). The magnitude of the transient displacement during the first year in I was approximately < − 1.1 cm (where the negative sign indicates the direction away from the satellite), and that in II was approximately < − 2.0 cm (Figs. 4a,   5a). Seasonal surface displacements with a magnitude of ~ 0.5 cm were identified in the northern and central parts of the study area (III in Figs. 4 and 5) in both the PALSAR-2 and Sentinel-1 results. The area of the relatively large seasonal displacement is almost the same in the descending PALSAR-2 and Sentinel-1 images; the areal distances involved are several kilometers (Figs. 4b,  5b). The seasonal displacement with a magnitude of less than 0.25 cm was estimated between January 7, 2007 and March 5, 2011; however, the magnitude of the seasonal displacement in the period following the 2016 Kumamoto earthquake sequence was approximately 0.5 cm (Fig. 5b), which is larger than these earlier estimates and beyond the error derived from the comparison between the GNSS and PSI displacements (0.15, PALSAR images; 0.072, PALSAR-2 images; 0.027 cm Sentinel-1 descending images).
We then combined the line-of-sight surface displacements derived from the Sentinel-1 descending and ascending images using 2.5-D analysis. As described in "Persistent scatterer interferometry and characterization of time-series displacement pattern" section, because the temporal span of the Sentinel-1 ascending images used in this study (November 16, 2016-June 15, 2018 was shorter than that of the Sentinel-1 descending images (April 18, 2016-December 10, 2018), we used 43 Sentinel-1 descending images acquired from November 22, 2016 to June 9, 2018, which correspond to the time span of the ascending images, to compare with the Sentinel-1 ascending image analysis. The transient displacements during the first year are shown in Fig. 6a for the Sentinel-1 descending images and Fig. 6b for the Sentinel-1 ascending images. Although the first SAR image used in Fig. 6a, b was acquired approximately 7 months after the 2016 Kumamoto earthquake, the surface displacements away from the satellite around the center of Kumamoto area (I in Fig. 6a, b) and coastal area (II in Fig. 6a, b) were found as well as in Figs. 4a, b, 5a and b. These displacements from the Sentinel-1 ascending images were validated by a comparison with the GNSS data, which yielded an absolute difference of the seasonal displacement magnitude of 0.038 cm and of the transient displacement during the first year of 0.11 cm (Additional file 1: Fig. S1). Figure 6c, d illustrate the transient surface displacements in the quasi-vertical and quasi-EW directions estimated using the analysis. We found that the transient displacements in I and II of Figs. 4 and 5 occurred in the vertical downward (subsidence) direction (Fig. 6d). The surface displacements along the Futagawa fault included both the vertical and the EW directions (Fig. 6c, d).
We further conducted cross-comparisons of our PSI processing of the PALSAR-2 and Sentinel-1 descending and ascending images at pixels where the horizontal displacement was negligible based on the quasi-EW components. For pixels with a quasi-EW displacement of < 0.05 cm, we assumed a negligible horizonal displacement and averaged the line-of-sight displacement within a 100-m radius of these pixels. The averaged line-of-sight

Implications of surface displacement and groundwater level change
After a large earthquake, surface displacements along the mainshock fault include post-seismic afterslip, viscoelastic relaxation, and poroelastic displacement (e.g., Feigl and Thatcher 2006). Post-seismic afterslip occurs because a mainshock fault releases strain energy, even after the main earthquake. Post-seismic afterslip, therefore, generally occurs along the mainshock fault. Conversely, viscoelastic relaxation occurs because of stress release in the lower crust and upper mantle associated with mainshock fault slip. Moore et al. (2017) and Politz et al. (2017) estimated surface displacement owing to the viscoelastic relaxation in Kumamoto using surface displacement data for 91 days and 8.5 months following the earthquake. They assumed that the depth of the ductile region was approximately 20-80 km based on the regional Moho depth pattern (e.g., Katsumata, 2010). Although there may be spatial heterogeneity in the viscoelastic structures, this viscoelastic relaxation depth causes high wavelength surface displacement (> 20 km). Poroelastic displacement by earthquakes has been observed around mainshock faults for the 1992 Landers earthquake (Johnston et al. 1995;Fialko 2004) and 2000 southern Iceland earthquake (Jonsson et al. 2003). The poroelastic displacement model shows that groundwater migrates because of volumetric stress/strain (pore pressure) changes caused by an earthquake. This allows the spatial patterns in the displacement model to be correlated with the volumetric stress/strain associated with the earthquake. We identified post-seismic surface displacements in the Kumamoto area that cannot be explained by the above model (e.g., afterslip, viscoelastic relaxation, and poroelastic displacement owing to static stress changes). These surface displacements show that the spatial patterns do not agree with the spatial characteristics expected from afterslip and viscoelastic relaxation, and do not correlate with the coseismic volumetric strain shown in Fig. 7. Such areas include the transient subsidence around the central part (I in Figs. 4,6) and western coastal area (II in Figs. 4,5,6) of the study area, and the seasonal displacement around the northern and central parts of the study area (III in Figs. 4,5). We subsequently discuss how these displacements differ from the above post-seismic displacements (e.g., afterslip, viscoelastic relaxation, poroelastic displacement) and their implications. Hosono et al. (2019) analyzed the groundwater drawdown immediately after the 2016 Kumamoto earthquake sequence around the central part of Kumamoto (I in Figs. 4,5,6), and a newly formed rupture zone was identified by the InSAR analysis of Fujiwara et al. (2016) (Fig. 7). Groundwater drawdown was observed in both confined and unconfined aquifers with a maximum decline of 4.74 m. Groundwater levels recovered and returned to the background levels within 45 days. Hosono et al. (2019) demonstrated that the groundwater level reductions can be explained by downward groundwater transfer (for a volume of 2.7 × 10 7 m 3 ) through the newly formed rupture system, which extended downward to 6.8 km. The area of the ground subsidence around the central part of Kumamoto (I in Figs. 4,5,6) corresponds to the newly formed rupture zone and post-seismic groundwater drawdown area. Because the magnitude of the groundwater level drop cannot be explained by volumetric stress changes caused by the 2016 Kumamoto earthquake sequence   (Fig. 7), the process of poroelastic deformation is not a viable mechanism for the ground subsidence estimated in this study.

Transient subsidence around the central and western areas
In the western coastal areas where ground subsidence was identified, groundwater level increases of up to 1 m were observed after the 2016 Kumamoto earthquake sequence . The subsidence and associated groundwater level changes can be explained by liquefaction and consolidation of the sediment. Large-scale liquefaction after an earthquake has often been reported by InSAR, for example, in the 2011 Tohoku earthquake, Japan (M w 9.0) (Ishitsuka et al. 2012) and the Fairview (M w 5.0), Pawnee (M w 5.8), and Cushing (M w 5.0) earthquakes in Oklahoma (Barnhart et al. 2018). The alluvial delta in the western parts of Kumamoto consists of interbedded clay, sand, and gravel layers, and flood deposits  Fujiwara et al. (2016). Blue lines indicate surface ruptures associated with the 2016 Kumamoto earthquake based on Toda et al. (2016) overlying the pyroclastic deposits are partly distributed along the rivers. The aforementioned explanation is supported by the work of Wakamatsu et al. (2017) who carried out a field reconnaissance study, analyzed aerial photographs, and found that soil liquefaction had occurred at numerous locations in the western coastal area.

Seasonal displacements in the northern part of the Kumamoto area
Groundwater levels in the Kumamoto area show a seasonal pattern in response to seasonal rainfall variations. The rainfall in the study area generally increases from June to August (Fig. 8) (JMA 2020). Correspondingly, the groundwater level starts to increase in June and decrease in October or November (Fig. 8). Figure 8c shows the monthly precipitation and groundwater levels around the seasonal surface displacement areas between January 2016 and December 2018. The seasonal surface displacement in the line-of-sight direction identified by the PSI analysis are in Fig. 8a, b, and show that the seasonal surface displacement correlates with the groundwater level changes. To examine the possible amount of precipitation induced the seasonal displacement, we calculated the average monthly precipitation during the period and obtained a value of 174.3 mm/year, which is similar to the average precipitation during the PALSAR period (159.6 mm/year). This possibility of seasonal displacement owing to the increased precipitation can be dismissed from the observed groundwater levels. Compared with the seasonal groundwater level changes during the PALSAR-2, Sentinel-1, and PALSAR periods, the general pattern of seasonal changes in groundwater levels did not vary significantly. We therefore suggest that there are two possible scenarios that can explain the increase in seasonal surface displacement. One is that the magnitude of the seasonal groundwater levels increased in a deeper part of the area where groundwater levels have not been measured. The other possibility is that the sensitivity of the surface displacements to groundwater level changes was altered by the earthquake.
With regard to the first possibility of groundwater level variations in the deeper section, the groundwater variation after an earthquake can be explained by a permeability enhancement. For example, seismic vibration due to an earthquake might remove colloids and bubbles that clog pores, thereby changing the permeability. Such changes in permeability have been previously observed in association with the 1952 Kern Country earthquake (Manga et al. ), 1989 Loma Prieta earthquake (Rojstaczer and Wolf 1992), and 1995 Kobe earthquake (Tokunaga 1999;Sato et al. 2000).
The storage coefficient must be taken into consideration to investigate the second possibility with regard to the alteration of the surface displacement sensitivity to groundwater level change. The storage coefficient ( S ), or storativity of a confined aquifer or aquitard, is linked to the volume of water released/accumulated from a unit area when the hydraulic head declines/increases by a unit amount. The coefficient can be measured using surface displacement and groundwater level changes. The coefficient is a dimensionless quantity that ranges between zero and the effective aquifer porosity (Chaussard et al. 2014), and is given by: where d and h represent the surface vertical displacement and hydraulic head change, respectively. If it is assumed that no gas phase is present in the pores of a (5) S = d h , Fig. 8 a, b Examples of time-series displacements and best-fit models from a PALSAR-2 (Fig. 4e) and b Sentinel-1 descending (Fig. 5e) images. c Monthly precipitation in the Kumamoto area (JMA 2020) during January 2016-December 2018 (the period corresponding to ALOS-2/PALSAR-2 and Sentinel-1 images used in this study). Dashed lines are the groundwater level data with respect to the average level. The locations of G1, G2 and G3 are mapped in Fig. 1b confined aquifer, then S is simply a function of the matrix compressibility ( c m ) and porosity ( ν ) (Riley 1969), where: where ρ w is the pore fluid density, c f is the pore fluid compressibility, g is the acceleration of gravity, and b is the thickness of the saturated system. If it is assumed that pore fluid compression/expansion can be ignored, the S represents the purely skeletal characteristics of an aquifer system. It is plausible to assume that ρ w , g , c f , and b are identical before and after an earthquake. Consequently, assuming that the enhanced permeability described above did not occur and the porosity ν did not significantly vary, one plausible mechanism to explain why the seasonal displacement occurred is the change in compressibility ( c m ) of the aquifer. Alteration of the characteristics of an aquifer by an earthquake generally occurs in unconsolidated materials (Wang et al. 2004;Cox et al. 2012;Quigley et al. 2013;Manga and Wang 2015). This suggests that areas with high values of α in Eq. (4) generally represent underlying terrace deposits, which might contain large amounts of soft material such as unconsolidated sediments. Furthermore, areas with a larger seasonal displacement are found mainly in the northern section of the analyzed area. This occurs because these areas are covered mostly by terrace deposits, whereas the southern and southeastern parts of the analyzed area contain volcanic and metamorphic rocks (Fig. 2a). Resistivity maps obtained using magnetotelluric methods show that the northern areas had a lower resistivity than the southern and southeastern areas (Asaue et al. 2012). Therefore, it is likely that the geological conditions led to the variation in the aquifer properties, which are related to the characteristics of the observed seasonal surface displacement.

Conclusions
We mapped the surface displacements in the Kumamoto area using PALSAR, PALSAR-2, and Sentinel-1 images. Transient and seasonal surface displacements were observed following the 2016 Kumamoto earthquake sequence. In addition to the transient displacements found along the mainshock fault, subsidence was also detected around the central part of the study area where groundwater drawdown has previously been reported to be associated with coseismic ruptures. Moreover, post-seismic seasonal displacement was identified in the northern section of the Kumamoto area, showing a higher magnitude than that estimated for the periods between January 2007 and March 2011. We suggest that the seasonal displacements are correlated with the seasonal variations of groundwater (6) S = ρ w g(c m + νc f )b, level, and can have thus been enhanced by the 2016 Kumamoto earthquake due to changes in aquifer permeability and/or compressibility. Surface displacement associated with coseismic ruptures and the alteration of permeability/compressibility has rarely been reported. Our estimation of surface displacement, therefore, provides important guidelines to better understand surface displacement following a large earthquake, especially in relation to changes in groundwater level.

Acknowledgements
The ALOS/PALSAR and ALOS-2/PALSAR-2 level 1.1 images used in this study were obtained through cooperative research with the Japan Aerospace Exploration Agency (JAXA), the PALSAR Interferometry Consortium to Study our Evolving Land Surface (PIXEL), and the Earthquake Research Institute of the University of Tokyo. A part of the ALOS-2/PALSAR-2 images used in this study were provided as part of the 2nd research announcement of the earth observations by JAXA. The groundwater level data used in this study were measured by several public institutions and corrected by the working group of the Japan Association of Groundwater Hydrology. Parts of the figures were created using Generic Mapping Tools (Wessel et al. 2013). We also express our appreciation to an editor and two anonymous reviewers for their insightful comments.
Authors' contributions KI analyzed the SAR images and modeled the seasonal and the transient displacement. KI, TT, and WL interpreted the estimated surface displacement. MK and JS analyzed the groundwater level data. KI produced the initial draft of the manuscript, which was improved by TT, WL, MK, and JS. All authors discussed the research issues and results. All authors read and approved the final manuscript.

Funding
This study was supported in part by the Japan Society for the Promotion of Science (JSPS) KAKENHI (Grant Nos. 17K14724 and 20K15219) and the Japan Construction Information Center Foundation (JACIC) (Grant No. 2019-2).

Availability of data
The PSI results analyzed during the current study are available from the corresponding author on reasonable request. The PALSAR and PALSAR-2 SAR images used in this study are available from the Japan Aerospace Exploration Agency, but restrictions apply to their availability, which were used under license for the current study and are therefore not publicly available. The PAL-SAR and PALSAR-2 SAR images are available from the authors upon reasonable request and with permission from the Japan Aerospace Exploration Agency. The Sentinel-1 SAR images were provided by the European Space Agency and can be freely downloaded from https ://scihu b.coper nicus .eu/.