Investigating Sediment Dynamics in a Landslide-Dominated Catchment by Modeling Landslide Area and Fluvial Sediment Export

Abstract

Few models are capable of simultaneously simulating the sequences of landslide occurrence and sediment export. Quantification of the temporarily stored sediment within the watershed plays a key role to link hillslope landslides with fluvial sediment export. In this study, two coupled models were proposed to simulate time-series total landslide area and the subsequent sediment export on a daily basis with only the inputs of rainfall and runoff. The landslide model considers per-existing and models new landslide, and the sediment transport model incorporates a sediment storage variable. The landslide and sediment transport model were well evaluated with Nash-Sutcliffe efficiency (EC) of 0.89 and logarithmic Nash-Sutcliffe efficiency (EC_log) of 0.90, respectively, in the Tsengwen Reservoir watershed in southern Taiwan by using long-term observed data (2005–2015). It is found that reactivated landslides were up to 72% of the pre-landslide area, which contributed sediment comparable to the new landslide. Besides, the landslide model indicates that pre-landslide area controls the total landslide area but when rainfall is large it takes control in turn. With the simulation of sediment storage, the sediment transport model can well simulate the sediment export after the catastrophic event (typhoon Morakot in 2009). During the post-Morakot period, small rainfall and runoff can lead to high sediment export owing to the storage of Morakot-triggered landslide. This model will be a useful tool to diagnose the sediment dynamics in the watershed.

1. Introduction

Landslides and the consequent fluvial sediment export play an important role in landscape evolution [1,2], physical and chemical weathering processes [3], and costal ecosystem [4,5]. Recent studies also highlighted that the land-to-ocean sediment export conveying terrestrial carbon stock to the ocean [6,7] is tightly related to the fate of organic carbon [8,9], which is the key to understand global carbon cycle [10,11]. Erosion rate, which is commonly estimated by either landslide area or sediment export, is usually taken as a surrogate linking to the above-mentioned phenomena. However, there is an implicit linkage between landslide (sediment supply) and sediment export. Sediment export depends on the interaction between sediment supply and stream power or the magnitude of rainstorms [12,13]. Sediment supply is mainly dominated by slope failure during extreme events [14] and stream runoff determines the amount of sediment that can be flushed off. Oceania, featured by high-standing mountainous watersheds, has the highest erosion rate and is a hotspot of global land-to-ocean sediment export [15,16]. Taiwanese rivers, have long been considered analogous to the rivers of Oceania [17,18], are the best sites to resolve the relations between landslides and sediment export.

The chain reactions, from rainfall to sediment export, are such complex processes that few models could have simulated. Some studies have demonstrated that downstream sediment export does not necessarily synchronize with the concurrent upstream landslide but depends on the interaction of residual sediment and stream runoff. Sediment export should be more directly controlled by the accumulation and depletion of the stored sediment in the watershed. Fuller et al. [19] introduced a stochastic model to simulate the temporal variation of suspended sediment by using synthetic rainfall to determine river runoff, landslide, and sediment export. However, no observed data are used to verify the model, which makes the application of the model skeptical. Furthermore, estimating sediment export should not solely consider the concurrent contribution of landslides but also the lingering effects of previous landslides. Milliman et al. [20] demonstrated the effect of Tsaoling landslide in the central of Taiwan, in which they indicated that the residual sediment did influence the subsequent sediment export. Van Sickle and Beschta [21] have proposed a sediment storage variable in the channel system to augment the sediment concentration that is estimated by the conventional rating curve in the form of power function. They concluded that considering the sediment storage could improve predicting the sediment concentration during storm events. These findings suggested that successful simulation of sediment export should consider not only sediment supply from landslides but also the transient sediment storage.

In this study, a conceptual model was proposed to simulate both the total landslide area and the sediment export. Firstly, the landslide model, which applied daily rainfall amount to model new landslides and the reoccurrence of landslide scars, estimated the total landslide area. Secondly, the time-series sediment storage was obtained by an empirical landslide volume-area formula proposed by Hovius et al. [22] and Malamud et al. [23]. Thirdly, the sediment export was calculated using a modified sediment rating curve in which the sediment storage was incorporated [21,24]. The purpose of this study is to propose a model that is capable of estimating time-series total landslide area and sediment export by solely inputting daily rainfall data and river runoff.

2. Materials and Methods

2.1. Model Development

In this study, a time-series mass balance model is proposed to simultaneously consider rainfall, rainfall-induced landslide, and subsequent sediment export. The model calculates the sediment input from landslide and the output through river runoff in the channel system (Figure 1). First of all, rainfall-induced landslide area is simulated by inputting rainfall data. Landslide material (i.e., available sediment or sediment storage in this study) estimated from the landslide area is then stored in a sediment storage variable, which represents a temporary stock for the generated sediment in the catchment, and the storage is modified as landslide supplies and river delivers. Next in order, sediment export is simulated using a varied rating curve in which the coefficients are adjusted according to the current sediment storage.

2.1.1. Landslide Model

