Nationwide urban ground deformation monitoring in Japan using Sentinel-1 LiCSAR products and LiCSBAS

Ground subsidence in urban areas is a significant problem because it increases flood risk, damages buildings and infrastructure, and results in economic loss. Continual monitoring of ground deformation is important for early detection, mechanism understanding, countermeasure implementation, and deformation prediction. The Sentinel-1 satellite constellation has globally and freely provided frequent and abundant SAR data and enabled nationwide deformation monitoring through InSAR time series analysis. LiCSAR, an automatic Sentinel-1 interferometric processing system, has produced abundant interferograms with global coverage, and the products are freely accessible and downloadable through a web portal. LiCSBAS, an open source InSAR time series analysis package integrated with LiCSAR, enables users to obtain the deformation time series easily and quickly. In this study, spatially and temporally detailed deformation time series and velocities from the LiCSAR products using LiCSBAS for 73 major urban areas in Japan during 2014–2020 were derived. All LiCSBAS processing was automatically performed using predefined parameters. Many deformation signals with various temporal and spatial features, such as linear subsidence in Hirosaki, Kujyukuri, Niigata, and Kanazawa, episodic subsidence in Sanjo, annual vertical fluctuation in Hirosaki, Yamagata, Yonezawa, Ojiya, and Nogi, and linear uplift in Chofu were detected. Unknown small nonlinear uplift signals were found in Nara and Osaka in 2018. Complex postseismic deformations from the 2016 Kumamoto earthquake were also revealed. All the deformation data obtained in this study are available on an open repository and are expected to be used for further research, investigation, or interpretation. This nationwide monitoring approach using the LiCSAR products and LiCSBAS is easy to implement and applicable to other areas worldwide.


Introduction
Ground subsidence poses significant threats, especially in highly populated areas, because it increases flood risk, damages buildings and infrastructure, and results in economic loss. Many cities have experienced anthropogenic subsidence due to the fluid extraction of groundwater, oil, or natural gas, combined with rapid urbanization, or industrialization. For example, in Japan, significant subsidence has been observed in Tokyo (> 4 m in total) since the beginning of the twentieth century and spread to alluvial plains across the country in the 1950s (Yamamoto 1995; Tokyo Metropolitan Government 2019; Ministry of the Environment 2020). To prevent rapid anthropogenic subsidence, laws to control groundwater usage were enacted in 1956. As a result of new policies, ground subsidence has remarkably decreased since the 1970s.
Continual ground deformation monitoring is important for early detection of subsidence, understanding the deformation mechanisms, implementing countermeasures, and predicting future deformation. Leveling surveys and the Global Navigation Satellite System (GNSS) have traditionally been used to monitor ground subsidence with high accuracy (~mm). However, these techniques obtain only pointwise deformation data and are costly and time consuming. Synthetic Aperture Radar (SAR) Interferometry (InSAR) has also been used for monitoring since the 1990s, and the InSAR time series analysis that can achieve higher accuracy (~mm/ year) than that of conventional InSAR (~cm) has substantially enhanced the capability of slow subsidence monitoring since the 2000s (e.g., Bell et al. 2008;Chaussard et al. 2013). In Japan, many case studies using conventional InSAR and the time series analysis with Japanese Earth Resources Satellite-1 (JERS-1), Advanced Land Observing Satellite (ALOS), and ALOS-2 data have detected various subsidence signals (e.g., Nakagawa et al. 2000;Morishita et al. 2010;Nonaka et al. 2020). However, the temporal sampling frequency was generally low, and therefore, it was difficult to determine the (nonlinear) temporal deformation signal features.
The C-band Sentinel-1 constellation launched by the European Space Agency (ESA) in 2014 and 2016 has been able to provide frequent and abundant data globally and freely (Potin et al. 2019). Several public nationwide deformation monitoring services using Sentinel-1 data are in operation or planned in Denmark, Germany, Italy, the Netherlands, and Norway (Crosetto et al. 2020), as well as commercial services for Denmark, Japan, France, and California in the USA (Bischoff et al. 2020).
The Looking Inside the Continents from Space (LiCS) project by the Centre for the Observation and Modelling of Earthquakes, Volcanoes, and Tectonics (COMET) also uses Sentinel-1 data for large-scale deformation monitoring. The Looking into Continents from Space with Synthetic Aperture Radar (LiCSAR) system, an automatic Sentinel-1 interferometric processing system, makes the products freely accessible and downloadable through a web portal (COMET 2020;Lazecký et al. 2020). LiCSBAS, an open source InSAR time series analysis package integrated with LiCSAR, can almost automatically conduct time series analysis using freely available LiCSAR products . Weiss et al. (2020) derived the velocity and strain rate fields over the approximately 800,000 km 2 Anatolian region from roughly 300,000 LiCSAR products (i.e., unwrapped interferograms) using LiCSBAS and revealed large-scale tectonic deformation as well as localized strain accumulation and subsidence with high spatial resolution. The LiCSAR products and LiCSBAS enable users to obtain the deformation time series quickly and easily, using only free and open resources without a highperformance computing facility that is required for general InSAR processing.
In this study, the deformation time series and velocities in major urban areas in Japan were derived from the LiCSAR products using LiCSBAS. The majority of the processing was performed automatically with predefined parameters. Many subsidence signals, as well as other types of deformation with various spatial and temporal features, were detected with high accuracy. All processed data (i.e., deformation time series, velocities, and noise indices) are available on an open repository (see "Availability of data and material"), can be easily displayed by the time series viewer in LiCSBAS or geographic information system (GIS) software, and are expected to be used for further research, investigation, or interpretation of the deformation signals. This study also demonstrates the feasibility of nationwide deformation monitoring using the LiCSAR products and LiCSBAS and their applicability to other areas worldwide.

