Spatiotemporal behavior of a large-scale landslide at Mt. Onnebetsu-dake, Japan, detected by three L-band SAR satellites

We applied differential InSAR analysis to the Shiretoko Peninsula, northeastern Hokkaido, Japan. All the interferograms of long temporal baseline (~ 3 years) processed from SAR data of three L-band satellites (JERS-1, ALOS, ALOS-2) commonly indicate remarkable phase changes due to the landslide movement at the southeastern flank of Mt. Onnebetsu-dake, a Quaternary stratovolcano. The area of interferometric phase change matches to known landslide morphologies. Judging from the timing of the SAR image acquisitions, this landslide has been moving at least from 1993 to the present. Successive interferograms of 1-year temporal baseline indicate the temporal fluctuation of the landslide velocity. Especially for the descending interferograms, the positive line-of-sight (LOS) length change, which indicates large subsidence relative to the horizontal movement, is observed in the upslope section of the landslide during 1993–1998, while the negative LOS change is observed in the middle and the downslope section after 2007 indicating less subsidence. The landslide activity culminates from 2014 to 2017: the eastward and the vertical displacement rates reach ~ 6 and ~ 2 cm/yr, respectively. Utilizing high spatial resolution of ALOS and ALOS-2 data, we investigated velocity distribution inside the landslide. During 2007–2010, the eastward component of surface displacement increases toward the east, implying that the landslide extends toward the east. During 2014–2017, the vertical displacement profile exhibits spatially periodic uplift and subsidence consistent with surface gradient, which indicates the ongoing deformation driven by gravitational force. Heavy rainfall associated with three typhoons in August 2016 might have brought about an increase in the landslide velocity, possibly due to elevated pore-fluid pressure within and/or at the base of the landslide material. Also, annual rainfall would be an important factor that prescribes the landslide velocity averaged over 3 years.