Total landslide in this study is divided into two categories, i.e., new landslide and reactivated landslide, analogous to Lin et al. [25]. At each snapshot of landslide map at time i (usually after typhoon events in this study), the new landslide area is defined as the product of the total landslide area and the new generation ratio (Equation (1)). Additionally, the reactivated ratio is defined as the ratio of reactivated landslide area at time i to the old landslide area derived at time i − 1 (Equation (2)). New landslide and reactivated landslide could be delineated by overlaying the landslide map taken at time i on the one taken at time i − 1. The overlapped landslide area is defined as the reactivated one; on the other hand, the dissimilarity at time i is outlined and defined as new landslide. The new generation ratio and reactivated ratio of landslide was defined as in Equations (1) and (2). Therefore, the total landslide area at time i can be calculated by Equation (3), where A_i, A_i,re, A_i,new stand for the total, reactivated, new landslide area (km²) at time i, respectively. a_LA and A_i−1 are reactivated ratio of landslide and total landslide area at time i − 1, respectively.

$$\mathrm{New}\mathrm{generation}\mathrm{ratio}=\frac{\mathrm{New}\mathrm{landslide}\mathrm{area}\mathrm{after}\mathrm{the}\mathrm{event}}{\mathrm{Total}\mathrm{landslide}\mathrm{area}\mathrm{after}\mathrm{event}}\times 100\%=\frac{A_{i,new}}{A_i}\times 100\%$$

(1)

$$\mathrm{Reactivated}\mathrm{ratio}=\frac{\mathrm{Reactivated}\mathrm{landslide}\mathrm{area}\mathrm{after}\mathrm{the}\mathrm{event}}{\mathrm{Total}\mathrm{landslide}\mathrm{area}\mathrm{before}\mathrm{the}\mathrm{event}}\times 100\%=\frac{A_{i,re}}{A_{i-1}}\times 100\%.$$

(2)

$$A_i=A_{i,re}+A_{i,new}=a_{LA}{A}_{i-1}+A_{i,new}$$

(3)

In Equation (3), new landslide area is the function of rainfall [26]. According to Chen et al. [26], maximum 24-h rainfall is the most representative rainfall characteristic to determine landslide erosion in the study watershed. Unfortunately, due to our rainfall data limitation, maximum daily rainfall during the event was used instead. The relation of new landslide area and rainfall is shown as Equation (4), where $R_i$ is the maximum daily rainfall and $R^{*}$ (=382.5 (mm)) is the long-term average of annual maximum daily rainfall at Alishan (X46753) and Tsaoling (C1M40) gauge. According to Soil and Water Conservation Bureau, Taiwan, rainfall triggers landslides in mountain areas at a minimum threshold amount of 200 (mm). We hereby define new landslide can only be generated when R_i/R* > 0.52 (i.e., 200 (mm)/382.5 (mm)). b_LA and c_LA are model parameters that need calibration. Combining Equations (3) and (4), the total landslide area at time i can be rewritten as Equation (5).

$$A_{i,new}=f\left(\mathrm{Rainfall}\right)=b_{LA}{\left(\frac{R_i}{R*}\right)}^{c_{LA}}$$

(4)

$$A_i=a_{LA}{A}_{i-1}+b_{LA}{\left(\frac{R_i}{R*}\right)}^{c_{LA}}$$

(5)

The total landslide area is then converted to the amount of the sediment storage that could be transported to downstream by river runoff. However, the conversion ratios of landslide area to sediment storage should be different for new landslide and reactivated landslide. For new landslide, Malamud et al. [23] allows the calculation of the average volume of one landslide (in the curly brackets of Equation (6)) and we multiply by the number of landslides, N_L, to get the total volume of the triggered landslide population, V (km³). In Equation (6), ε and γ are intercept and scaling exponent of landslide V-A relation (see more details in [23]). According to our observed landslide data, fairly good relation ($N_{L,i}=186.4ln\left(A_{i,new}\right)+350.6$, R² = 0.78) exists between number (N_L) and area of new landslide. In Equation (6), the values of inverse-gamma parameters ρ (–) and a (km⁻³) are 0.87 and 1.3 × 10⁻³ derived from observed landslide data, and ε and γ are 0.0082 and 1.268 adjusted from the results of Chen et al. [26].

$${\mathrm{V}}_{i,new}=\left\{\frac{\epsilon a^{\gamma }}{\Gamma \left(\rho \right)}\left[\frac{1}{\gamma -\rho }{\left(\frac{N_{L,i}}{\rho \Gamma \left(\rho \right)}\right)}^{\frac{\gamma -\rho }{\rho }}-10.69\right]\right\}\times N_{L,i}$$

(6)

For reactivated landslide, the landslide V-A relation is not used for volume estimation but the thickness of residual sediment. Reactivated landslide volume, $V_{i,re}$ (km³), is presented as Equation (7), where H is the average landslide depth, i.e., 1.7 (m) from field survey [27], and C (<1) is a residual ratio adjusting the reactivated landslide volume that needs calibration.

$$V_{i,re}=a_{LA}{A}_{i-1}\times H\times C/1000$$

(7)

2.1.2. Sediment Transport Model