Data
Sentinel-1 unwrapped and geocoded interferograms and corresponding coherence data (0.001°resolution) were processed by LiCSAR and published on the COMET-LiCSAR web portal (COMET 2020). There are 34 LiCSAR frames (18 and 16 in ascending and descending orbits, respectively) covering Japan (Fig. 1, Table 1) over an observation period from November 2014 to April 2020. The naming convention of the LiCSAR frame ID is OOOP_AAAAA_BBBBBB, where OOO denotes the number of the relative orbit, P denotes the orbital direction (i.e., descending (D) or ascending (A)), AAAAA is a colatitude identifier, and BBBBBB identifies the number of included bursts (Lazecký et al. 2020). While the general acquisition interval was primarily 24 days before the end of 2016, it shortened to 12 days after 2017 because of the start of the Sentinel-1B operation. Each frame has 240-430 (350 on average) interferograms derived from 70-140 (110 on average) acquisitions (i.e., three or four interferometric pairs for each acquisition with the preceding acquisitions; Lazecký et al. 2020). In total,~12, 000 interferograms and~4,000 acquisitions were used. The spatial baselines are almost always within 300 m.
Notably, the acquisition start dates are not consistent for all frames because of the non-uniform observation strategy of Sentinel-1. For example, the earliest acquisition in 127A_06256_070603 and 170A_04675_131008 is April 30, 2016 and January 4, 2016, respectively, which means the deformation from 2014-2015 cannot be measured in these frames. There are also long acquisition gaps halfway through the observation period in some tracks and frames. For example, almost no data were ac-

Target areas
Most (~70%) national lands in Japan are covered with forests (Figs. 2 and S1). Dense vegetation prevents Cband SAR data from obtaining sufficient coherence in an interferogram, even with short temporal baselines (Rosen et al. 1996;Morishita et al. 2020), and tends to yield unwanted unwrapping errors. To avoid the impact of unwrapping errors and save the processing time and storage, this study focuses on urban areas. The 73 target urban areas are defined based on the densely inhabited district (DID; Ministry of Land, Infrastructure, Transport and Tourism 2015) metric, and compliance with at least one of the following conditions (Table 2, Figs. 2 and S1): (1) Having a population of ≥ 100,000.
(4) Leveling surveys have been continuously conducted (Ministry of the Environment 2020).
The area ID with triple digits consists of a double digits prefectural ID (01-47, e.g., 01 Hokkaido) and a single digit sub-number. As clear deformation has been detected outside of the DID in Hirosaki (ID: 022), Kanto (ID: 131), and Kumamoto (ID: 431) by past studies (e.g., Morishita et al. 2010;Nonaka et al. 2020;Hashimoto 2020), their target areas were expanded to measure complete deformation signals. All the defined areas except Naha (ID: 471) had at least one ascending and descending data set each, and some had two data sets due to overlapping areas between adjacent tracks (Fig. 1). All 191 available data sets for the defined clipped areas were processed (Table 2). ishita et al. 2020) was used in this work. The LiCSBAS processing flow involves five steps for data preparation (step 0-1 to 0-5) and six steps for time series analysis (step 1-1 to 1-6). Each step is briefly discussed below. For further detailed information, please refer to Morishita et al. (2020).