Introduction
In monsoonal Asia, landslides due to heavy rainfall frequently caused severe damages. To mitigate landslide damages, it is essentially important to map the landslide-prone areas in advance. The Japan Seismic Hazard Information Station (J-SHIS) provided by the National Research Institute for Earthquake Science and Disaster Resilience (NIED) extensively integrated such potential landslide information in Japan (National Research Institute for Earth Science and Disaster Prevention 2020). To step forward, two pieces of additional information are required: (1) current velocity of the landslide movement and its temporal change; (2) internal deformation of the landslide mass. Space-borne Synthetic Aperture Radar (SAR) is an ideal tool to investigate the above-mentioned points. By taking phase difference between temporally separated two SAR images, called interferometric SAR (InSAR), we can obtain surface velocity maps with a very high spatial resolution (e.g., Massonnet et al. 1993). More specifically, what we can measure with InSAR is called as line-of-sight (LOS) displacements: the inner product of Takada and Motono Earth, Planets and Space (2020) 72:131 surface displacement vectors and unit vectors of radar incidence direction. So far, many previous studies have applied the InSAR technique to landslide movements using C-band satellites (e.g., Hilley et al. 2004;Colesanti and Wasowski, 2006). The landslide in densely vegetated areas have also been detected by InSAR analysis of L-band SAR data acquired by ALOS (e.g., Zhao et al. 2012;Handwerger et al. 2013). Recently, the high spatial resolution of ALOS-2, the latest L-band SAR satellite, has been utilized to detect landslide movements (e.g., Sato and Une, 2016;Nishiguchi et al. 2017;Fujiwara et al. 2017). Furthermore, Hu et al. (2019) applied the InSAR technique using both ALOS and ALOS-2 data to a slowmoving landslide in California, USA. In this study, we applied InSAR analysis using three L-band satellites to investigate a large-scale slow-moving landslide in the Shiretoko Peninsula, Japan.
The Shiretoko Peninsula, one of the UNESCO world heritages, locates at the northeastern part of Hokkaido ( Fig. 1), where many landslide geomorphologies have been identified by aerial photographs (Ito, 1994(Ito, , 1996. One of the largest landslide morphologies in this region has been found at the southeastern flank of the Mt. Onnebetsu-dake, a Quaternary stratovolcano (height above sea level: 1330 m). The current volcanic activity of Mt. Onnebestu-dake is very low, and its edifice is highly dissected. Figure 2 shows the topographic map of the proximity of Mt. Onnebetsu-dake, in which a gentle slope extending to the southeast is the landslide morphology studied in this paper. This landslide area is bounded by two tributaries of the Shunkarikotan River as shown in Fig. 2. The length and width of this landslide are 4.5 and 1.6 km, respectively, which associates the main scarp with a relative height of 250 m and several lateral scarps especially along its southwest margin (Ito, 1994). Within the landslide mass, we can observe minor scarps, pressure ridges, and depressions (small lakes) formed by the landslide activities (Ito, 1994). The current morphology was formed by large-scale rock-slump-type movements and subsequent rock-slides and debris flows (Ito, 1994). The current velocity of this landslide is poorly known, mainly due to many obstacles to the field observations. The Shiretoko Peninsula is covered by primeval forest (no trail in the Mt. Onnebetsu-dake) and well known as a dense population of brown bears. Also, the Shiretoko Peninsula is characterized by heavy snow-fall in winter and spring, which makes it difficult to conduct ground-based geodetic observations. To overcome such difficulties, in this study, we use a series of InSAR images acquired by L-band SAR satellites taking advantage of high interferometric coherence. We use three successive L-band satellites, JERS-1 (1992-1998), ALOS (2006-2011, and ALOS-2 (2014 to present) to investigate the spatiotemporal behavior of the landslide at the southeastern flank of Mt. Onnebetsu-dake. We further investigate its internal velocity field by making the most of ALOS and ALOS-2 high spatial resolution. Finally, we discuss the relationship between landslide activity and the amount of precipitation, especially under abnormal weather conditions.

Method
We processed interferograms of SAR data acquired by JERS-1, ALOS, and ALOS-2 with the GAMMA software suite (Wegmüller and Werner, 1997). To remove topographic fringes, we used a 10-m mesh digital elevation height model (DEHM) released from the Geospatial Information Authority of Japan (GSI). The original data in this area have been acquired on 1976, and improved on 1995 and 2017, which is digitized and provided by GSI. We adopted the minimum cost flow method (Costantini, 1998) for phase unwrapping. The interferometric pairs are listed in Table 1. The SAR data coverages are illustrated in Fig. 1, where the paths and frames ('lows' for JERS-1) of each satellite acquisition are indicated as well.
All the ALOS-2 data used in this study were acquired in right-looking mode.
It should be noted that the SAR images acquired in winter and spring cannot be used for interferometry because the snow cover causes coherence loss. Considering such limitation of useful InSAR images, we decided not to apply any advanced time-series analysis at this stage (e.g., Feretti et al. 2001;Berardino et al. 2002). Also, if the landslide movement changes with time drastically, the minimum norm solutions may smooth out the temporal change in velocities. Thus, in this study, we opt to compare a series of interferograms with each other to extract the temporal behavior of the landslide.
We aim at flattening the systematic noise in the proximity of the landslide. For this purpose, each interferogram was corrected by a two-step approach. First, we modeled the long-wavelength phase trend by a quadratic function, and removed it. The removed components are

Table 1 Interferometric pairs used in this study
Asc and Dsc are the abbreviations of ascending and descending orbit, respectively, Bp means perpendicular baseline, Inc means incidence angle at scene center, SS beam mode of master and slave data are FBS (fine beam single polarization), DD beam mode of master and slave is FBD (fine bean dual polarization), SD combination between FBS and FBD, U beam mode of master and slave is ultra-fine. Beam mode of JERS-1 SAR is unique typically due to orbit estimation error, ionospheric and/ or tropospheric disturbances with length-scale much larger than that of landslides. Second, we cropped out the proximal area of Mt. Onnebestu-dake. In the cropped area, we still find a remnant of the long wave-length noise as well as a height dependence of the interferometric phase due to an elevation-correlated tropospheric noise. For all the interferograms, we modeled such systematic noise in the cropped area by a linear function of height along with linear function in space (ramp), and removed it. In the following part of this paper, the cropped areas are shown as blow-up figures.

Overall activities
First, to grasp overall landslide activities, we focus on the long temporal baseline pairs by which the signal-tonoise ratio can be enhanced. Figure 3 indicates the interferograms taken from descending orbit, spanning about 3 years (2 years for Fig. 3g and h). At the southeast flank of the Mt. Onnebestsu-dake, we see a remarkable line-ofsight (LOS) length change in an area elongated toward the southeast in the interferograms of JERS-1 (Fig. 3a), ALOS ( Fig. 3c) and ALOS-2 ( Fig.3e and g). In this paper, we define positive LOS length changes as surface movements away from the satellite. We set the reference point (black dot in Fig. 3) within 1 km distance from the landslide. The blow-up figures of the landslide area ( Fig. 3b, d, f, h) indicate that the LOS displacement inside the landslide is not uniform, which will be discussed later in detail. Using ALOS and ALOS-2 images, we also processed ascending interferograms spanning about 3 years ( Fig. 4, 2 years for Fig.4e and f ). All the ascending interferograms clearly display the positive LOS displacement due to the landslide movement at the southeastern flank of the Mt. Onnebetsu-dake. Thus, averaging over a long (2-3 year) time-scale, it is highly likely that the landslide has been moving persistently from 1995 (actually 1993 as shown in the next section) to the present at a rate slower than 7.3 cm/yr (Fig. 5). Such temporal behavior is called as "slow-moving landslide" (e.g., Hilley et al. 2004;Handwerger et al. 2013;Hu et al. 2019). Incidentally, if the direction of surface movement coincides to that of steepest gradient (roughly 9 degrees down to N147E in the middle stream area), the LOS velocity of 7.3 cm/yr corresponds to the surface velocity of 20.0 cm/yr. Figure 6 indicates the landslide morphologies provided by J-SHIS (https ://www.j-shis.bosai .go.jp/en/) draped over the ALOS-2 interferogram indicated in Fig. 3f (2014/09/15-2017/08/28). It should be noted that the reference point does not locate within any landslide morphology. From Fig. 6, we understand that (1) the area of LOS length-change coincides with mapped landslide morphologies (Fig. 6b), and (2) many other landslide morphologies in and around Mt. Onnebetsu-dake remain stable (Fig. 6a). The mapped landslide morphologies indicate that the landslide area is expected to flow down in the southeast direction as a whole (Fig. 2), which is consistent with the overall trend of Figs. 3 and 4: negative and positive LOS changes from descending and ascending orbits, respectively.

Internal velocity structure
In this section, we focus on the internal deformation pattern of the landslide by decomposing the surface displacement vectors into the quasi-east and the quasivertical components lying in the plane including 2 LOS vectors for ascending and descending interferograms (Fujiwara et al. 2000). Hereafter, we call the quasi-east and the quasi-vertical simply as east (or eastward) and vertical, because the former is almost the same as the east and the latter deviates less than 8 degrees from vertical. We selected ALOS and ALOS-2 interferograms with more than 2-year temporal baseline to enhance the signal-to-noise ratio. Figure 7 indicates the eastward and the vertical components of the surface velocities for 2007-2010, 2014-2017, and 2017-2019. Overall, the eastward component is much larger than vertical component; for 2014-2017, the former reaches 6 cm/yr while the latter is about − 2 cm/yr at most. The vertical velocities in Fig. 7 display remarkable temporal change in the spatial patterns as described below.
During 2007-2010 (Fig. 7b), the uplift rate increases from west to east, and culminates along the eastern edge bounded by a stream flowing down to the southeast. Combined with the eastward increase in the eastward component (Fig. 7a), the landslide mass is expanding toward the east and would be burying the stream.
During 2014-2017, the spatially periodic movement can be identified in the vertical component (Fig. 7d). Figure 8 indicates the profile of the surface velocities as well as a surface gradient along the section line AA' in Fig. 7d. Along the section AA' , the vertical velocity and the surface gradient changes in harmony (Fig. 8a): rapid subsidence where relatively steep surface slope, slow subsidence or slow uplift where very shallow to upslope, which implies that the internal deformation of the landslide is essentially driven by gravitational force. In Fig. 7d, each area of relative subsidence trends in the transverse direction (i.e., NE-SW to ENE-WSW), that is, perpendicular to the entire landslide morphology, which implies the southeast movements of landslide mass as a whole. Figure 8b indicates the eastward components of surface velocities along the section AA' , which dominantly shows eastward movement with a long wave-length undulation ranging ~ 2 to ~ 4.5 cm/yr. The slow eastward movements around 145.03 E and 145.046 E (Fig. 8b) are in good accordance with the fast uplifts there (Fig. 8a). For the latter (145.046 E), the surface gradient is the largest (slightly upslope) implying that the downslope mass transportation was obstructed by the upslope morphology, and the transported material was clogged up to go upward.
From 2017 to 2019, Fig. 7f does not display distinct vertical motion in comparison to the background noise level over the entire region. The eastward component of surface velocity indicates a gradual increase toward the internal part, but internal deformation is very slow during this period.

Temporal change of landslide activities
In this section, we focus on the temporal change of the landslide activities. Considering no ascending data acquisition by JERS-1, we begin with comparing descending interferograms that cover a longer period than ascending ones. Although the landslide movements are commonly detected by three satellites (Fig. 3), the LOS velocities are not constant in time. The maximum LOS velocities in the processed interferograms are compiled in Fig. 5. The JERS-1 interferograms indicate length increase in the LOS direction at rates from 3.1 to 4.1 cm/yr (Figs. 3b, 5, 9a), while ALOS and ALOS-2 descending interferograms show clear length shortening (Fig. 3d, f, h). The interpretation of this flip in LOS sign will be described in a later section. Figure 5 indicates that the landslide velocity is relatively low (− 3.0 cm/yr) from 2007 to 2010, then is raised to − 4.7 cm/yr from 2014 to 2017, and declines to − 0.7 cm/yr from 2017 to 2019. The same trend is found in the ascending interferograms (Figs. 4 and 5); the velocity is low (1.6 cm/yr) during 2007-2010, then is increased to 7.3 cm/yr during 2014-2017, and declines to 1.8 cm/yr from 2017 to 2019.
To investigate the temporal change much in detail, for both of the ascending and descending orbits, we split the ALOS and ALOS-2 interferograms spanning 3 years (2 years after 2017) into successive interferograms over 1 year, respectively (Figs. 9 and 10). We included a JERS-1 interferogram spanning 1993-1995 in Fig. 9, for reference. The maximum LOS velocities are compiled in Fig. 5 as well, where it should be noted that the maximum velocity of each 3-year interferogram does not agree with an average of maximum velocities over 3 interferograms of 1-year temporal baseline because the fastest moving point is changing with time. The fastest moving point strongly reflects the local deformation especially at the toe and along the stream running northeast of the landslide area (Figs. 3, 4, 9, 10). To grasp the entire activities, we calculated the LOS velocities averaged over the landslide area as defined in Fig. 3b, by which the difference between the long-term and the short-term maximum (minimum) velocities are mitigated (Fig. 5). For ascending orbit, we could not make 1-  Fig. 9d), which, however, does not indicate remarkable reactivation of whole landslide area as represented by the averaged rates of − 1.3 cm/yr (Fig. 5). From 2014 to 2017, the overall velocities get higher than 2007-2010 for both ascending and descending orbit that reach 15.9 cm/ yr and − 14.9 cm/yr, respectively (Figs. 5,9,10). During this period, the LOS velocity seems to fluctuate. In concrete, the descending image exhibits sudden speed-up during 2015/09/28-2016/08/29 (Fig. 9f ). On the other hand, the ascending image displays relatively low velocity in 2015/06/29-2016/06/27 (Fig. 10e, but higher than 2007-2010). Such differences between the ascending and the descending images may be explained by (1) temporal change in relative amount of horizontal and vertical component (i.e., direction of surface motion), or by (2) an abnormal increase in precipitation in July and/or August of 2016 which is reflected only on the descending image.

Discussion
In this section, first, we focus on the cause of positive LOS displacement in the descending interferograms taken by JERS-1 during 1993-1998 (Figs. 3b and 9a). In principle, LOS displacements are strongly dependent on radar incidence angles. For example, smaller incidence angles make interferometric phase changes more sensitive to vertical components of surface displacements. However, the incidence angle of JERS-1 SAR data is almost the same as that of ALOS and ALOS-2 (Table 1)   which, therefore, does not account for the change in LOS displacement from positive (1993-1998) to negative (2007-2010 and 2014-2019). The positive LOS displacements in the JERS-1 interferograms are observed in the upslope section of the landslide (Figs. 3b and 9a), while the negative LOS displacements in ALOS and ALOS-2 are mainly observed in the middle to the downslope section of the landslide (e.g., Fig. 3d and f ). Following the mass conservation law, it is widely observed that the upstream area of active landslides subsides due to material depletion, and that the provided material flows down along gentle slope of the downstream area (e.g., Varnes, 1978;Hu et al. 2018). Let us consider the geometric relationship between the LOS direction and topographic features under such situation (Fig. 11). The upslope section is characterized by steep scarps that enhance the surface subsidence rather than horizontal motion. The downslope section, on the other hand, is characterized by shallow deposit and gentle slopes that prompt the horizontal motion rather than surface subsidence. Thus, the positive and negative LOS displacements would be preferable at the upslope and the downslope sections, respectively, as illustrated in Fig. 11. Not to mention, however, the spatiotemporal change of the LOS displacement that we observed cannot be explained solely by above-mentioned geometrical relationship. We should keep in mind that complex hydraulic systems can cause such enigmatic behavior of slow-moving landslides (e.g., Handwerger et al. 2013;Hu et al. 2019). For example, heterogeneously high pore-fluid pressure at a limited area may reduce the local surface slope angle supported by internal and basal friction (e.g., Davis et al. 1983). Furthermore, an abnormal increase in pore-fluid pressure may reduce the frictional strength (e.g., Terzaghi 1950), which may account for such remarkable temporal change. Thus, to better understand such spatiotemporal behavior, we have to constrain the hydraulic properties by field observations as well as geodetic surveys.
Next, we discuss the effect of height error of the DEHM on the periodic uplift shown in Figs. 7d and 8a. The apparent LOS displacement u caused by the DEHM error h is written as:  Fig. 8 where B ⊥ , R , and θ denote the perpendicular baseline, distance between the satellite and the ground, and look angle, respectively (Hanssen, 2001). This equation states that the effect of DEHM error is proportional to the baseline length. In Fig. 7b, we cannot find clear periodic uplift unlike Fig. 7d, while the baseline length of the former is much larger (> 1000 m) than the latter (Table 1). Therefore, the periodic uplift pattern in Figs. 7d and 8a should not be caused by the DEHM error. We further verify this conclusion quantitatively. We set B ⊥ , R , and θ as 134 m, 8.02 × 10 5 m, and 35.4 degrees, respectively, considering the interferogram 2014/09/15-2017/08/28, which is used for calculation of the uplift rate in Fig. 7d. Assuming that the periodic uplift rate in Figs. 7d and 8a were caused by the height error of the DEHM, we projected amplitude of the spatial change in the uplift rate (0.5-1.5 cm/yr as indicated in Fig. 8a) to the LOS direction, and multiplied temporal baseline of the interferogram, by which we obtain u as 1.2-3.6 cm. Then, following Eq. (1), we can estimate h as 41.6-124.8 m, which is far larger than the DEHM error officially announced by GSI (5 m). Thus, the periodic uplift rate in Figs. 7d and 8a is not an artifact due to the DEHM error. To confirm advantages of the DEHM, we made a map of height difference between the DEHM and SRTM (Fig. 12).
Here, we removed the geoid height of 30 m around the Shiretoko Peninsular (Kodama et al. 2014) from the DEHM. Figure 12 shows small height difference in the landslide area, where surface slope is gentle. The height difference is large in the area of steep surface slope especially around the main scarp. Figure 12b shows the height difference along the section line in Figs. 7d and 12a, which indicates periodic undulation with amplitude mostly ranging 2-15 m, which is much smaller than h that we estimated above. In conclusion, the spatial resolution of SRTM and the DEHM provided by GSI is 90 and 10 m, respectively, and the latter is advantageous for removing topographic fringes associated with the small-scale surface topographies in the study area. Next, we further discuss on the source of errors included in the interferograms that we processed. As we mentioned before, we have removed (1) long wave-length phase trend assuming a quadratic function in space, and (2) height-dependent phase assuming a linear function of altitude as well as linear phase trend in space. In spite of these efforts, we still find phase changes that seem not related to the landslide movements especially near the main ridge running along the highest part of Shiretoko Peninsular (Fig. 2). In Fig. 9f, for example, we find negative LOS displacements at the uppermost part of the landslide and to the northeast of the landslide (around 44.02 N, 145.05 E). Such noise is remarkable only at high altitude area (> 1000 m), which does not strongly affect the LOS change detected in the middle to downstream area with relatively low elevation (Figs. 2 and 8). In this study, we did not apply advanced time-series analyses (e.g., Feretti et al. 2001;Berardino et al. 2002)   Finally, we examine the effects of precipitation on landslide activities. Figure 13a indicates monthly precipitation at Utoro (Fig. 1), the nearest meteorological station ~ 9 km NNW of the Mt. Onnebestu-dake, from January 1993 to December 2019 (data provided by Japan Meteorological Agency). The data sampling interval is 1 month, and the gauge measures the precipitation as snow and rainwater combined. We also plotted the cumulative precipitation which is reset to zero every October 1 (Fig. 13b). An abnormally large amount of rainfall (~ 800 mm per month) was recorded during August 2016, which is caused by three typhoons passing by the Shiretoko Peninsula (e.g., Kitano et al. 2017) (Additional file 1. Fig. S1). The interferograms over August 2016 display fast movements of the entire landslide region (Figs. 9f and 10f 11 Geometric relationship between radar incidence angle and surface topography of rotational landslides ( modified from Cruden and Varnes, 1996). Red and green arrows indicate surface movement of upper and lower portion of landslide, respectively. Blue line indicates radar line-of-sight direction. Incidence angle is set to be 39 degrees, as is the case of JERS-1 (Table 1) would increase pore-fluid pressure of the landslide material and reduce the frictional strength, which may increase the surface velocity of the landslide at the southeast flank of Mt. Onnebetsu-dake, as reported by previous studies (e.g., Iverson and Major, 1987;Handwerger et al. 2013;Hu et al. 2019). Unfortunately, the ALOS-2 data acquisition is not frequent enough to resolve the effect of typhoon passages in a few days (Additional file 1. Fig. S1). This problem should be investigated with the aid of more frequent Sentinel-1 data acquisition in the future. Next, we compare annual precipitation and surface velocity detected by long temporal baseline interferograms (Figs. 3 and 4). As indicated in Figs. 13b and 14, the annual precipitation at Utoro station seems slightly increasing with time. Total rainfall during 1995Total rainfall during -1998Total rainfall during , 2007Total rainfall during -2010Total rainfall during , and 2014Total rainfall during -2017.5 mm, respectively, which would account for faster landslide velocity during 2014-2017 (Figs. 3f and 4d) relative to that during 2007-2010 (Figs. 3d and4b). But the difference in the deformation pattern between 1995-1998 and 2007-2010 ( Fig. 3b and d) cannot be explained from the total amount of precipitation.

Conclusion
We investigated the spatiotemporal behavior of a large-scale landslide at the southeastern flank of Mt. Onnebestu-dake, in the Shiretoko Peninsula, from 1993 to the present using L-band SAR images acquired by three satellites. Many coherent interferograms over densely vegetated areas commonly indicate persistent movement of the landslide during this period, although SAR data have not been acquired during 1999-2005 and 2012-2013. We found pronounced temporal changes in the landslide velocity. In concrete, the positive LOS displacement is dominant in the upslope section during 1993-1998, and the negative LOS displacement is dominant in the middle and downslope section from 2007 to 2019. The landslide velocity culminates during 2014-2017 where the eastward component reaches 6 cm/yr while the vertical subsidence rate is 2 cm/yr at most. The internal deformation pattern exhibits periodic uplift and/or subsidence consistent with topographic gradient, which can be observed only during 2014-2017. Heavy rainfall caused by three typhoons in August 2016 might have brought about an increase in the landslide velocity. Also, the annual amount of precipitation would be an important factor that prescribes the landslide velocity averaged over 3 years. To understand the physical mechanisms driving such remarkable spatiotemporal variations, further acquisition of L-band SAR images with shorter time intervals and ground-based field surveys with the aid of the information reported in this paper would be necessary.