Fluvial sediment export, which is controlled by sediment supply (i.e., sediment storage) and carrying capacity of the stream (i.e., river runoff), is usually calculated as a function of runoff (Q) in the form of a power function, i.e., aQ^b, without explicitly considering the sediment supply. According to Van Sickle and Beschta [21], the conventional sediment rating curve is adjusted with considering sediment supply, i.e., S_i, as shown in Equation (8), where $Q_{s,i}$ (ton day⁻¹) represents sediment export at time i in the condition of the corresponding river runoff $Q_i$ (m³ s⁻¹) and sediment storage $S_i$ (mm). $S_{mean}$ is the long-term average erosion depth, 4.78 (mm), derived from reservoir siltation record in the downstream of the study watershed [26]. $a_{SC}$ and $c_{SC}$ are analogous to the intercept and slope parameter, respectively, in the conventional rating curve, and need calibration in this study. $b_{SC}$ is a model parameter regulating the supply rate of the sediment storage and is derived by calibration as well. It should be noted that rating curve method is applied to estimate suspended load. Moreover, previous studies show that the fraction of suspended load in Tsengwen Reservoir watershed is close to 1 [28]. Therefore, the estimation of bedload export is neglected.

$$Q_{s,i}=a_{SC}{Q}_i^{c_{sc}+1}{\left(\frac{S_i}{S_{mean}}\right)}^{b_{SC}}$$

(8)

2.1.3. Sediment Storage

The amount of sediment storage is modulated by the input from landslide at time i (S_in,i) and the output through river at time i − 1 (S_out,i−1), which can be expressed by Equation (9).

$$S_i=S_{i-1}-S_{out,i-1}+S_{in,i}$$

(9)

S_in is derived from landslide volume estimation (Equations (6) and (7) in the unit of km³) and S_out is derived from fluvial sediment export calculation (Equation (8) in the unit of ton day⁻¹). S_in and S_out should have identical unit before implementing Equation (9). Therefore, S_in is calculated by dividing the summation of Equations (6) and (7) by the watershed area (Area (km²) in Equation (10)) and the unit is converted to mm. S_in,i is the sediment contributing to sediment storage from landslides at time i in the condition of rainfall larger than 200 (mm), as shown in Equation (10). In Equation (9), the initial value of sediment storage, S₀, is determined as the landslide volume after typhoon Mindulle in 2004, which is the first influential typhoon after 2001.

$$S_{in,i}=\frac{{\mathrm{V}}_{i,new}+{\mathrm{V}}_{i,re}}{Area}\times {10}^6$$

(10)

On the other hand, the amount of fluvial sediment export at time i, S_out,i (mm), can be calculated by cumulative volume of exported sediment during the time interval [0,T] (days) (Equation (11)). In Equation (11), D is sediment density ranging between 1.6–2.5 (ton m⁻³) (unpublished data).

$$S_{out,i}=\frac{{\int }_{t=0}^T{Q}_{s,i,t}dt}{D\times Area\times {10}^3}$$

(11)

2.1.4. Calibration of Landslide and Sediment Transport Model

The integrated model proposed in this study is designed to simulate time-series of total landslide area (and thus volume) and sediment export on a daily basis. The observed data of total landslide area and event-based sediment export in 2005–2014 were used to estimate the parameters by simultaneously calibrating the two models. Over 100,000 parameter sets were generated by the uniform distribution in fixed ranges for each parameter (a_LA: 0.6–0.8; b_LA: 0.3–0.5; c_LA: 3–5; a_SC: 0.1–0.3; b_SC: 0.7–0.9; c_SC: 0.3–0.5; C: 0.8–1). Nash-Sutcliffe efficiency coefficient (EC) was used as the performance measure for landslide model. The value of EC ranges from negative infinity to 1.0 which represents perfect simulation. Nash-Sutcliffe efficiency coefficient with logarithmic value (EC_log) was used for sediment transport model to assure the simulation for the low values which is overlooked by EC. The non-inferior solutions, based on EC and EC_log for the two models, were selected by Pareto front [29]. As a result, with the help of this integrated model, one can obtain time-series simulated landslide volume and fluvial sediment export by only inputting the daily rainfall and river runoff.

2.2. Study Area

The study watershed is in the upstream of the Tsengwen Reservoir in southwestern Taiwan, with the Dapu gauging station as the watershed outlet (Figure 2). Tsengwen Creek drains a watershed of 314 (km²), which originates from the Ali Mountain. The watershed elevation ranges from 126 to 2609 m a.s.l., and the average slope is 26.1°. The land use type of Tsengwen Reservoir watershed is mostly composed of forest and agriculture, which accounts for 86% and 8%, respectively. The average annual rainfall in the study period (2005–2015) was around 3516 (mm) and 78% of the annual rainfall occur in wet season (June to October). Geological setting mainly comprises the alternation of sandstone and shale, making the Tsengwen Reservoir watershed highly erodible. The average landslide ratio in the study period is ~1.1%, ranging from 0.4–3.3%.

Given the conditions of frequent typhoons, highly erodible geology, high relief and steep gradient, Taiwan is featured by extremely high sediment production rates [30,31] and one of the most applicable places to landslide and sediment issues. Of which, reservoir siltation is one of the most serious problems to the reservoir management, particularly after Typhoon Morakot bringing a record-breaking rainfall in 2009. More than 3500 (mm) rainfall poured down in 5 days, resulting in more than 700 landslides and a 14% reduction in the reservoir capacity. The amount of siltation is equal to the total amount accumulated over the past 35 years.