LiCSBAS (v1.3) for the InSAR time series analysis (Mor
Step 0-1: Download LiCSAR products Step 0-2: Convert file format Step 0-3: Tropospheric noise correction (optional) Step 0-4: Mask low coherence areas in the unwrapped interferograms (optional) Step 0-5: Clip a specified rectangular area of interest (optional) Step 1-1: Quality check and identify bad interferograms Step 1-2: Loop closure check and identify bad interferograms Step 1-3: Small baseline inversion Step 1-4: Calculate standard deviations of the velocity Step 1-5: Mask time series Step 1-6: Filter (and deramp) time series The tropospheric noise correction using the Generic Atmospheric Correction Online Service (GACOS) data (Yu et al. 2018) in step 0-3 was not applied except to data from Kanto (ID: 131) and Kumamoto (ID: 431) because the clipped areas were small enough to neglect the long-wavelength tropospheric noise and focus on flat areas where the topography-correlated stratified tropospheric noise is minimal ( Figure S1). To reduce the impact of unwrapping errors in non-urban areas, pixels with an average coherence of ≤ 0.1 in the unwrapped data were masked in step 0-4 before the subsequent time series analysis. Thereafter, all the unwrapped data and corresponding coherence images were clipped to the predefined areas in step 0-5 (Figs. 2 and S1; Table  2). The unwrapped coverage and coherence threshold in step 1-1 were set to 0.5 and 0.06, respectively. To mask noisy pixels based on several noise indices in step 1-5, strict thresholds to retain reliable pixels were set. The gap number threshold in the small baseline (SB) network was set to zero so that no gap was allowed in the network and the obtained time series was highly robust. The threshold numbers for loop error, maximum connected time length, and interferograms without loops were set to 3, 3 years, and 10, respectively. The temporal filter width in step 1-6 was set to 0.15 years (55 days), which is long enough compared to the sampling interval and short enough to retain seasonal components of the deformation time series. The default settings were used for other parameters. All parameter settings are listed in Table S1.
Kanto (ID: 131) is very wide (~160 km) and was significantly affected by long-wavelength postseismic deformation from the 2011 Mw 9.0 Tohoku earthquake (Suito 2018;Morishita et al. 2020). Because the purpose of this study is to measure localized urban ground deformation, deramping was applied with the best-fitting second order polynomial to remove the long-wavelength deformation component in step 1-6. The Mw 7.0 Kumamoto earthquake occurred in Kumamoto (ID: 431) on April 16, 2016. Coseismic interferograms tend to be decorrelated due to very large surface deformation (> 2 m) along the seismogenic Futagawa Fault, and they cause a gap in the SB network (Fujiwara et al. 2016;Himematsu and Furuya 2016;Kobayashi et al. 2019). To avoid the gap, interferograms from after the earthquake were used, as they detect postseismic deformation. The number of acquisitions before the earthquake was not sufficient for deriving reliable velocities and time series; thus, the pre-seismic data were not used in this study.
During the LiCSBAS processing, a primary stable reference point is automatically selected in step 1-6 as follows. The root mean square (RMS) of cumulative deformation time series after the spatio-temporal filter in step 1-6 at each pixel with reference to the median deformation for each epoch is calculated. Subsequently, a pixel with the minimum RMS is selected as a primary stable reference point. In this study, a common reference point for an area among different frames was further set to compare or combine the results from multiple frames by computing the average RMS among all available frames for an area and selecting a pixel with the minimum averaged RMS. Locations for a stable GNSS station outside of the significant deformation area were manually selected as the common reference point in Sanjo (ID: 153) and Kumamoto (ID: 431), because most of the processed area was deformed, and the automatically selected point was not stable. The amplitude and time offset (i.e., phase) of annual sinusoidal components from the filtered deformation time series with the common reference point for each area were also estimated simultaneously with the velocity and constant term using least-squares . Moreover, the multiple line-of-sight (LOS) velocities were decomposed into vertical and east-west (EW) components, neglecting north-south (NS) deformation because of low sensitivity to the NS component (Wright et al. 2004;Motagh et al. 2017;Fuhrmann and Garthwaite 2019). Morishita et al. (2020) estimated the uncertainty of the deformation time series and velocities derived from LiCSBAS by comparing to GNSS and leveling and concluded that the accuracy was < 1 cm/epoch and~2 mm/year with a 24-day interval and a time period longer than 2 years. Because the targets in this study were urban areas where interferometric coherence is generally high and severe thresholds are used for masking, the expected accuracy would be comparable to or better than 1 cm/epoch and 2 mm/year.
In this section, subsidence signals seen in rice paddy fields are not discussed because they may not represent true deformation and result from systematic bias caused by fading signals (Ansari et al. 2020; see Section 4.1). Many subsidence signals were observed in reclaimed sites and mining areas as expected (e.g., 231 Nagoya, 341    (Suito 2018). Subsidence, uplift, and annual fluctuations in DID presumably associated with groundwater or natural gas were observed in approximately 10, 10, and 3 areas, respectively, and the details of some of the significant cases are described in the following subsections.
As the observed subsidence in urban areas is slow (≤ 15 mm/year except for Sanjo), and only four areas have the rate of ≥ 10 mm/year, overall subsidence in Japan seems to be quite inactive due to strict regulations on groundwater usage (Table 2). Other nationwide studies using InSAR detected rapid subsidence in developing countries, such as Indonesia (8-22 cm/year; Chaussard   (Table 2). Therefore, the seasonal subsidence is likely attributed to the  In contrast, Nogi had the opposite seasonal component; negative (subsidence) peaks in summer. Settlement and groundwater gauges also measured the consistent seasonal fluctuations, and the subsidence signals in summer were considered to be caused by high summer agricultural groundwater usage (Tochigi Prefecture 2011). Considerably large amounts of subsidence (~70 cm) occurred from 1977 to 1997 in Nogi, although the subsidence speed has since decelerated, and almost no subsidence has been detected recently (Tochigi Prefecture 2011). In 1996, this area had a severe drought and experienced the largest amount of subsidence in Japan measured by leveling surveys and JERS-1 InSAR (Nakagawa et al. 2000). This implies that this area has contraction-prone soil layers.

Fast linear subsidence and segregated annual fluctuation in Hirosaki
Fast subsidence (~40 mm/year) was detected by the ALOS data in 2007-2009 in the Hirosaki plain (ID: 022; Morishita et al. 2010). Most of this area is covered with light vegetation, and several unwrapping errors exist especially in the descending LiCSAR products. The masking step in LiCSBAS (step 1-5) returned too sparse a result to capture the entire deformation signal ( Figure  S10a); therefore, unmasked results were used here. They contain more noise than the masked ones; however, the subsidence signals were still larger than the unmasked noise, as shown below.
Two rapid subsidence basins were found outside of the DID in the Hirosaki plain (~40 mm/year; Fig.  10a). They had contractional EW components as well as large vertical components, as expected for a subsidence basin ( Fig. 10b; Samieie-Esfahany et al. 2009;Fuhrmann and Garthwaite 2019). The deformation time series at P0221 and P0222 indicated almost linear subsidence (Fig. 10e). Interestingly, a significant annual component (amplitude of~15 mm) was detected in the northern DID (Goshogawara City, P0223) adjacent to the northern subsidence basin (Fig. 10c, d). The time series showed that the negative (subsidence) peaks occurred in winter, as in cities shown in the previous subsection, which implies this annual deformation is also attributed to the groundwater usage for snow melting.

Episodic subsidence in Sanjo
Morishita et al. (2020) detected~4 cm episodic subsidence in the winter of 2017/2018 in Sanjo (ID: 153), likely caused by unusually heavy snowfall and related over pumping of groundwater, as well as annual components of negative (subsidence) peaks in winter. Here, three independent data sets (one ascending and two descending) are available and showed consistent measurements of the episodic subsidence and annual components (Figs. 11 and S11). Notably, the time series at P1531 potentially indicated a slight uplift tendency after the episodic subsidence in the winter of 2017/2018. Continual monitoring will reveal whether the subsided ground surface is rebounding or not over a longer timeframe.

Linear uplift in Chofu, Kanto
While uplift signals were rarely observed in this study, Chofu (ID: 131, Tokyo in Kanto; Figs. 12a and S4) data showed slow uplift (~5 mm/year). To minimize the atmospheric noise contribution and obtain the precise time series, I manually set the reference point near the uplift region for the time series plot (P1313ref). The time series at P1313 indicates that the uplift rate was nearly constant (Fig. 12b). The Tokyo Metropolitan Government (2019) reported that the groundwater level is rising, and uplift has been observed by leveling in this area. The shape of the uplift was not isotropic; a V-shaped edge with relatively large deformation gradient was seen in the southern uplift area, which might be controlled by subsurface geological boundaries.

Nonlinear uplift in Nara
Two distinct small uplift signals with a distance of~10 km (P2911 and P2912) were detected in the Nara prefecture (ID: 291; Figs. 13a and S12). Because the EW deformation component also showed small horizontal extension (Fig. 13b), and the deformation areas were not correlated with Normalized Difference Vegetation Index (NDVI; i.e., the land use; Fig. 13c) or topography ( Figure  S12b), the deformation source was not on the surface but at depth. The deformation time series showed approximately 1 cm of nonlinear uplift in 2018 in both areas (Fig. 13d). The uplift timing at the western P2912 seemed to be a few months earlier than that of the eastern P2911. These facts suggest that the uplift signal sources may not be common, but are related. In this area, no earthquakes with > Mw 3.5 and depth of < 50 km occurred during the observation period of Sentinel-1 (Japan Meteorological Agency 2020). No active volcanoes or known faults were around. It is unlikely that these uplift signals stemmed from atmospheric noise, because the spatial extent is only~2 km and~5 km, and the deformation time series showed clear temporal correlation. Currently, the sources of these nonlinear uplift signals are unknown, and further investigation is required.

Nonlinear uplift and annual fluctuation between faults in Osaka
In the northern part of the Osaka prefecture, the Arima-Takatsuki fault zone and Rokko-Awaji fault zone are nearly parallel and form a triangular zone between them (The Headquarters for Earthquake Research Promotion 2020a, 2020b; Fig. 14a). In this study, a small uplift signal with an extent of~3 km was found in the western edge of the triangle zone, accompanying small westward deformation (Figs. 14a, b and S13). The deformation time series at P2711 showed~1.5 cm nonlinear uplift in 2018, which resembles signals from a slow slip event on a fault (Fig. 14f). Adjacent to the uplift zone, an area with relatively large annual deformation exists in the east (Fig. 14c, d, f). As these areas are highly urbanized (Fig. 14e), the deformation time series should be very accurate. Hashimoto (2016) revealed the deformation history in the triangle zone from 1992 to 2010 using various SAR data. After the 1995 Mw 6.8 Kobe earthquake, subsidence started throughout the zone. The subsidence became concentrated in two spots, accelerating after 2003, and lasting at least until 2010. The western spot is adjacent to the deforming zones detected in this study (Fig.  14). The deformation style changed with each decade, among small blocks with the width of a few km in the triangle zone. The mechanism of these varying deformations, whether related to groundwater or tectonic activity, is not clear. On June 18, 2018, a Mw 5.6 earthquake with a depth of 13 km occurred 25 km east of the triangle zone along the Arima-Takatsuki fault zone (Fig. 14a). While significant coseismic deformation was not detected, a small and shallow slip may have been triggered along the fault zone (Fujiwara et al. 2020). However, considering the temporal feature of the uplift in the triangle zone and its long distance from the epicenter and shallow slip, there seems to be no direct relationship between the earthquake and the uplift.

Postseismic deformation of the 2016 Kumamoto earthquake
The Mw 7.0 Kumamoto earthquake occurred on April 16, 2016, and associated large coseismic deformation (> 2 m) and small displacement lineaments were observed around the Futagawa and Hinagu faults (e.g., Himematsu and Furuya 2016;Fujiwara et al. 2016; Fig  S14b). While postseismic deformation was detected by ALOS-2 InSAR (Fujiwara et al. 2016;Kobayashi 2018;Hashimoto 2020), the detailed temporal evolution was not revealed due to infrequent observations (2-4 times/year in general). This study shows the spatially and temporally detailed postseismic deformation revealed by the Sentinel-1 LiCSAR products and LiCSBAS, which is essential for understanding the complex postseismic deformation mechanism (Figs. 15 and S14).
Only the postseismic acquisitions (i.e., after April 16, 2016) of Sentinel-1 were used. The reference point was manually selected at the 950464 GNSS station, which is far from the deformed region (Fig. 15a, b). Two ascending data sets had long acquisition gaps between the summer of 2018 and 2019, which prevented the observation of detailed temporal behavior. However, the long interferograms spanning the gap had sufficient coherence, and therefore no gap in the SB network (i.e., no effect of temporal constraint in the SB inversion, see Morishita et al. 2020) was generated.
Subsequently, the LOS deformation time series was compared with the GNSS data at two GNSS stations (950465 and 021071; Fig. 15f). Each time series largely agreed with each other, and all standard deviations of the difference between InSAR and GNSS were less than 1 cm, which is consistent with the uncertainty estimated in Morishita et al. (2020). The EW component (Fig. 15b) showed that the entire region, except the west side of the Hinagu Fault, moved westward. This is opposite from the coseismic deformation caused by the right-lateral strike-slip and is considered to be attributed to viscoelastic relaxation of the lower crust and upper mantle (Pollitz et al. 2017;Kobayashi 2018).
The vertical component (Fig. 15a) showed more complex signals than the EW component. Around the west coast of the western extension of the Futagawa Fault (P4311), significant subsidence was detected. The temporal behavior showed exponential decay, which implies that afterslip is the dominant source of this subsidence (Fig. 15c). Similar downward velocities were also observed by ALOS-2 (Kobayashi 2018;Hashimoto 2020). Because most of this region is covered with rice paddy fields ( Figure S14a), this subsidence signal may include the bias related to fading signals (see the Section 4.1). Around the Suizenji Park (P4312), exponentially decaying localized subsidence was detected (Fig. 15d). Dip-slip coseismic deformation occurred in this region and is considered to have been triggered by the motion of the main seismogenic Futagawa Fault and associated extensional strain field (Fujiwara et al. 2016(Fujiwara et al. , 2020. The time constant of the exponential decay seemed to be shorter than that of P4311, and the subsidence almost stopped around 2018. The consistent temporal behavior using only the Sentinel-1 descending data and infrequent ALOS-2 data was also reported in Bischoff et al. (2020) and Hashimoto (2020), respectively.
Large postseismic deformations were observed across the Hinagu Fault (Fig. 15e), while the coseismic deformation was much smaller than that of the Futagawa Fault (Himematsu and Furuya 2016;He et al. 2019;Kobayashi et al. 2019). Because the ascending and descending LOS deformations showed opposite trends, the horizontal, rather than the vertical, component was dominant. The detected EW deformation is consistent with right-lateral slip along the Hinagu Fault. Considering the northnortheast strike direction, NS deformation would be much larger than EW deformation, though InSAR has almost no sensitivity to the NS component.
Almost no measurement points were obtained in the southeast side of the faults because the region is covered with dense vegetation. Although L-band ALOS-2 can retain sufficient coherence in these areas, the temporal sampling frequency is much lower and the detailed temporal evolution of the postseismic deformation cannot be revealed (Hashimoto 2020). The next-generation Lband SAR missions with high observation frequency such as ALOS-4 (Motohka et al. 2019) or NASA-ISRO SAR (NISAR; Jet Propulsion Laboratory 2020) will help clarify the complex postseismic deformation, even in such densely vegetated regions. . Because LiCSAR products consist of temporally short interferometric pairs, the impact of the systematic bias could be significant. In this study, several subsidence signals were observed outside of the DID, and the distribution of the subsidence signals was consistent with rice paddy fields (Figs. 16a, b and S5-S14), which almost all show some degree of subsidence. Most of the rice paddy fields had no significant annual deformation component. It might be true that the amount of soil in rice paddy fields is slowly decreasing or underground layers are slowly compacting due to groundwater usage; however, it is also possible that common subsidence signals are created by the systematic bias from the fading signal. Soil moisture change and biomass growth in the rice paddy fields are more drastic than natural herbaceous fields. However, not all subsidence signals are attributed to systematic bias. In the urbanized DID, the impact of the fading signal is negligible, and the detected signals represent true ground deformation.  Independent data sets observed from different geometries also help distinguish the true deformation signals from the systematic bias. Because the source of the systematic bias should be on the ground surface, multiple data sets would show the biased signal at the same locations. However, a deformation source related to groundwater or tectonic activity is underground and generates a horizontal component, which shifts the peak location and distorts the shape of the deforming region when viewed from an oblique LOS (Fuhrmann and Garthwaite 2019). Therefore, LOS deformations observed from ascending and descending orbits have different peak locations (Fig. 16c, d). Observations from multiple geometries are useful for both decomposition to the horizontal and vertical components and for distinguishing the true deformation with the source at depth from the surficial systematic bias.

Applicability to other areas worldwide
Because all the LiCSAR products and software used in this study are open and freely available, anyone can perform the same processing and obtain the same results. LiCSBAS automatically downloads the LiCSAR products, clips the area of interest, and performs the time series analysis with predefined parameters.
LiCSBAS processing after clipping in step 0-5 can be done quickly if the area of interest is narrow, while the downloading time of the LiCSAR products in the step 0-1 highly depends on the user's internet speed . For example, the Nara area (ID: 291) has a dimension with 36 × 20 km (324 × 216 pixel), and 81 and 133 epochs in the ascending and descending orbits, respectively. A low-cost small single-board computer Raspberry Pi 4 (Quad core Cortex-A72 1.5 GHz, 4 GB Random Access Memory) was used for a processing test. The processing time was~10 and~20 min for data conversion (step 0-2) and~20 and~35 min for clipping (step 0-5) to the time series analysis (step 1), for the ascending and descending data, respectively. The file size was 3.6 and 6.3 GB for the downloaded LiCSAR products, 9.8 and 17 GB after data conversion, and 0.3 and 0.4 GB after clipping, for the ascending and descending data, respectively. These facts indicate that users need not have high-performance computing resources or large amounts of storage, which are necessary for general InSAR processing from large low-level data (e.g., Single Look Complex). Although this study focuses on Japan, this method is easily applicable to other areas wherever sufficient LiC-SAR products are available. The LiCSAR system has produced its products mainly in high-priority areas (i.e., the Alpine-Himalayan Belt, East African Rift System, and > 1000 global volcanoes); 300,000 interferograms are available as of July 2020, and the coverage is being expanded (Lazecký et al. 2020). Even if the LiCSAR products have not been processed in the frame of interest, users can request processing from the LiCSAR frame request form on the COMET-LiCSAR web portal (COMET 2020).
The approach used in this study is also applicable to volcanoes, unless they are covered with dense vegetation. Unfortunately, most Japanese volcanoes are covered by forests, from which LiCSAR products derived from the C-band Sentinel-1 data cannot retain adequate coherence for the LiCSBAS time series analysis, even with short temporal baselines. To reveal deformation in such densely vegetated areas, L-band SAR data are preferred. It is difficult to detect the detailed deformation time series using existing L-band SAR data because the observation frequency for interferometry of the unique operating L-band satellite ALOS-2 is much lower (< 10 times/year) than that of Sentinel-1 (30-60 times/year). However, nextgeneration L-band SAR satellites such as ALOS-4 (Motohka et al. 2019) and NISAR (Jet Propulsion Laboratory 2020), which have higher frequency observation capabilities comparable to Sentinel-1, will enable detailed deformation monitoring even in forested areas.

Toward continuous operational monitoring
The deformation time series and velocities should be updated with newly acquired data by reprocessing for continuous monitoring. The LiCSAR products in prioritized regions are continuously updated (Lazecký et al. 2020). In the current LiCSBAS processing chain, all steps in data preparation (step 0) process only new data (i.e., do not reprocess the existing data); thus, new data can be efficiently added. However, all steps in the time series analysis (step 1) must be reprocessed to update the results. If the area of interest is small, the reprocessing is relatively quick (see the Section 4.2). Nevertheless, reprocessing all the time series analysis steps with every new acquisition is not efficient, because the existing results are hardly changed. A more efficient approach to updating the results (e.g., processing only a small part related to new data) should be implemented. Additionally, the preferred update frequency depends on the needs of the end users. The update strategy should be determined considering the objective of the monitoring activity.
In this study, the deformation signals are identified by visual inspection of the velocities. Small, nonlinear deformation, especially that occurring near the edges of the observation period, might be missed. Introducing automatic change detection methods (e.g., Gaddes et al. 2019) could reduce missing true deformation signals and manual operations.
Some areas are masked due to unwrapping errors included in the LiCSAR products. If the number of the unwrapping errors is small, the errors can be corrected using the redundant network of the SB interferograms (Yunjun et al. 2019) and the masked deformation data may be usable. Although currently LiCSBAS does not support the unwrapping error correction, it would be worth experimenting and implementing.

Conclusions
This study performed an InSAR time series analysis using the LiCSAR products and LiCSBAS for 73 major urban areas in Japan and revealed the spatially and temporally detailed deformation time series and velocities during 2014-2020. All LiCSBAS processing was automatically performed using predefined parameters. Many deformation signals including linear/nonlinear subsidence/uplift and annual fluctuations were comprehensively detected with high accuracy. Additionally, the detailed postseismic response to the 2016 Kumamoto earthquake was detected using the approach developed in this study. All deformation data obtained in this study, including those not mentioned above, are available on an open repository (see "Availability of data and material") and are expected to be used for further research or interpretation. This nationwide monitoring approach using the LiCSAR products and LiCSBAS is easy to implement and applicable to other areas globally.
Additional file 1: Figure S1. Map of the clipped urban areas and topography of Japan. Figure S2.    Figure S12. (a) Vertical deformation velocities overlaid with the green rice paddy fields in Nara (ID: 291). The black lines denote the outlines of the DID. The interval of the red (uplift) and blue (subsidence) contour lines is 5 mm/yr without the 0 mm/yr lines. The black square denotes the reference point. (b) Topography. (c,d) LOS deformation velocities from the 010A_05549_121311 and 017D_05514_131312, respectively. Figure  S13. (a) Vertical deformation velocities overlaid with the green rice paddy fields in Osaka (ID: 271). The black lines denote the outlines of the DID. The interval of the red (uplift) and blue (subsidence) contour lines is 5 mm/yr without the 0 mm/yr lines. The black square denotes the reference point. (b) Topography. (c,d) LOS deformation velocities from the 010A_05549_121311 and 017D_05514_131312, respectively. Figure   S14. (a) Vertical deformation velocities overlaid with the green rice paddy fields in Kumamoto (ID: 431). The black lines denote the outlines of the DID. The interval of the red (uplift) and blue (subsidence) contour lines is 5 mm/yr without the 0 mm/yr lines. The black square denotes the reference point. (b) A coseismic ALOS-2 SAR interferogram between March 7 and April 18, 2016, available on GSIMaps (https://maps.gsi.go.jp/). (c) Topography. (d-f) LOS deformation velocities from the 054A_05688_051313, 156A_05639_121313, and 163D_05736_131313, respectively. Table S1. LiCSBAS parameter settings used in this study. The name of the parameters are explained in the LiCSBAS script (https:// github.com/yumorishita/LiCSBAS)