Occurrences of Deep-Seated Creeping Landslides in Accordance with Hydrological Water Storage in Catchments

Rainfall is one of the most important triggers of both shallow debris flows and deep-seated landslides. The triggering mechanism involves the process of water infiltration into the failure zone. For deep-seated landslides, the deeper and more extensive failure surfaces delay the effect of the process and thus delay landslide initiations. The delay is difficult to assess, especially if the sites only have scarce or insufficient monitoring data. Under these circumstances, we illustrate that the occurrences of landslides can be estimated by their correlations with the phenomenological water storage index (WSI) of a given catchment. In the present study, a total of five deep-seated landslides in TienChih (4) and SiangYang (1) are investigated. The displacements of the landslides were recorded by global positioning system (GPS) and the WSI was modelled using the tank model. The result demonstrates that the WSI correlates closer in time to the landslide motion than the rainfall, and the WSI thresholds for the landslides are inferred. Thus, this technique can be applied as an associated method to evaluate landslide initiation.


INTRODUCTION
Recently, in Taiwan, with continuous global positioning system (GPS) network data, eight deepseated creeping landslide sites were identified and their causes were associated with hydrological loading, Hsu et al. (2014). The evidence was based on relating the creeping landslide events to the rainfall intensity (RI) or accumulated precipitation. In this paper, we further inspect the landslide events with a simple phenomenological quantitative hydrological model. The model, known as the tank model of Sugawara and Funiyuki (1956) and Sugawara (1961Sugawara ( , 1995, can provide a simple idealized sketch for a regional water cycle: including precipitation, infiltration and runoff discharge. Using the model, we attempt to identify model thresholds for landslide occurrences and observe the sequential relationship between the landslides and regional water cycle. Because of the wide extent area and deep failure surface of deep-seated landslides, the rainfall triggering mechanism involve delicate infiltration processes with site geological conditions. Detail investigations of the infiltrated water flows, however, can involve sophisticated techniques. For example, Cappa et al. (2004) traced sulfate, nitrate concentration and spring waters on the slope to identify the characteristics of precipitation infiltration and concluded that landslide motions correlate to precipitation. Padilla et al. (2014) measured electrical conductivity, oxygen and hydrogen isotopic concentration of the groundwater and correlated the recharged groundwater to rainfall. The infiltration is found to be associated with high hydraulic conductivity area in the bedrock caused by the gravitational deformation in the investigated hillslope. Lv et al. (2019) applied tracer test, injection test and electrical resistivity tomography to resolve the underground water flow in their study area. The landslide is found mainly controlled by a vertical preferential flow pathway, one of the two major types of underground water flow patterns of the site. Zêzere et al. (2005) analyzed a series of long-term monitoring data since 1950s on the occurrences of rainfall induced shallow and deep landslides. The deep-seated landslides correlate to long period but less intensive rainfalls and it indicates that the infiltration causes the rise of groundwater table and the deep failures through reduction of shear strength of the slope materials. Nie et al. (2017) extended a tank model to include effects of snow accumulation, snowmelt and infiltration. These effects resulted in time lags and the model improved the groundwater level predictions in deep-seated landslides.
To adapt a simple tank model in the present study is because the landslides were in remote areas where detailed geological surveys and field observation data required to constrain sophisticated models are scarce or insufficient. In addition, because the lack of field data is common in the mountainous areas of Taiwan Island, a simple model can be a valuable tool that can be incorporated into hazard mitigation systems if its applicability is validated. For this purpose, we selected two landslide sites: the TienChih and SiangYang areas. Provided that a minimum set of hydrological data is available, the selection criteria are described as follows: the sites are in the catchments of streams of the lowest-possible order (furthest upstream); the catchments are as small as possible; and the landslide is of the translational type, according to the landslide classification in Varnes (1978). The concept of applications of the tank model was first tested feasible in a preliminary study of the TienChih area, Lin et al. (2017), and is extended in the present paper.
Both sites are on the southern transverse roadway. The roadway was built in the 1970s to connect the southwest (Kaohsiung City) and southeast (Taitung County) parts of the island. The total length of the roadway is approximately 209 km and most of the mountainous section of the roadway is in Yushan National Park. The highest elevation of the roadway is 2,722 m at approximately 147.5 km along the roadway. On the immediate west side (Kaohsiung side) of the peak point, Yakou, the terrain is in the catchment of the Laonung River and, on the east (Taitung side), the terrain is the Beinan River catchment. The Laonung River flows to the southwest, and the Beinan River flows to the southeast. The upstream-most flow-rate-monitored catchments of the studied regions are plotted in Figure 1A. The river flow rate observation stations for the two catchments are A-Chi-Ba Bridge and Hsin-Wu-Lu. Their monitoring catchment areas are 411 and 628 km 2 , respectively. In the catchments, there are two and four rainfall monitoring stations, as marked in the figure. However, for simplicity, only the data recorded by the rainfall station closest to the landslide site are used to represent precipitation in each catchment in the current analysis. The adapted rainfall stations are marked in Figure 2.
The geological settings are shown in Figure 1B. The roadway cuts through two primary faults, the Singuan and Keiku Faults, in the south Central Mountain Range, adjoined with abundant secondary named or unnamed faults, according to Lee and Wang (1985); Lai et al. (1989) and Lin (2012). These faults are mainly of the thrust type and are arranged in the NNW-SSE direction. The geological layer at both sites is the Pilu Formation. The formation extends from NNE to SSW along the ridges of the Central Mountains, as marked in the light blue colour in Figure 1B. The formation mainly consists of Eocene rocks: slate and phyllite. Moving westward from the TienChih site, the surface layer changes to the Miocene Lushan Formation (argillite, slate, phyllite), and towards the east from the SiangYang site, it becomes Late Paleozoic and Mesozoic Tananao schists (black/green/siliceous schist). These rocks all have rich joint, cleavage and fold structures which indicates an active orogenous mechanism in the region. The common filling material in the Pilu Formation is quartz-rich calcite Lai et al. (1989), which generally has high water permeability.
The locations of the TienChih and SiangYang creeping landslide sites are marked in Figure 1 and they are enlarged in Figure 2. Both sites were designed as resting and sightseeing areas on the roadway such that facilities, including a few monitoring instruments, were constructed. Therefore, creeping landslides were found by the installed instruments. The displacement measurements and field observations will be reported in the next section.

GPS NETWORK AND IDENTIFICATION OF CREEPING LANDSLIDES
GPSs have been gaining widespread use in real-time monitoring of surface displacements. With proper postprocessing procedures, continuous GPS (CGPS) data can achieve an accuracy of several millimetre. The data presented in this section were acquired from the island-wide GPS network constructed by the Institute of Earth Sciences, Academia Sinica and Central Weather Bureau since 1989, Yu et al. (1997Yu et al. ( , 2003; Chen et al. (2015). At the TienChih site, a CGPS station that is codenamed TENC, was set up, and recording started in 2002. At the SiangYang site, a similar station known as SIANY, started operating in 2009.
The postprocessing procedures of the collected GPS data were performed in the framework of GAMIT/GLOBK software, version 10.3, Herring et al. (2002). Double-difference phase observables were used and the influences of tropospheric and ionospheric effects were corrected. To obtain an accurate and consistent regional deformation pattern in Taiwan, a total data of 362 Taiwan and 17 IGS Asia-Pacific CGPS sites were used. The least squares linear fit was used to identify the tectonic linear FIGURE 1 | (A) The upstream catchments and the section of the south transverse highway in the study area (thick solid black line). The highlighted northwest area is the 411 km 2 TienChih catchment and the southeast area is the 628 km 2 SiangYang catchment. The landslide sites, the TienChih and SiangYang areas, are marked by thick hollow red boxes. The coordinate system is TWD97. (B) The regional fault system and lithology. MI, Miocene argillite, slate, phyllite; Ep, Eocene slate, phyllite.
Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 743669 3 station velocity, annual and semiannual periodic motions, postseismic relaxation, and offsets caused by coseismic jumps and instrument changes in the station position time series, Altamimi et al. (2007). Then, the statistical outliers, linear velocities, postseismic relaxation, and instrumental offsets were removed.
The residual GPS position time series of both stations are plotted in Figure 3. They exhibited several clearly identifiable intermittent displacements, of which TENC recorded four because of its longer period in operation. These intermittent displacements were correlated with extreme precipitation events:  Figure 4 which correspond exactly to the descending direction of the slopes, i.e., translational landslide. Because of poor climate conditions (typhoons, heavy clouds and rainfall) when landslides are in motion, GPS stations may not obtain data, site accessibility is lost or instrumental damage occurs. The breaks in data impede hydrological analysis, except for the unique recordings taken in the Morakot event, cf. Figure 10.
The cumulative rainfall of Morakot typhoon was reached 2,700 mm in 5 days. The amount of precipitation was almost equal to the average annual rainfall on the island. This typhoon resulted in numerous landslides, excessive floods and severe fatalities in the southern part of Taiwan and after the strike, a series of light detection and ranging (LiDAR) surveys was performed over that area in 2010 by the Central Geological Survey. The LiDAR technique produces high-resolution digital elevation models (DEMs) and is superior in terms of highprecision and vegetation-penetrating capacities. The derived DEM of the LiDAR survey had a grid size with a 1 × 1 m resolution based on the average point cloud density of 0.5 point/m 2 . For the two landslide sites, the DEMs are shown in Figure 4.
With high-resolution DEMs, surface cleavages, scarps and outcrops can be visually identified in the indoor preparation phase prior to field geological surveys. As shown in Figure 4, typical sackungen structures are found along ridge crests in massive competent rocks, running parallel to the slope contours and to the strikes of main rock discontinuities (foliation). Surface failure usually develops steep into the mountainside by gravitational stresses; therefore, scarps are created. For example, Figure 5A presents the slope profile of a shallow landslide area after the Morakot typhoon (the central scarp near the rainfall gauge marked in Figure 2A). A typical cleavage structure (fault gouge zone) in the area and underground water seepage are shown. Figure 5B illustrates the crack openings  Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 743669 5 found in the main building site of the roadway rest area. The movements, although slow (long-term gravity-driven creep), can continue for long periods in principle, producing significant cumulative displacements.
In the preliminary study? Chen et al. (2015), we visually examine the correlation between the temporal variation in GPS displacement time series and rainfall data (intensity and accumulated). The displacement takes place with a time lag compared with the rainfall hydrographs. This phenomenon suggests that the motion relates to the effect of water infiltration which alters the gravitational loads and effective stresses during typhoons and heavy precipitation. Combining the LiDAR and field observations, the two creeping landslides are concluded to be of a deep-seated cleavage and groundwatercontrolled type, Chen et al. (2015). In what follows, the relationship between the displacement of the creeping landslide and the regional water cycle is investigated. The regional water cycle will be modelled using the hydrological tank model.

TANK MODEL AND PARAMETER CALIBRATION
The tank model is a simple phenomenological quantitative hydrological model that models the procedure of water flow through soil using a series of linear reservoirs, or tanks Sugawara (1961Sugawara ( , 1995. Conceptually, each tank has several orifice exits to model discharging of water. The water in each tank represents the underground water storage and the discharges through the exits are the runoff flow to the river or the interflow to the lower cascading tanks. The computation of the resultant model is easy and fast once the model parameters are obtained. As a result, the tank model has been widely applied to flood forecasting in Asian countries, such as in Japan and Taiwan Osanai et al. (2010); Chen et al. (2013).
For example, Ishihara and Kobatake (1979) developed a synthetic flood forecasting model that combined the tank model for the rainfall-runoff conversion process of the subbasins of the studied catchments and the time-areaconcentration diagram for the runoff concentration process. Lee (1993) identified that tanks provide essential continuous damping and lagging effects to model water infiltration and discharge in catchments. He also demonstrated that the tank model can be equivalently represented by a variety of linear reservoir models with different connection configurations. Yokoo et al. (2001) illustrated that to some extent the tank model parameters can be evaluated based on the geographical characteristics of the basin. Takahashi et al. (2008) used the dynamic dimensional search algorithm for parameter calibration and applied the model to forecast the groundwater table on slopes. Liu et al. (2010) applied the Kalman filter to assimilate the tank model prediction and measurements for the Singapore River catchment in Singapore. Chen et al. (2013) analysed hundreds of mudflow events and proposed a debris warning system in Taiwan, based on the soil water index (SWI), which was the total depth of the tank model. Figure 6 depicts the tank model used in the present study. The simplified model contains only two serially connected tanks and the water flow in the system describes the water infiltration and discharge in a unit area of the catchment. In the present paper, the model is separately applied to the TienChih and SiangYang catchments in Sections 4, 5. The upper tank takes rainfall intensity, RI or R in mm/hr, as the water input and its water Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 743669 6 storage is represented by water depth h 1 (t) in mm in the tank. In the tank, there are three exits (top-side, bottom-side and bottom orifices) to discharge its water storage. These exits are used to model the surface runoff discharge, seepage of the surface layer and infiltration into the deep bedrock layer. The tank is assumed to be linear, i.e., the discharges of the three exits are proportional to the storage in the tank where A denotes the area of the catchment to which the model is applied, and H (·) is the Heaviside function to ensure unidirectional water discharge and infiltration. For the upper tank, there are five constant parameters: H 1 and H 2 are the controlling threshold depths, a 1 and a 2 are the discharge coefficients of the two side exits at 1/hr, and b 1 is the infiltration coefficient for flow infiltrating into the second tank at 1/hr. The governing equation of the storage of the tank is d h 1 (t)/dt R − (q 1 + q 2 + q 12 )/A, with the initial storage h 1 (T 0 ) at the initial time T 0 . The initial storage, h 1 (T 0 ), is denoted as S 1 .
The lower tank in Figure 6 models the water flow and storage in the deep bedrock layer. The water storage of the lower tank is defined as h 2 (t) and it has two discharge exits such that their discharges are q 3 a 3 A (h 2 − H 3 )H (h 2 − H 3 ) and q b b 2 Ah 2 H (h 2 ). Similar to the upper tank counterpart, the governing equation for the lower tank is d h 2 (t)/dt (q 12 − q 3 − q b )/A, with the initial value h 2 (T 0 ) S 2 . There are three additional constant parameters defined, a 3 , b 2 and H 3 , and their implications are the same as those of the upper counterpart. To summarize, there are two model equations to describe the water storage in each tank. The side outflows, q 1 , q 2 and q 3 , contribute to observable flow discharge and the discharge through the bottom-most exit, q b , represents the discharge through an unobservable underground path. Therefore, the sum of the side exits is the volume flow rate discharged from the catchment. Rainfall in the model is the sole water input to the system.
The key triggering mechanism of landslides is the water content in slopes. During rainfall, the increase in water content increases the pore water pressure such that the effective stresses of slope materials decrease; hence, landslide motion is triggered. With the tank model, the water content in slopes is thought to be proportional to the total water storage of the two tanks, h 1 + h 2 . Hence, we validate this assumption by inspecting the temporal correlation of the GPS surface displacement records with the water storage. This quantity, h 1 + h 2 , is referred to as the WSI, and, we assume that it is proportional to the water content per unit area of the catchment.
To adapt the model to the present study areas, the model parameters must be calibrated. This calibration is achieved by minimizing the mean square error between the model discharge (the observable total side exit discharge) and the recorded volume flow rate. Mathematically, it reads where T 0 , T f are the initial and terminal time of the inspected period and p is the parameters to be calibrated. Symbol q obs is the monitored volume flow rate that flows out of the catchment. In our study areas, q obs was the volume flow rate measured at A-Chy-Ba Bridge for the TienChih catchment and at Hsin-Wu-Lu for the SiangYang catchment. Precaution must be taken because even though the present tank model is simplified from a standard three-tank model, there are still ten parameters to be calibrated. These parameters, represented by p in the above expression, include three discharge coefficients (a 1 , a 2 , a 3 ), two infiltration coefficients (b 1 , b 2 ), three controlling threshold depths (H 1 , H 2 , H 3 ), and two initial water storage levels (S 1 , S 2 ). The extensive number of parameters makes the calibration a ten-dimensional optimization process and can be computationally expensive. To improve the efficiency of calibration, a global optimization method called the shuffled complex evolution method (SCE-UA) is used to search the best set of parameters, Chen et al. (2013). The SCE-UA method was developed in the 1990s, Duan et al. (1994), and has been widely applied to many hydrological-related studies, for example, characterization of catchment properties by Nasonova et al. (2009), Funke (1999, and Rustanto et al. (2017); algorithmic and mathematical aspects of the method and similar aspects by Kuczera (1997), Singh and Bárdossy (2015) and Seong FIGURE 6 | Conceptual sketch of tank model. The model has two tanks with an identical area A. Rainfall intensity (RI), given by R(t), is the water input source (h 1 , h 2 ) are the tank water levels, (q 1 , q 2 , q 3 ) are discharge rates, and (q 12 , q b ) are infiltration rates. Parameters (a 1 , a 2 , a 3 ) are discharge coefficients and (b 1 , b 2 ) are infiltration coefficients. The initial values of (h 1 , h 2 ) are denoted as (S 1 , S 2 ) (not shown in the figure). The total water depth, h 1 + h 2 is referred to as the water storage index (WSI) in the present study. H (·) is the Heaviside function to ensure uni-directional water discharge and infiltration. In applications, A is the area of the catchment.
Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 743669 et al. (2015); and a comprehensive review of its adaptations by Rahnamay Naeini et al. (2019). The search for the global optimum is based on the concept of competitive evolution of clusters of parameter sets, and is intrinsically a combination of probabilistic and deterministic approaches. In the present calibration scheme, we use a simple least square error estimator for a better resolution of the peak discharge rates between the model and measurement. In the following sections, the calibration and temporal correlation analysis of the TienChih landslide will first be described in detail, and then, the correlation analysis of the SiangYang landslide will be presented.
With the simple phenomenological quantitative tank model, the following two sections show that promising results can be obtained. However, we need to point out a few important precautions that should be taken when utilizing hydrological models. First, the tank model is suitable in principle for small-to medium-sized catchments. According to Okada et al. (2001), the maximum catchment size is conservatively suggested to be 10 km 2 , although occasional applications to large catchments have been seen, e.g., ∼ 600 km 2 in Chen et al. (2005). In this regard, the application of the model to the two areas is a somewhat exaggerated overextension. Under these circumstances, the spatial transports of surface runoff and underground water are omitted. Exaggerated applications simply arise because of the lack of field data. In the mountain ranges of the island, only a few rain gauges and flow metre are installed to monitor the catchments. Therefore, aggressive tradeoffs have been made among the model accuracy, complexity and acquisition of phenomenological information.
Second, we tested the parameter calibration using standard three-tank models but the improvements were negligibly marginal. A heteroscedastic maximum likelihood estimator [HMLE; Sorooshian (1981)], has also been tested with the calibration procedure but this scheme usually led to underestimation of the peak discharge rates. Last, WSI, as the sum of water storage in the tanks, is expected to correlate with the underground water level. Unfortunately, there is no underground water information available for the two presented sites. It remains to be answered in future studies whether underground water information can be adopted as an additional constraint for parameter calibration and what efficacy can be obtained by doing so.

ANALYSIS FOR THE TIENCHIH LANDSLIDE
For the TienChih landslide, rainfall records are taken from the Nan-Tien-Chih rainfall station, and the flow rate of the Laonung River is measured at the A-Chyi-Ba Bridge station. The locations of the two hydrological stations can be found in Figure 2 (caption) and 1. The catchment area in the calibration process is 411 km 2 .
A typical result of the optimized tank model is shown in Figure 7. The result is for the year 2009. The discharge contribution of each tank in the model to the river flow rate is plotted in Figure 8. The discharge from the top tank is responsible for the peak river flow rate immediately after each rainfall event while the discharge from the bottom tank contributes to the lagged recessing river discharge. This  To streamline the analysis, we do not select any particular time series data for any individual landslide event. That is, the same procedure is applied annually to 2004, 2005, 2006 and 2009 (provided that 2009 has only 9 months of data). The optimized parameter sets for the years in which the TienChih GPS-recorded landslide motion is tabulated in Table 1. Because of the nonstationary regional climate system and inhomogeneous natural environment, there are variations among the four sets of model parameters. The unavoidable scarcity of data sources (one rainfall gauge and one flow station for each catchment) renders that model calibration suffers larger parameter variations among different sampling periods. Nevertheless, the results of a year-on-year comparison of the calibrated parameters are described as follows: 1) The discharge thresholds and initial water storage are on the orders of a unit magnitude, i.e., between 0 and a few tens, and 2) H 2 is close to zero. Figure 9 shows a simple sensitivity test for the model using the 2009 parameters for the 2006 data. The simple sensitivity test was performed to check if the characteristics of the river flow hydrograph remain. The hydrograph characteristics are the essential information for this model to work in the limit of present data scarcity. The sensitivity results show that the timing of peak river discharges is influenced very little, but the values of peak flow rates are underestimated in the case of the 2009 parameters. By inspecting the model sensitivity with other cases, we find that the variations in the peak flow rates are limited within an order of magnitude, or more specifically, within a maximum relative error of 300%. Therefore, caution must be used in regard to the parameter effects if the model parameters are not well calibrated. To illustrate the model sensibility further, we will reveal the result of the tank model with the 2009 mismatched parameter WSI-2009 in the following discussion.
Next, we inspect the relationship between the WSI and GPSrecorded surface displacement. For the four distinct landslide events, the related time series data, RI, WSI, and displacement, are sketched in Figure 10. Notably, the WSI in 2004 had a much longer water recession period which means a much lower discharge rate than the others. This difference occurs because the second tank has very low discharging rates (a 3 and b 2 ) after calibration for this year. Without other available auxiliary measurement data to verify this WSI estimation, we need to be cautious when interpreting its implications or adopting the model towards regional mitigation systems. For simplicity, we neglect the effects of the water recess periods but address the rising WSI on the landslide motion. A straightforward visual inspection suggests that the triggering threshold of the landslide motion is approximately WSI 350 ∼ 400 mm, as indicated by the light-blue shaded marks in Figure 10. The main criterion is to fit the WSI and all of the four aforementioned landslides: the higher bound of the threshold range is about the lowest possible value that the landslide occurs; and the lower bound is about the highest possible value that the landslide does not occur. The WSI-2009 curves are also superposed with their corresponding events for sensitivity comparison, and the same threshold is applicable. A similar RI threshold can be obtained, and the threshold is approximately 40 mm/h. However, based on the RI, the 6/19/2005 event is not detected, as shown in Figure 10B,  Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 743669 9 because precipitation duration is not properly weighted in this approach.
We have to keep in mind that the WSI threshold should be treated as a suggested level and is in general site dependent. It is because the determination of the threshold is by a simple and subjective comparison between the historical landslide events and WSI hydrograph, on the assumption that the WSI in slopes is directly associated with the variations of the material strengths without time lags. The threshold may be different if the model detail or monitoring network is changed. The accuracy of the threshold is also hard to assess in the present conditions because the number of the verified landslides does not reach a level of statistical significance. In addition, we also neglected the possibility that the threshold may change after each landslide because its change may involve subtle changes in geological conditions and they are not quantifiable in the current framework. Nevertheless, we illustrated the possibility to construct a WSI threshold for a landslide site from the aspects of presented model and data.
Most of the mainstream mitigation warning systems are based on regional statistical regression analyses using rainfall data, for example the peak RI, and duration. In this type of model, the infiltration of precipitation is usually neglected. With the present tank model, infiltration can be taken into consideration. Therefore for comparison, we evaluate the predictions between the RI and WSI approaches. We use the 2009 event for this assessment because it has a continuous displacement record to cover the entire motion. We assume that the rate of displacement is proportional to the water storage, i.e., WSI in the slope because water directly reduces the strengths of the materials; thus, it is expected that the displacement rate (DR) is synchronized with the WSI. To investigate this assumption, we apply a correlation Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 743669 function to calculate the phase difference, or time-lag, between the WSI and DR. The correlation function is defined as follows: as corr( DR, WSI ) E( DR(t) · WSI(t + τ) ), which is the expectancy value, E (·), of the product of the DR and a time-τ-shifted WSI. For comparison, the correlation functions between the DR, first tank water depth h 1 and RI are also calculated by replacing WSI in the above expression with R(t) and h 1 (t). In the calculation, we use 1-month time series data starting from 7/28/2009 and the close-up time series, including RI, displacement, DR, h 1 and WSI are shown in Figure 11. The three correlation functions are plotted in Figure 12. In the plot, the maximum values of the correlation functions indicate the similarities in shapes and timescales of the correlated signals and their corresponding time shifts imply time lags between the correlated signals. From the time shift, we find that the peak of the WSI correlation function occurs at approximately 13 h while those of h 1 and RI correlation functions are at − 22 and −34 h. The negative time shifts indicate that h i and RI events temporally occur prior to the DR. The alternation of the time shifts, from RI, h 1 to WSI, reveals the role of water infiltration to the landslide. The timing sequentially approaches the landslide displacement from the precursor rainfall, water accumulation in the shallow surface layer, and to the deep layer. This sequence confirms not only the assumption that the peak DR occurs closer to the predicted peak of WSI than h 1 and RI but also the phenomenological fact of the tank model that the deep-seated landslide involves water in the deep bed rock layer. Last, because of the water infiltration mechanism, the WSI has a much longer water recession period than rainfall, and subsequently, this recession period results in a lowering and broadening peak of the correlation function.

THE SIANGYANG LANDSLIDE CASE
The SiangYang landslide was triggered by the 2009 Morakot typhoon. By using the same technique described in Section 3, the tank model parameters for 2009 are computed and listed in Table 2. The calibrated runoff hydrograph compared with the monitoring data of the Hsin-Wu-Lu River flow station is shown in Figure 13. To investigate the occurrence sequence of the landslide more precisely, close-up hydrographs of the RI, displacement, DR, h 1 and WSI are plotted in Figure 14.
Comparing the calibrated parameters to those of the 2009 TienChih catchment, the discharge coefficients, a 1,3 , are roughly an order of magnitude larger (10×), indicating that the region has a higher water mobility, i.e., lower water retention (holding) capacity. Consequently, the discharge recession period is short and the hydrological characteristics of the catchment are dominated by those of the first tank. Verification is shown by the fact that the WSI is almost identical to h 1 , as shown in the bottom panel of Figure 14. The same WSI threshold in the TienChih landslide (350 ∼400 mm) is also found applicable to the SiangYang landslide. However, we cannot yet conclude if the threshold fits tightly to the landslide initiation under the present scope of data. The constraints of the threshold discussed in Section 4 still hold to the SiangYang site. Under this circumstance, we should treat the threshold conservatively as a necessary but not a sufficient condition for triggering the SiangYang landslide motion. Nevertheless, the coincidence of the threshold in both the SiangYang and TienChih landslide sites may imply that it may be insensitive to a wide range of geological conditions. This implication is important for  Similar to the TienChih case, the time series of the GPS displacement-related data are also missing for approximately one and a half months after the Morakot typhoon, as shown in the middle panel of Figure 14. To proceed with the calculation of the temporal correlation functions, we assume that the slope starts to move immediately after the time when the GPS data are missing, i.e., on September 8, 2009, and the motion lasts for approximately 5 days by referring to the TienChih landslide, cf. Figure 11. Over the 5-day duration, the area catches most of the precipitation from the Morakot typhoon. Under this assumption, we use two constant position values and a smooth sinusoidal function to synthesize the landslide displacement, which is a quarter wave length sinusoidal function with both ends joining the constants on its both sides, as shown by the blue line in the middle panel of Figure 14. The constant position values are the average position of the GPS station before and after the landslide. Based on the synthesized displacement, the DR is calculated, as shown by the red line in the same figure panel. Because the displacement of the SiangYang landslide is only approximately one-tenth of that in the TienChih case, the GPS displacement data have a low signal-to-noise ratio. To avoid employing further sophisticated signal processing techniques, we use the synthesized displacement for the temporal correlation analysis. Similar to Section 4, 1-month time series data starting from 7/28/2009 are used in the calculation.
The correlation functions for the SiangYang landslide are calculated and presented in Figure 15. Our attention is drawn to the time shifts between the hydrological variables and DR, which are −4.7, −6.3 and −19 h for the WSI, h 1 and RI, respectively. All negative time shift values indicate that the    Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 743669 hydrological hydrographs are ahead of the DR. The alternation of the time shifts follows the same order we observed in the previous section, from RI, h 1 to WSI, revealing the local water cycle movement. Similar to the TienChih landslide, the DR is closest to the WSI hydrograph, implying that the landslide involves water in the deep layer. Because the hydrological characteristics of the site are primarily dependent on the first tank, the WSI and h 1 correlation functions are also nearly identical and their time difference is only approximately 1.6 h. Finally for this case, some artefacts have been introduced into the correlation functions because of the use of the synthesized displacement. These effects include the noise-free well-defined shape of the function and the nearly perfect correlation ( ≈ 0.9) at the function peak.

DISCUSSION
Our attention was drawn to the temporal correlation between the creeping landslide motion and water content in slopes. The TienChih and SiangYang landslides were chosen because they are rainfall-triggered, deep-seated creeping landslides, and their locations are in the upstream-most catchment areas. Because of the accessibility of the sites and the associated mountain range, field measurements for calibrating sophisticated field models are scarce. As a trade-off, we adopted a phenomenological quantitative hydrological model, the simplified tank model, to model the local water cycles. The water storage in the solid slope was assumed to be proportional to the sum of the water depths of the two seriallyconnected tanks, h 1 + h 2 , which is referred to as the WSI. In the classical literature of the model, the WSI is called the soil water index, or SWI. We deliberately avoided using the term because it is intuitively associated with shallow soil-based debris landslides. By virtue of the phenomenological tank model, no restrictions are imposed on the material types, detailed distribution configuration of the materials, or physical dimension of the slopes. These physical and geological effects are indirectly incorporated into the parameters of the model through the parameter calibration processes.
Because of the active tectonic motion Lee and Wang (1985) and Lai et al. (1989), rich networks of fracture and cleavage structures are developed in the mountain bedrock. Borehole investigations and rock core samples in the neibouring area with similar geological conditions reveal that fractured geostructures persist at least to one hundred metre from the surface Lin and Lin (1999); Chinese American Tech. Corp. et al. (1992). Therefore, it is assumed that water can infiltrate deeply into bedrock through these open disjointed networks during heavy rainfall periods. As a result, water accumulation in the slope lifts the groundwater table, leading to an increase in the pore pressure, a reduction in the effective stresses and consequently triggering landslide motion. Based on this hypothesis, the occurrence of landslide motion is expected to correlate closely with the WSI and has been verified with the two presented landslides.
In the tank model, the arrangement of the side water exits is conceptually associated with the discharges from the surface runoff, surface and deep bedrock (or subsoil) layers according to Sugawara (1995). In this arrangement, H 1 must be higher than H 2 , as shown by the calibrated parameters of the two presented catchments, Tables 1, 2. Having confirmed the configuration of the discharge exits of the first tank, it is noticed that the discharge coefficient of the top-most exit, a 1 , is smaller overall than a 2 . This result implies that the discharge of the first tank is primarily controlled by the surface layer, i.e., the discharge through the lower a 2 exit.
Comparing the 2009 calibration results, we find that the two sets of parameters lead to one distinct difference in the local water cycles. The TienChih catchment has a smaller b 2 , such that it keeps water above the second tank for a longer period and causes a longer recession period of the flow rate of the Laonung River. In contrast, the SiangYang catchment has a larger infiltration capability such that water can quickly infiltrate through the second tank into the deeper bedrock layer. Despite the aforementioned distinct hydrological characteristics, the WSI threshold of 350 ∼400 mm is applicable for both landslides.

CONCLUSION
In this paper, the rainfall-triggered TienChih and SiangYang deep-seated creeping landslides were compared with the regional water cycles. The landslide displacement motion was recorded by continuous GPS stations, and the regional water cycles included hydrological variables, such as precipitation, water storage in the catchment and runoff flow. The aim of this study is to illustrate that the infiltration of surface water and water storage in the catchment were associated with landslide motion. The water cycles were modelled using the phenomenological tank model, which was composed of two serially-connected tanks. The model contained a total of ten parameters to be calibrated against the rainfall and river runoff flows. The water content in the catchments, i.e., on the solid slopes, was quantified by using the water depth of the first tank (h 1 ) and the total water depth, (WSI, h 1 + h 2 ).
The association between landslide motion and water cycles was investigated by their correlations. The correlation functions, and thus the time lags, between the DR and the hydrological variables (RI, h 1 and WSI) were calculated. Comparing the time lags, we found that the RI lead ahead, followed by h 1 , and finally, the WSI correlated the closest in time to landslide motion. This sequence clearly indicated the infiltration processes where the water was transported from the surface, through the shallow soil layer, and into the deep bedrock layer. When a certain depth was reached, the water weakened the slope material and triggered landslide motion. Therefore, we can expect that the WSI and landslide motion correlates the best in time but not their appearing order because the detail underground water flows and material weakening mechanism were not modelled. On the direct relationship between the WSI and landslide motion, we suggested a WSI threshold for the investigated deep-seated landslides that was between 350 and 400 mm. The same WSI threshold is applicable to both landslides lead to an important Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 743669 implication that it may be insensitive to a wide range of geological conditions. This proposition needs to be carefully examined with more landslide cases in future. To summarize, the hydrological technique presented in the paper provides the interrelation among the hydrological variables, and in turn, the WSI hydrograph and its threshold can act as valuable accessorial quantities that are beneficial to accurately estimating deepseated landslide initiation in future hazard mitigation warning system.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
C-YK, R-FC and Y-JH contributed to conception and design of the study. C-YK, S-EL and Y-HC developed and executed the hydrological model. S-PL, R-FC and Y-JH provided GPS and hydrological data. C-WL, K-JC and R-YW provided geological investigation details. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This work is supported in parts by the Bureau of Soil Water Conservation, Taiwan, under grants SWCB-105-128 and SWCB-110-028.