2.3. Landslide Data

The yearly landslide maps taken by Formosa-2 (FS-2) are available on the website (http://246gis.swcb.gov.tw) provided by Soil and Water Conservation Bureau, Taiwan. Those maps reveal the landslide inventory for every year, which basically represent the collective impact from typhoons in a year. The FS-2 was launched and placed in a daily revisiting orbit in 2004. It was equipped with a high-spatial-resolution sensor providing panchromatic images at the resolution of 2 m and multispectral band images at 8 m. Therefore, FS-2 potentially can provide daily time-series satellite images that are capable of showing the instantaneous changes of land surfaces, e.g., the evolution of landslide after typhoon events. The image processing procedures and quality control regarding landslide delineation can be referred to Liu [32] for details. The following is a brief introduction to it. First, the raw FS-2 images are preprocessed automatically by the “Formosat-2 image processing system” (F2 AIPS), including band-to-band coregistration, orthorectification and geometrical registration of multi-temporal imagery. Second, each landslide patch was identified by integrating the Getis statistic, spectral index, and unsupervised K-means classification analysis [33]. The overall accuracy of automatically delineated landslide area can achieve 99.04% compared to those manually processed ones. It should be noted that the yearly landslide map of the study watershed was a compilation of several satellite images taken on adjacent dates to minimize cloud cover effects. After examining the date of a used image on a map, it is assumed that the landslide inventory in that year (e.g., at time i) is attributed to the preceding maximal rainfall event occurring after the date of image in the previous year (i.e., at time i − 1). There are 8 typhoon events and 2 rainstorm events shown in

Table 1. Basic information for the observed total landslide area and their triggers, i.e., typhoons or rainstorms.

Rainstorm	Date	Max Daily Rainfall (mm)	Max Daily Discharge (m3 s−1)	Date of Landslide Map	Total Landslide Area (km2)
Haitang	20/7/2005	437.4	2475.2	27/6/2006	2.70
Bilis	15/7/2006	427.1	2698.1	3/7/2007	0.94
Krosa	7/10/2007	544.1	2770.1	19/1/2008	2.11
Sinlakut	15/9/2008	554.6	2572.5	10/122008	1.47
Morakot	9/8/2009	887.9	6814.9	1/4/2010	12.53
Rainstorm	19/9/2010	183.0	979.1	16/4/2011	7.65
Rainstorm	31/8/2011	315.5	1175.3	28/6/2012	4.83
Saola	3/8/2012	338.4	1491.1	3/6/2013	3.83
Soulik	14/7/2013	312.4	1615.7	10/4/2014	4.21
Matmo	25/7/2014	272.3	808.3	9/6/2015	2.97

.

2.4. Hydrometric Data and Calculation of Sediment Transport

Long-term sediment concentration and runoff data (1957–2014) were derived from Hydrotech Research Institute, National Taiwan University and Tainan Hydraulics Laboratory, National Cheng Kung University. Of which, event-based sediment samples were taken every 1 to 3 h at Dapu gauging station during 9 rainstorm events (

Table 2. Information for the sampled sediment data from the 9 rainfall events.

ID	Typhoon Name	Sampling Period		Numbers of Sample	Maximum Sampled Sediment Concentration (g L−1)	Maximum Sampled Discharge (m3 s−1)	Rating Curve F(g s−1) = aQb	R2	Bias Correction Factor	Sediment Export (Mt)
		From	To
2005	Haitang	N.A.	N.A.	61	16.32	1208	24.61Q1.89	0.83	0.21	14.97
2007	Sepat	N.A.	N.A.	42	9.01	1307	766.32Q1.27	0.76	0.10	2.21
2008	Sinlakut	14/9/2008 01:00	15/9/2008 18:00	43	14.91	4251	78.19Q1.61	0.93	0.02	2.97
2009	Morakot	N.A.	N.A.	49	68.20	5797	5011.1Q1.28	0.95	0.01	67.57
2012	Rainstorm	11/6/2012 17:00	14/6/2012 14:00	31	38.72	2314	1.92Q2.28	0.91	0.08	9.67
2012	Saola	1/8/2012 14:00	3/8/2012 12:00	19	25.54	2634	142.57Q1.63	0.99	0.11	2.72
2013	Sulik	13/7/2013 06:00	14/7/2013 09:00	15	56.71	2712	520.28Q1.55	0.94	0.03	5.34
2013	Trami	22/8/2013 00:00	23/8/2013 10:00	17	33.60	2303	4.00Q2.13	0.96	0.01	3.33
2013	Kong-rey	29/8/2013 10:00	31/8/2013 12:00	24	18.33	2743	24.91Q1.81	0.94	0.01	2.03

). Depth integrated water samples were obtained using a vertically mounted 1 L bottle attached to a weighted metal frame that was gradually lowered from the bridge. The US Geological Survey DH-48 sampler was not used because of its difficulty in sinking in turbulent flows > 2.5 (m s⁻¹). Total suspended matter (TSM) was determined gravimetrically using pre-combusted 0.7 (μm) pore size GF/F filters [34]. The measured sediment concentration varied from 1331 to 68,197 (mg L⁻¹). The highest measurement occurred when typhoon Morakot invaded. Hourly river runoff was also measured at Dapu station, which is maintained by Water Resources Agency, Ministry of Economic Affairs, Taiwan. The average runoff of the study period was around 48 (m³ s⁻¹), but the daily runoff fluctuated from ~0.1 to 6815 (m³ s⁻¹). Typhoon Morakot led to an observed maximal runoff on 9 August 2009.

Flux estimator is essential in this study to calculate sediment export during a rainstorm event given the conditions of discrete sediment concentration measurement and continuous river runoff data. The rating curve method is one of the most appropriate flux estimation methods and has been widely applied to rivers in Taiwan, particularly suitable for the sediment flux estimation [13,17]. This method presumes that a power function (i.e., F = aQ^b) exists between the observed runoff (Q) and sediment flux (F) that is the product of observed sediment concentration and the instantaneous river runoff. The coefficients of the power function, a and b, can be derived from the observed sediment fluxes and the river runoff by the log-linear least-square method [17]. For the non-typhoon periods, the rating curve, constructed by the long-term sediment concentration and runoff data, was used to calculate daily sediment export by substituting daily runoff into the rating curve (i.e., F(g s⁻¹) = 33.94Q^2.90, R² = 0.84). Moreover, 9 event-based rating curves were established for each individual rainstorm event. Each rating curve was then applied to hourly runoff during the event to calculate hourly sediment export. The summation of the hourly sediment export was the total sediment being flushed off the watershed in the event. Besides, it is well known that the rating curve method inherently underestimates the summation of the sampled sediment export within a concerned period, particularly for the fluxes at the highest runoff. Therefore, a bias-correction factor proposed by Kao et al. [17] was introduced to reduce residuals.

3. Results

3.1. Observed Total Landslide Area and Sediment Export in Tsengwen Reservoir Watershed

Figure 3 illustrates the observations in this study, including daily rainfall, river runoff, sampled sediment data, and yearly total landslide area (named after the 8 typhoon and 2 rainstorm events). Among the events, the observed minimum daily rainfall and runoff occurred during the 2010 rainstorm, which were 183 (mm day⁻¹) and 979 (m³ s⁻¹), respectively, while the highest daily rainfall and runoff reached 878 (mm day⁻¹) and 6815 (m³ s⁻¹) during typhoon Morakot in 2009 (Table 1). As for observed total landslide area, before typhoon Morakot on 9/8/2009, the total landslide area ranged between ~1 to ~3 (km²). The maximum area, 12.53 (km²), was triggered by typhoon Morakot. Afterwards, the landslide area gradually decreased from ~12 (km²) to ~2 (km²). However, the smallest total landslide area was found after typhoon Bilis rather than the smallest rainfall event, i.e., the 2010 rainstorm. Typhoon Bilis had the 5th highest daily rainfall, 427 (mm day⁻¹), but resulted in merely 0.94 (km²) landslide area. Poor rainfall–landslide area relationship (r = 0.38) in Tsengwen Reservoir watershed hampers the application of the conventional regression models [19,26].

Table 2 shows the sediment export during the 9 sampled rainfall events. Maximum sampled sediment concentration and runoff were 68.20 (g L⁻¹) and 5797 (m³ s⁻¹), respectively, during Typhoon Morakot. Event-based rating curved were established for each individual event to estimate sediment export, and the R² ranged from 0.83–0.99. Therefore, sediment export of the 9 sampled rainfall events fell between 2.03–67.57 (Mt).

3.2. Calibration of Model Parameters

The landslide model and the sediment transport model were calibrated simultaneously but were evaluated by different performance measures, i.e., EC and EC_log, respectively. Figure 4 illustrates the Pareto front of EC against EC_log with the 23 red dots representing the non-inferior parameter sets (out of the 100,000 parameter sets). The values of the calibrated parameters are shown in

Table 3. Best-fit parameters for the landslide model and the sediment transport model.

	aLA	bLA	cLA	aSC	bSC	cSC	C	EC	EClog
Range	0.6–0.8	0.3–0.5	3–5	0.1–0.3	0.7–0.9	0.3–0.5	0.8–1	-	-
1	0.70	0.50	3.26	0.20	0.71	0.33	0.88	0.88	0.91
2	0.72	0.39	3.57	0.11	0.72	0.40	0.91	0.89	0.90
3	0.68	0.46	3.46	0.13	0.71	0.38	0.94	0.88	0.91
4*	0.72	0.39	3.62	0.15	0.71	0.36	0.86	0.89	0.90
5	0.73	0.32	3.86	0.19	0.75	0.33	0.96	0.89	0.88
6	0.69	0.48	3.42	0.20	0.70	0.33	0.82	0.88	0.91
7	0.72	0.35	3.73	0.12	0.73	0.37	0.87	0.89	0.89
8	0.72	0.36	3.70	0.10	0.76	0.39	0.93	0.89	0.89
9	0.69	0.42	3.58	0.11	0.72	0.40	0.99	0.89	0.90
10	0.74	0.31	3.89	0.17	0.78	0.33	0.99	0.89	0.87
11	0.73	0.30	3.92	0.16	0.80	0.30	0.98	0.89	0.86
12	0.70	0.43	3.56	0.11	0.73	0.39	1.00	0.89	0.90
13	0.72	0.48	3.19	0.15	0.70	0.37	0.84	0.87	0.91
14	0.68	0.44	3.53	0.13	0.71	0.38	0.81	0.88	0.91
15	0.68	0.50	3.27	0.24	0.71	0.31	0.81	0.88	0.91
16	0.73	0.37	3.65	0.17	0.70	0.35	0.84	0.89	0.89
17	0.72	0.38	3.64	0.14	0.70	0.37	0.89	0.89	0.90
18	0.74	0.35	3.70	0.14	0.70	0.38	0.93	0.89	0.90
19	0.71	0.45	3.44	0.13	0.71	0.38	0.88	0.89	0.90
20	0.70	0.46	3.38	0.18	0.71	0.34	0.91	0.88	0.91
21	0.74	0.30	3.88	0.10	0.70	0.42	1.00	0.89	0.89
22	0.70	0.42	3.55	0.17	0.71	0.35	0.91	0.89	0.90
23	0.72	0.38	3.64	0.14	0.75	0.35	0.87	0.89	0.89
Avg.	0.71	0.40	3.58	0.15	0.72	0.36	0.91	0.88	0.90
Std.	0.02	0.06	0.20	0.03	0.03	0.03	0.06	0.01	0.01

Note: Shaded cells represent the p-value of the parameter is <0.05. *4th parameter set is defined as the optimal parameter set.

with the shaded ones representing its p-value less than 0.05. The EC values of the 23 parameter sets are between 0.87–0.89 for the landslide model, and EC_log are between 0.86–0.91 for the sediment transport model.

For the landslide model, the landslide reactivated ratio, $a_{L\mathrm{A}}$, is between 0.68–0.74 among the 23 parameter sets, which means that more than half of the previous total landslide area would be retriggered. Representing new landslide, the parameter $b_{L\mathrm{A}}$ and $c_{L\mathrm{A}}$ are between 0.30–0.50 and 3.19–3.92, respectively. As for the sediment transport model, the intercept, $a_{SC}$, and slope, $c_{SC}$, in rating curve are 0.10–0.24 and 0.30–0.42. The scaling parameter for sediment storage, $b_{SC}$, is between 0.70–0.80. The contributing ratio, C, is between 0.81–1.00, indicating that the sediment production rate on reactivated landslide is nearly at the same order as the new landslide. Overall, the 4th parameter set is defined as the optimal one because 6 out of 7 parameters are identified as significant and parameter C is not as sensitive as the other parameters (see Section 4.2). The optimal parameters are a_LA = 0.72, b_LA = 0.39, c_LA = 3.62, a_SC = 0.15, b_SC = 0.71, c_SC = 0.36 and C = 0.86, and the EC and EC_log are 0.89 and 0.90, respectively.

3.3. Simulated Time-Series Total Landslide Area

The simulated time-series total landslide area is shown in Figure 5. The gray zone in Figure 5 embraced the simulations from the 23 parameter sets. In the landslide model, landslide was triggered only when daily rainfall exceeded the designate threshold, i.e., 200 (mm day⁻¹). Without considering the recovery of landslide scars, the simulated total landslide area remained a constant before the new landslide generated, resulting in the stepwise curve of simulated landslide area. During the period of 1/1/2005–31/12/2008 (Period 1), the simulated total landslide area was between 1–4 (km²) while the observed total landslide area was between 1–3 (km²). During Period 2 (1/1/2009–31/12/2012) with the invasion of typhoon Morakot, the simulated total landslide area ranged from ~10 to 3 (km²). After typhoon Morakot, the total landslide model well simulated the decreasing trend of total landslide area. However, it slightly overestimated the total landslide area triggered by 9/9/2010 and 31/8/2011 rainstorm events. During the Period 3 (1/1/2013–31/12/2014), the simulated total landslide area gradually decreased from 2 (km²) to 1 (km²), which was close to the level before typhoon Morakot.

3.4. Simulated Sediment Export

Figure 6a illustrates the simulated time-series sediment export (gray line) as well as the observed one (black triangles). The daily sediment export varied at 7-order of magnitude (~10⁻⁴–~10¹ (Mt)), revealing the flushing nature of Taiwan rivers. The low-end values of sediment export shifted from ~10⁻⁵ (Mt) in Period 1 to ~10⁻³ (Mt) in Period 2, and then to ~10⁻⁴ (Mt) in Period 3. The sediment transport model simulated the observed sediment export well with EC_log = 0.90 (Figure 6b). The error bars in Figure 6b represent the 2-fold standard deviation of the simulated ranges among the 23 non-inferior parameter sets. To illustrate the changes of carrying capacity of the river, the simulated daily sediment export was normalized by daily river runoff (i.e., sediment export/runoff, Figure 6c). It is found that the carrying capacity was at the lowest in Period 1 (mean = 0.003 t m⁻³), and it increased to the highest in Period 2 (0.006 t m⁻³), and then decreased to a middle value in Period 3 (0.005 t m⁻³). It indicates that one unit river runoff could deliver more sediment during the post-Morakot period, owing to more sediment storage.

4. Discussion

4.1. Best-Fit Parameter Sets

The 23 non-inferior parameter sets derived from the Pareto front were determined as the best-fit parameter values (Table 3). The reactivated ratio ($a_{LA}$), ranging between 0.68–0.74, is close to the upper range of the observed reactivated ratio (Table 1). Previous studies have addressed the reactivated ratio for watersheds in Taiwan, showing values between 0.40–0.93 in Chenyoulan watershed in central Taiwan [25], and 0.52–0.67 in Sinwulu and Luye watersheds in Eastern Taiwan [35]. It can be concluded that more than half of the total landslide area would be retriggered in Taiwan watersheds where typhoons disturb recurrently every year.

The rating curve intercept parameter (Equation (8)), which is the product of a_sc and ${\left(\frac{S_i}{S_{mean}}\right)}^{b_{sc}}$, ranges from 0.12–4.35 and is compared with the intercept in the conventional rating curve, i.e., a in the equation of aQ^b. The intercept derived in this study is much higher than that of other world rivers, ranging from 0.01–0.04 [36], which indicates that there is considerable sediment being transported during low-to-medium flow period. In other words, in conditions of small rainfall events, the river in this study has sufficient sediment storage to support higher sediment export in comparison with other rivers in the world. Correspondingly, it is essential to consider sediment storage in the model. As for the rating curve slope, parameter c_SC, between 1.3–1.42, agrees with the results of 16 primary rivers in Taiwan (ranging mostly from 1 to 2.5) [17]. The values of slope parameter for rivers worldwide also fall within 1 to 2 [36].

The parameter C (in Equation (7)) ranging from 0.81–1.00 means that the average landslide depth (or erosion rate) for the reactivated landslide is nearly 1.4–1.7 (m). Tsai et al. [37] investigated the erosion rate on 13 reactivated landslides in Shimen Reservoir watershed in north Taiwan and obtained the erosion rate at 0.78–1.06 (m year⁻¹) which was shallower than that of our study watershed, perhaps owing to different geological conditions. The geology of Shimen Reservoir watershed comprises mainly schist and phyllite, while Tsengwen Reservoir watershed consists of more erodible alternations of sandstone and shale (Figure 2).

4.2. Sensitivity Analysis

To understand the model responses to the change of each parameter, the sensitivity analysis is presented here. The sensitivity analysis was implemented via adjusting one of the parameters in the optimal parameter set from −20% to +20% while leaving the rest parameters unchanged, to examine model responses to each respective parameter changes (Figure 7). The model responses are presented by the changes of EC and EC_log.

Figure 7 illustrates the results of sensitivity analysis for the landslide model (Figure 7a–c) and the sediment transport model (Figure 7d–g). The landslide model is most sensitive to the adjustments (+/- 20% of the optimal value) of a_LA and c_LA (black circle in Figure 7a,c), resulting in −55%/−21% and −33%/−100% changes, respectively, on the EC of the landslide model. The model is less sensitive to parameter b_LA, while the EC varies within −25% of the optimal one, i.e., 0.89. The sediment transport model is relatively insensitive to the landslide model parameters. The adjustments would lower the EC_log to 0.86 (i.e., −4.4%) at most.

The results from landslide model are the inputs of the sediment transport model. The four parameters (Figure 7d–g) are irrelevant to the landslide model; hence the black line remains unchanged. The most sensitive parameters to the sediment transport model (red circle) are the slope parameter, c_SC of the rating curve, with the EC_log decreasing to 0.75 (−16%) and to 0.51 (−43%) from 0.90, respectively, in respond to the −20% and 20% adjustments of the optimal value. The parameter bsc that regulates the effect of sediment storage on sediment export is sensitive as well, with the EC_log deceasing to 0.78 (−13%)/0.62 (-31%) for the −20%/+20% adjustments.

4.3. The Controlling Factors of Total Landslide Area

Figure 8 illustrates a contour diagram of the estimated total landslide area as a function of rain factor (the x-axis, i.e., R/R* in Equation (5)) and pre-landslide area (the y-axis, total landslide area in the previous year). The estimated total landslide area was estimated by substituting the optimal parameter sets in the landslide model (Equation (6)). At a given pre-landslide area, higher rainfall would trigger larger total landslide area, particularly when rain factor is larger than 1.5. When a fixed rain factor is given, the total landslide area is positively correlated to the one in the previous year. However, the correlation becomes weaker when rain factor gets larger. Overall, it indicates that pre-landslide area controls the total landslide area when rainfall is small, like the cases in 2010–2013. Rainfall controls the total landslide area when rainfall is large, as seen in the cases in 2005–2008.

Rainfall is long regarded as the key factor to trigger landslide and is used as the only independent variable to estimate landslide area by many previous studies [19,26]. We suggest that in landslide-dominated region considering not only the rainfall but also previous landslide status as controlling factors in simulating total landslide area is necessary.

4.4. The Sediment Storage and Controlling Factors of Sediment Export

The best-fit parameter sets (in Table 3) were used to estimate a range of sediment storage as the gray zone in Figure 9a. It shows that the study period (1/1/2005–31/12/2014) was divided into three periods (i.e., Period 1, Period 2, and Period 3 as mentioned in Section 3.3) to clearly explain the effects of sediment storage on sediment export. The sediment storage ranged from 0–150 (mm) in Period 1 and suddenly increased to as much as 250 (mm) in Period 2, and then gradually decreased to ~100 (mm) in Period 3. The range of estimated sediment storage became large in Period 2 and Period 3, indicating the relatively uncertain estimation of the available sediment after the extreme weather event, i.e., typhoon Morakot. The ranges of the sediment storage also reflect the uncertain estimation of the total landslide area (Figure 5).

Figure 9b illustrates the rating curve of the simulated sediment export against river runoff for the three periods. The exponents of the three rating curves remain at ~1.38–1.39, indicating the similar effect of stream power on sediment export. However, the coefficients of the rating curves vary among the three periods. At a given river runoff, Period 2 would have the highest sediment export, and then Period 3 and Period 1 would have the lowest. The coefficient in the rating curve reflects the sediment supply in the watershed. The coefficients for Period 2 (i.e., 2 × 10⁻⁴) and Period 3 (1 × 10⁻⁴) are one-order of magnitude larger than Period 1 (6 × 10⁻⁵), indicating substantial sediment supply after typhoon Morakot (Figure 6c). The impact of residual sediment was also found in Choshui watershed in central Taiwan after a severe landslide event (i.e., ChiChi earthquake in 1999), which can last for ~6 years [20]. The recession time is comparable to the Morakot case in this study. We believed that the model with the idea of sediment storage can refine the simulation of sediment export in a more validated way.

Sediment transfer, from hillslope to channel and off watershed, is a global issue wherever there is soil erosion, particularly crucial in the landslide-prone region. Sediment deposition will lead to flood plain and channel aggradation and further result in water level rise and sediment-related hazards. The estimation of sediment storage might help to understand the relevant risk. Period 3 in this study is prone to the risk compared to Period 1 owing to more readily-eroded sediment in the hillslope and channel in Period 3. Moreover, landscape evolution can be inferred from watershed-wide erosion rate [1,2]. However, discrepancy exists between the rates inferred from landslide volume (e.g., digital elevation models (DEMs) derived from photogrammetry of aerial photographs) and fluvial sediment export [38]. The sediment storage might help to explain the discrepancy and prolonged/delayed response of sediment transport to landslide [39,40,41]. Zhang et al. [42] demonstrated that intense monsoonal runoff drove sediment transfer and depositional processes until 2 years after the Wenchuan earthquake, implying the applicability of the sediment storage concept (if earthquake-induced landslide model was incorporated). The sediment storage could also help to explain sediment transport is transport-limited or supply-limited [43].

5. Conclusions

Linkage between sediment supply from landslides and export by fluvial system is unclear due to insufficient knowledge about the spatial and temporal pattern of the sediment deposition within the watersheds, especially in landslide-dominated regions, where highly frequent occurrence of landslide affirms the effect of sediment storage not to be ignored. In this study, we established a conceptual model that simulates the chain reaction from rainfall, landslide to sediment export by considering the effects of pre-landslide area and the sediment storage. Rainfall at a higher resolution temporal scale (e.g., hours) could strongly modulate the timing landslide debris reach the rivers from the slopes, affecting both landslide (e.g., depth of the failure surface; velocity; movement typologies; run-out distance) and water erosion (either splash, sheet, rill, or gully erosion) processes. Nevertheless, the sediment storage variable, to some extent, cleverly and implicitly considers the complexity of temporal and spatial distribution of sediment.

In our landslide model, the excellent performance, EC = 0.89, indicates that the incorporation of the reactivated landslide is crucial in the landslide-dominated region like Taiwan who faces 3–5 typhoon invasions every year. The reactivation of previous landslides significantly contributes available sediment as much as the new landslides. However, the reactivation ratio can be improved and should vary among landslide occurrences, seasons, and even among topographic conditions. As for the sediment transport model, it is found that higher sediment storage after substantial landslide input induced by catastrophic typhoon events could elevate sediment export 1-order of magnitude higher than that before the input. Several studies have indicated the effects of temporarily stored sediment on sediment export from the observed C-Q relation [13,19,20]. In this study, we further quantify the amount of sediment storage and its consequent effects. Furthermore, we suggest that the concept of storage could be applied to other stock relevant issues, such as the lingering effects of litter fall on nutrient release, to quantify the effects on biogeochemical cycle [44]. However, bedload transport is one of the limitations in this study, which can be improved by incorporating bedload transportation estimation equations in the model [45,46]. Nevertheless, the validation of bedload transport is always a problem in the field [47].

Global warming is proven in the ascendant and is resulting in more unprecedented weather extremes. The increased rainfall intensity has led to magnified sediment export from Taiwan rivers [13,26]. Our proposed model could benefit the understanding of the sediment dynamics in the watersheds. Pre-landslide area and the current sediment storage are both crucial in controlling sediment budget in the river.