Integrating Landscape Metrics and Hydrologic Modeling to Assess the Impact of Natural Disturbances on Ecohydrological Processes in the Chenyulan Watershed, Taiwan

The Chenyulan watershed, located in the central mountain area of Taiwan, has been suffering from earthquakes, typhoons, and heavy rainfalls in recent decades. These sequential natural disturbances have a cumulative impact on the watershed, leading to more fragile and fragmented land cover and loss of capacity of soil water conservation. In this study, the Soil and Water Assessment Tool (SWAT) and a landscape metrics tool (FRAGSTATS) were used to assess the direct impact (e.g., by annual rainfall) and indirect impact (e.g., by landscape configuration and composition) of natural disturbances on the ecohydrological processes of the Chenyulan watershed. Six SPOT satellite images from 2008 to 2013 were analyzed by using the nearest feature line embedding (NFLE) approach and reclassified into six land cover types: forest, cultivated land, grassland, river, landslide, and built-up. Forest was found to have the largest patch size, indicating that it is more resilient to disturbances, while agricultural land tended to expand from the river side toward the hill. Two land cover change scenarios were compared in the SWAT model. The results showed that there was no significant difference in simulated streamflow during 2004–2015 and sediment loading during 2004–2009; however, the model performed better for sediment loading during 2010–2015 with dynamic land cover change (coefficient of determination (R2) = 0.66, Nash-Sutcliffe efficiency coefficient (NSE) = 0.62, percent bias (PBIAS) = 10.5%, root mean square error observation standard deviation ratio (RSR) = 0.62) than with constant land cover (R2 = 0.61, NSE = 0.54, PBIAS = −17.3%, RSR = 0.68), indicating that long-term land cover change should be considered in hydrologic modeling. Changes in landslides during 2008–2013 were found to significantly affect ecohydrological processes, especially after 2011. In general, annual precipitation plays a dominant role, and landscape composition had by far the strongest influence on water yield and sediment yield compared to landscape configuration. The results can be useful for understanding the effects of land cover change on ecohydrological processes in the Chenyulan watershed and the potential impact of ecohydrological changes on the environment and public health.


Introduction
Land use and land cover are the results of interactions among the natural environment and human activities, and their distribution can reflect anthropogenic types and decision behaviors [1].
( Figure 1d). The average annual precipitation in the Chenyulan watershed is between 2000 and 4000 mm, and approximately 80% of annual rainfall is between May and October (typhoon season). During 2008-2013, there were five to eight typhoons each year, and the most severe were Typhoon Sinlaku (1468.5 mm of accumulated rainfall at Shen-Mu gauge during 11-16 September 2008

Data Source and Description
The key input data for hydrological modeling are the digital elevation model (DEM), soil data, weather data, and land use/cover data. The DEM data are at a 30 m resolution, processed by the Center for GIS, Research Center for Humanities and Social Sciences (RCHSS), Academia Sinica, Taiwan, in 2012. The soil data were collected from the Construction and Planning Agency, Ministry of the Interior, Taiwan. The surveyed soil data contain the soil erodibility factor (USLE-K); hydraulic conductivity; and percentages of silt, clay, and sand, while other information needed by the SWAT model (soil bulk density (SOL_BD), available water capacity of the soil layer (SOL_AWC), and saturated hydraulic conductivity (SOL_K)) were further calculated by using the soil-plant-air-water (SPAW) model developed by [35]. The basic soil parameters for the SWAT model are shown in Table 1.

Data Source and Description
The key input data for hydrological modeling are the digital elevation model (DEM), soil data, weather data, and land use/cover data. The DEM data are at a 30 m resolution, processed by the Center for GIS, Research Center for Humanities and Social Sciences (RCHSS), Academia Sinica, Taiwan, in 2012. The soil data were collected from the Construction and Planning Agency, Ministry of the Interior, Taiwan. The surveyed soil data contain the soil erodibility factor (USLE-K); hydraulic conductivity; and percentages of silt, clay, and sand, while other information needed by the SWAT model (soil bulk density (SOL_BD), available water capacity of the soil layer (SOL_AWC), and saturated hydraulic conductivity (SOL_K)) were further calculated by using the soil-plant-air-water (SPAW) model developed by [35]. The basic soil parameters for the SWAT model are shown in Table 1.  1 The parent material is not seen in this type of soil, which is usually thick-layered and compacted, and has poor tillage due to poor drainage. USLE-K, soil erodibility factor; SOL_BD, soil bulk density; SOL_AWC, available water capacity of soil layer; SOL_K, saturated hydraulic conductivity.
Daily weather parameters (precipitation, minimum and maximum air temperature, relative humidity, solar radiation, and wind speed) were collected from the Data Bank of Atmospheric and Hydrologic Research (DBAR), Taiwan. Daily streamflow and sediment loading were collected from  the hydrological yearbook of Taiwan from 2003 to 2015, published by Figure 2). The SPOT images were purchased from the Space and Remote Sensing Research Center (SRSRC), and further used for watershed land cover classification.

Selection of Ground-Truth Points
In order to select the ground-truth points for land cover classification, we first calculated the normalized difference vegetation index (NDVI) (Equation (1)). NDVI values range from -1 to +1. Negative values are mainly generated from clouds, water, and snow. A zero value means no vegetation (i.e., rocks and bare soil), and very low positive values (0.1 and below) represent barren areas. Moderate values (0.2-0.3) correspond to shrub and grassland, while values close to +1 indicate the highest possible density of green leaves (i.e., forests).
where NIR denotes near infrared reflectance and R denotes red (visible) reflectance.
Based on the NDVI maps (Figure 3), the watershed was divided into two groups: NDVI ≤ 0 and NDVI > 0. The areas of NDVI values smaller than 0 indicate that the possible land cover type is river, built-up, or landslide, while the areas of NDVI values greater than 0 indicate grassland, cultivated land, or forest. We additionally collected the land cover classifications of 1996 and 2005 from [33] and the landslide maps reclassified by satellite images derived from the Forestry Bureau, Council of Agriculture, Executive Yuan (FBCAEY) as reference maps. Areas where specific land classes were unchanged between 1996 and 2005 helped to narrow down the supervised boundary, and the FBCAEY landslide map helped to increase the accuracy of landslide area selection.
A number of points for each land cover type were selected based on the NDVI and reference maps ( Figure 4). For example, 200 points were selected as water within the area of NDVI ≤ 0 in the NDVI map with reference to the water area in the unchanged 1996-2005 map. According to the relative sizes of land cover types, 200 points were selected for built-up, landslide, and grassland, and 1000 and 1500 points were selected for cultivated land and forest, respectively.
Negative values are mainly generated from clouds, water, and snow. A zero value means no vegetation (i.e., rocks and bare soil), and very low positive values (0.1 and below) represent barren areas. Moderate values (0.2-0.3) correspond to shrub and grassland, while values close to +1 indicate the highest possible density of green leaves (i.e., forests).
where NIR denotes near infrared reflectance and R denotes red (visible) reflectance. Based on the NDVI maps (Figure 3), the watershed was divided into two groups: NDVI ≤ 0 and NDVI > 0. The areas of NDVI values smaller than 0 indicate that the possible land cover type is river, built-up, or landslide, while the areas of NDVI values greater than 0 indicate grassland, cultivated land, or forest. We additionally collected the land cover classifications of 1996 and 2005 from [33] and the landslide maps reclassified by satellite images derived from the Forestry Bureau, Council of Agriculture, Executive Yuan (FBCAEY) as reference maps. Areas where specific land classes were unchanged between 1996 and 2005 helped to narrow down the supervised boundary, and the FBCAEY landslide map helped to increase the accuracy of landslide area selection.
A number of points for each land cover type were selected based on the NDVI and reference maps ( Figure 4). For example, 200 points were selected as water within the area of NDVI ≤ 0 in the NDVI map with reference to the water area in the unchanged 1996-2005 map. According to the relative sizes of land cover types, 200 points were selected for built-up, landslide, and grassland, and 1000 and 1500 points were selected for cultivated land and forest, respectively.   In our previous study [36], we proposed a nearest feature line (NFL) embedding transformation for land cover classification. The feature points (prototypes) were manually collected and labelled for classifier training. Since there are very few collected feature points during the training phase, the discriminant power of the trained classifier is decreased in the classification phase. A feature line is

Nearest Feature Line Embedding
In our previous study [36], we proposed a nearest feature line (NFL) embedding transformation for land cover classification. The feature points (prototypes) were manually collected and labelled for classifier training. Since there are very few collected feature points during the training phase, the discriminant power of the trained classifier is decreased in the classification phase. A feature line is generated by two feature points and represents a linear interpolation or extrapolation of each pair of feature points within the same class. An extremely high number of pseudo-prototypes for each class are generated for training by linear interpolation, which enhances the classification performance. The NFL embedding strategy was used to construct the point-to-line adjacency matrix instead of the point-to-point matrix during training. This measurement was directly embedded in the transformation in the discriminant analysis, not in the classification phase. Class separability, neighborhood structure preservation, and nearest feature space (NFS) measurement were all considered to find the most effective and discriminating transformation in the eigenspaces for land cover classification. In this study, the images from 2008 to 2013 consist of four bands: green, red, near infrared reflectance (NIR), and shortwave infrared (SWIR). The criterion used for judging the accuracy of final SPOT images was an overall accuracy value exceeding 70%.

Landscape Metrics
Landscape metrics are usually used to describe the landscape ecosystem, format, and trend of landscape change to analyze the interactions among land uses and anthropogenic activities in watersheds [9,10]. We adopted FRAGSTATS software, developed by the United States Department of Agriculture (USDA), to quantify the composition and spatial configuration of land cover types [34]. Based on previous studies [37][38][39][40][41], we selected a subset of metrics that are commonly used and can affect ecohydrological processes to analyze landscape changes in the Chenyulan watershed from 2008 to 2013, and further evaluate how watershed responses and ecohydrological processes were affected by these changes. Landscape composition was quantified by the proportion of each land cover type. Configuration metrics included: (1) patch-based metrics: patch density (PD) and area-weighted mean patch area (AREA_AM); (2) shape metrics: edge density (ED), area-weighted mean radius of gyration (GYRATE_AM), and area-weighted mean shape index (SHAPE_AM); and (3) aggregation metrics: aggregation index (AI) and splitting index (SPLIT). The criteria for the landscape metrics were suggested by cases. Both PD and SPLIT describe the degree of subdivision of the class or landscape, and can be regarded as the degree of spread of disturbance in this study. AREA_AM, ED, and GYRATE_AM represent the physical continuity of the landscape, and can indirectly explain the influences on ecohydrologcial change. SHAPE measures the complexity of patch shape compared to a standard shape (square) of the same size. Thus, the index equals 1 for square patches of any size. AI refers to the tendency of patch types to be spatially aggregated. Detailed descriptions and equations of landscape metrics can be found in the FRAGSTATS documentation [34].

Model Description
The Soil and Water Assessment Tool (SWAT) was used to evaluate watershed responses to landscape changes induced by natural disturbances and anthropogenic activities in the Chenyulan watershed. The SWAT model was developed by the USDA Agricultural Research Service (USDA-ARS) in 1994, and it can predict long-term impacts of land use management on streamflow, sediment, and nutrient loadings in a watershed at different spatiotemporal scales [19]. In the SWAT model, a watershed is delineated into several subwatersheds, which are further portioned into homogeneous units (hydrologic response units, HRUs) with a unique combination of land use/cover, soil, and slope. For streamflow simulation, the surface runoff volume is computed using a modified Soil Conservation Service (SCS) curve number method [42]. The modified universal soil loss equation (MUSLE) was used to estimate soil loss at the HRUs [43]. More details on the theory can be found in the SWAT 2009 Theoretical Documentation [44].

Land Cover Update Module
In order to incorporate land cover changes during 2008-2013 into the SWAT model, the land cover update (LUP) module in SWAT was activated. Two land cover scenarios were simulated to quantify the impact of land cover changes on water yield and sediment yield. They are constant land cover (CLC), which assumes that land cover remains constant since 2005, and updated land cover (ULC), which represents the dynamic land cover during 2008-2013. To activate the LUP module, two types of files need to be prepared. One is an lup.dat file, which lists the order of changing dates of each land cover, and the other is the HRU fraction (HRU_FR) file of different land covers of concern. For the ULC scenario, the SWAT model starts to read the land cover data on the date when the SPOT image was taken, and stops reading on the previous date before the next SPOT image was taken.

Model Calibration and Validation
The sensitivity analysis, calibration, and validation of the SWAT model were done by using the SWAT Calibration Uncertainty Program (SWAT-CUP), which is open source software developed by [45]. In this study, Sequential Uncertainty Fitting version 2 (SUFI2) was selected for uncertainty analysis. The model performance was evaluated by using four statistical measures: coefficient of determination (R 2 ), Nash-Sutcliffe efficiency coefficient (NSE), percent bias (PBIAS), and root mean square error (RMSE)-observation standard deviation ratio (RSR), as suggested by [46].

Classification Results
The SPOT images were classified into six land cover types (river, grassland, built-up, cultivated land, landslide, and forest) ( Figure 5 and Table 2), with the overall classification accuracy ranging between 73.70% and 83.74% (Table 3). Forest was the major land cover, occupying 74.45-76.75% of the watershed. Cultivated lands are usually developed along streams, with an area between 11.87% and 14.05% of the watershed. However, cultivated lands tended to be smaller and aggregated during the study period. The image classification results showed that river, grassland, and built-up areas did not change much during 2008-2013, with ranges of 2.97-3.63%, 4.34-5.72%, and 0.44-0.76%, respectively. Due to several severe typhoons during 2008-2009 and the cumulative impact of typhoons, landslides increased from 2.00% in 2008 to 2.73-3.11% during 2010-2013.
It should be noted that the land cover data for 2008-2013 were adjusted based on the 2005 HRUs for updated land cover (ULC) scenario modeling. Compared to the original land cover areas, the adjusted land cover during 2008-2013 changed by −0.01 to 2.27 km 2 , −7.94 to 0.13 km 2 , −0.37 to 0.13 km 2 , 1.25-4.14 km 2 , 0.01-0.80 km 2 , and −4.99 to 4.30 km 2 for water, grassland, built-up, cultivated, landslide, and forest, respectively ( Table 2). In particular, grassland in 2005 was not evenly distributed in all subwatersheds. If grassland was identified in a subwatershed in any year during 2008-2013 where there was no grassland HRU in 2005, the grassland area would be replaced by cultivated land and forest. Therefore, the adjusted grassland area was generally smaller than the original area, while cultivated land and forest areas increased slightly after adjustment. distributed in all subwatersheds. If grassland was identified in a subwatershed in any year during 2008-2013 where there was no grassland HRU in 2005, the grassland area would be replaced by cultivated land and forest. Therefore, the adjusted grassland area was generally smaller than the original area, while cultivated land and forest areas increased slightly after adjustment.       Table 4 shows the Pearson's correlation for seven landscape metrics at the landscape level, and their values for 2008-2013 are shown in Table 5. There is a strong positive relationship between PD and ED, with more patches and longer edge lengths, indicating a high level of fragmentation. SHAPE_AM was found to be positively correlated with PD and ED. A high SHAPE_AM value indicates that patch shapes are less compacted. GYRATE is equal to the mean distance between each cell and the centroid of that patch. GYRATE has a zero value when the patch consists of only one pixel, and increases without limit as the patch grows. Therefore, GYRATE_AM is sensitive to patch area (AREA_AM). SPLIT is negatively correlated with GYRATE_AM and AREA_AM, while AI is negatively correlated with PD and SHAPE_AM.
A dramatic change was found in PD, AREA_AM, ED, and GYRATE_AM for the period 2010-2012 (Table 5) This finding is consistent with studies indicating that increased rainfall results in a higher percentage of shrub patches, with a higher shrub density and height [47,48]. Moreover, it was found that the GYRATE_AM value increased during 2008-2013, except in 2011, showing the process of fragmentation of a land cover patch beginning with a reduction in patch area and an increase in the proportion of edge-influenced patch area [49]. Table 4. Correlation matrix of landscape metrics at the landscape level. PD, patch density; AREA_AM, area-weighted mean patch area; ED, edge density; GYRATE_AM, area-weighted mean radius of gyration; SHAPE_AM, area-weighted mean shape index; AI, aggregation index; SPLIT, splitting index.   Table 6 shows the values of class-level landscape metrics and their Pearson's correlations (Table 7). Forest occupies more than 70% of the watershed and usually exhibits as clusters. It is expected that forest has a low PD; the lowest SPLIT; and the highest AREA_AM, ED, GYRATE_AM, SHAPE_AM, and AI. As built-up is the smallest land cover type in the watershed, it has the smallest AREA_AM, ED, and GYRATE_AM. A small AI and high SPLIT indicate that built-up has very low connectivity.

Metrics
Generally, ED is positively correlated with PD for all land cover types, except forest. Similar to what is found in landscape-level results (Table 4), SPLIT is negatively correlated with GYRATE_AM and AREA_AM for all land cover types, as there is a strong positive correlation between GYRATE_AM and AREA_AM (Table 7). Moreover, some relationships are not found at the landscape level, but are significant at the class level, except forest (positive relationship among SHAPE_AM, AREA_AM, and GYRATE_AM, and negative relationship between SPLIT and SHAPE_AM). Forest is the major land cover type of the Chenyulan watershed, thus the uncorrelated relationships between SHAPE_AM and AREA_AM, between SHAPE_AM and GYRATE_AM, and between SPLIT and SHAPE_AM for forest could significantly affect these relationships at the landscape level. However, the relationship between PD and SHAPE_AM is not consistent at the landscape and class levels. While PD was positively correlated with SHAPE_AM at the landscape level, the relationship between PD and SHAPE_AM was positive for river and landslide, and was negative for grassland, built-up, and cultivated land.   [50][51][52] were first examined for the sensitivity analysis. For streamflow parameters, a total of seven parameters with a p-value < 0.05 were selected for calibration. They were curve number (CN2), plant uptake compensation factor (EPCO), surface runoff lag time (SURLAG), baseflow alpha factor (ALPHA_BF), effective hydraulic conductivity in main channel alluvium (CH_K2), and Manning's "n" value for the main channel (CH_N2). In order to reflect the heterogeneity of parameters at different locations in the watershed, some parameters were calibrated separately for head streams (subwatershed nos. 17, 20-23), subwatersheds at a slope greater than 60% (nos. 4,8,12,16,18,19), and downstream subwatersheds (nos. 1-3, 5-7, 9-11, 13-15) ( Figure 6). Table 8 lists the model parameters along with their default values, calibrated ranges, and fitted values. Details of the model parameters and their functions can be found in the SWAT 2012 Input/Output documentation [53].

Model Calibration and Validation
CN2, which governs the surface runoff response, was calibrated for three land covers: forest (FRST), grassland (RNGE), and cultivated land (AGRL). The adjusted CN2 values indicate that the SWAT model with default parameters overestimated the daily streamflow. EPCO ranging from 0.01 to 1.00 was also calibrated in other studies [54,55]. When EPCO = 1.00, the model allows more of the water uptake demand to be met by lower layers in the soil. Therefore, a reduced EPCO indicates that the model allows less variation from the original depth distribution to take place. SURLAG controls the fraction of total available water that will be allowed to enter the reach on any one day [53]. As SURLAG increases, the streamflow hydrograph simulated in the reach will be smoother due to the delay in the release of surface runoff. ALPHA_BF can reflect the groundwater flow response to changes in recharge. A high value of CH_K2 indicates that the bed material has a very high loss rate and the stream is characterized as a flow-through stream that simultaneously receives and loses groundwater. As the head streams have more condense vegetation in the watershed, the CH_N2 value was calibrated higher than other reaches.
The calibration and validation results for daily streamflow showed that there was no significant difference between two land cover scenarios. The model performance was very good in terms of The calibration and validation results for daily streamflow showed that there was no significant difference between two land cover scenarios. The model performance was very good in terms of R 2 = 0.81, NSE = 0.81, PBIAS = −17.3%, and RSR = 0.44 for calibration, and R 2 = 0.71, NSE = 0.7, PBIAS = 0.2%, and RSR = 0.55 for validation (Figure 7).     Two sediment-related parameters, the peak rate adjustment factor (PRF) for sediment routing in the main channel and the linear parameter for calculating the maximum amount of re-entrained sediment in the channel (SPCON), were calibrated for daily sediment prediction (Table 9). Both values increased, so the daily measured versus simulated sediment agreed well (   Two sediment-related parameters, the peak rate adjustment factor (PRF) for sediment routing in the main channel and the linear parameter for calculating the maximum amount of re-entrained sediment in the channel (SPCON), were calibrated for daily sediment prediction (Table 9). Both values increased, so the daily measured versus simulated sediment agreed well (Figure 8

Swat Simulation Results
Annual flow simulation was dominated by rainfall, leading to a similar trend of sediment loadings during 2005-2015 ( Figure 9). For both land cover scenarios, annual flow ranged between approximately 3.3 × 10 11 m 3 and 6.6 × 10 11 m 3 . Annual sediment loading ranged between 246,000 and 589,100 tons for the CLC scenario, and between 245,700 and 589,300 tons for the ULC scenario. The differences in annual flow and sediment between the two scenarios are mainly due to land cover changes and the yields from different land cover types (Table 10)

Swat Simulation Results
Annual flow simulation was dominated by rainfall, leading to a similar trend of sediment loadings during 2005-2015 ( Figure 9). For both land cover scenarios, annual flow ranged between approximately 3.3 × 10 11 m 3 and 6.6 × 10 11 m 3 . Annual sediment loading ranged between 246,000 and 589,100 tons for the CLC scenario, and between 245,700 and 589,300 tons for the ULC scenario. The differences in annual flow and sediment between the two scenarios are mainly due to land cover changes and the yields from different land cover types (Table 10)

Impact of Land Cover Change on Ecohydrological Processes
For the two scenarios, the annual water yields generated from different land covers (cultivated land, landslide, forest, grassland, and built-up) during 2005 and 2008-2014 were compared ( Figure 10). The results showed that water yields were mainly affected by rainfall, so there was a similar trend of water yield generated and a similar average composition of water yield for all land cover types. For the constant land cover scenario, the water yield from cultivated land was composed of 50.92% surface runoff, 14.45% lateral flow, and 34.63% groundwater recharge; from forest, it was 34.69% surface runoff, 32.62% lateral flow, and 32.69% groundwater recharge; and from grassland, it was 44.30% surface runoff, 31.37% lateral flow, and 24.33% groundwater recharge. For the updated land cover scenario, the water yield from cultivated land was composed of 50.84% surface runoff, 13.46% lateral flow, and 35.70% groundwater recharge; from forest, it was 34.67% surface runoff, 32.32% lateral flow, and 33.00% groundwater recharge; and from grassland, it was 43.85% surface runoff, 30.63% lateral flow, and 25.51% groundwater recharge. However, land cover change had a slight impact on water yields generated from landslide and built-up. The composition of the water yield from landslide was 66.93% surface runoff, 13.52% lateral flow, and 19.55% groundwater recharge for the constant land cover scenario, and 67.72% surface runoff, 16.76% lateral flow, and 15.50% groundwater recharge for the updated land cover scenario. The contribution of lateral flow increased and groundwater recharge decreased for the land cover change scenario compared to the constant land cover scenario, indicating increasing pore water pressure, groundwater exfiltration from the bedrock, and hydraulic uplift pressure from below caused by landslides [56]. Similar changes in ecohydrological processes were found for built-up areas. The composition of water yield from built-up areas was 57.56% surface runoff, 4.46% lateral flow, and 37.98% groundwater recharge for the constant land cover scenario and 58.08% surface runoff, 6.74% lateral flow, and 35.17% groundwater recharge for the updated land cover scenario.
The difference in water yield was directly reflected by the change in land cover area during 2008-2014 (Figure 10f-j). Generally, due to decreasing areas of cultivated land, grassland, and built-up areas in the watershed, water yields decreased by 2.7 × 10 9 to 19.6 × 10 9 m 3 , 12.5 × 10 9 to 5.2 × 10 9 m 3 , and 3.5 × 10 9 to 1.8 × 10 9 m 3 , respectively, during 2008-2014 compared to the constant land cover scenario results.
Landslides slightly increased after 2010, resulting in greater water yields for the updated land cover scenario (Figure 10g)  Forest area was the smallest in 2005. Increasing forest area resulted in an increased water yield (Figure 10h). It was also found that the contribution of groundwater increased, indicating continuous groundwater recharge and replenishing rates. Built-up areas primarily consisting of impervious surfaces increase surface runoff and prevent groundwater from recharging to the land. Therefore, decreases in surface runoff (−46.45% to −65.80%) and groundwater recharge (−33.68% to −56.05%) were the two major contributors to the change in built-up water yield between land cover scenarios (Figure 10j). Table 11 shows the interactions between water yield, sediment yield, and landscape metrics at the class level. It should be noted that subwatershed landscape metrics at the class level were individually calculated for each year from 2008-2013, and water yield and sediment yield were simulated under the updated land cover scenario. Because most of the SPOT images were taken in the second half of the year, under simulation, the landscape pattern shows the impact the following year. Thus, simulated watershed responses during 2009-2014 were compared with landscape metrics during 2008-2013. Although the SWAT model is a semi-distributed hydrological model and the simulation results are dominated by landscape composition rather than landscape configuration, the relationship between landscape metrics and watershed responses could be regarded as a cause-and-effect relationship between landscape composition and configuration.

Relationship between Landscape Metrics and Watershed Responses
It was found that all built-up landscape metrics (configuration) were not significantly correlated with water yield and sediment yield, as built-up area is the smallest in the watershed. The patch density (PD) of grassland, cultivated land, and forest was negatively correlated with water yield, indicating that water yield could increase due to increasing greenland in amount and size. In a defined area, a higher patch density means smaller patch areas and more patch edges. The edge density (ED) of grassland, cultivated land, and forest was relatively higher (13.32 m, 21.20 m, 31.46 m, respectively) than the other land covers (built-up 1.81 m, landslide 5.22 m) ( Table 7). The PD and ED analyses proved that the higher the patch and edge density of grassland, the lower the generation of water in cultivated land and forest. This result is in agreement with the finding that increasing agricultural patch density leads to a decrease in surface runoff and sediment yield [17].
Through the fragmentation process, landscape composition and spatial configuration are affected (e.g., patch area, number of patches, patch shape complexity, number of patch edges, distances between patches), resulting in a change in landscape connectivity [49]. Thus, the aggregation index [58] and splitting index [59] are usually used to assess this change. AI reflects the physical aggregation of land covers within a watershed, and a higher value means more aggregation [34]. At the class level, the AI metric was positively correlated with water yield from forest and sediment yield from landslide, while it had a negative relationship with water and sediment yield from cultivated land ( Table 11). As baseflow decreases, an increase in surface runoff is induced by the increasing AI of forest [17]. Moreover, when the landscape is dominated (e.g., >60%) by a given cover type, it is nearly always well connected [60]. The different relationships between the AI metric and water/sediment yield indicate that more water and sediment are found when forest and landslide are more aggregated and cultivated lands are more scattered over the watershed. This finding also indicates that when those landscape areas expand, forest and landslide tend to be more aggregated and cultivated lands tend to be more scattered. Land cover patterns could affect ecohydrological processes and components of the water yield, but also control the amount of water and sediment yield within the watershed. Moreover, the shape metric (SHAPE_AM) of landslide had a positive trend with sediment yield, indicating that a large shape can intensify erosion [39].
The shape indices (ED, GYRATE_AM, and SHAPE_AM) of cultivated lands had negative relationships with water yield and sediment yield, while those indices of landslide had positive relationships with both yields. Usually, the flow rates between land cover types can be enhanced or disrupted, depending on the hardness of the edges [8]. Cultivated land that is disrupted by humans has straight and sharp edges [61], and natural landslides have more curvilinear borders. Thus, the edge characteristics may partially determine the erosion characteristics and sediment export [8]. Moreover, a lower aggregation index (AI) and higher splitting index (SPLIT) resulting in a greater water yield and sediment yield indicate that many small and interspersed land cover patches are more likely to accelerate soil erosion and increase sediment yields [8].

Conclusions
Landscape structures and patterns can affect runoff and non-point source pollution loadings in watersheds [58], and the impact of a single land cover can be affected by other land cover types at the watershed scale [39]. In this study, we first classified SPOT images into six land cover types by using the nearest feature line embedding (NFLE) method, and then quantified the landscape patterns by FRAGSTATS at the landscape and class levels. Our goals were to investigate how landscape patterns affect the water yield and sediment yield from different land cover types, and also to understand how ecohydrological processes are changed when updated land cover change is considered in the SWAT model. The results showed that SWAT could more accurately predict changes in streamflow and sediment exports under the updated land cover scenario. Although the SWAT model is a semi-distributed hydrological model, the relationship between landscape metrics and watershed responses could be regarded as a cause-and-effect relationship between landscape composition and configuration.
The indirect impact of natural disturbances was reflected in the change in landscape configuration, in terms of more fragmentation during 2009-2011 with increasing patch density (PD), edge density (ED), and area-weighted mean shape index (SHAPE_AM). Annual precipitation was the dominant influence on the amount of water yield, while the difference in water yield between land cover scenarios was led by the change in land cover area (landscape composition). The relationships between landscape metrics, water yield, and sediment yield were significant but with a relatively low statistical value, indicating that landscape composition had by far the stronger influence. This finding is in agreement with other studies [62][63][64]. Moreover, the shape indices (ED, GYRATE_AM, and SHAPE_AM) of cultivated lands had a negative relationship with water and sediment yield, while those indices of landslide had a positive relationship with water and sediment yield. However, [39] found that a proper fragile landscape status and more complicated patches can reduce the soil erosion yield and more patch edges can prevent soil erosion by disturbing formation and transportation; and [65] suggested that spatial distribution and the number of farmland areas need to be considered to reduce sediment yields. The contrasting results of previous studies and this study show the complicated relationship between landscape patterns and sediment yield.
By identifying the contributions of different hydrological components to water yield, we can understand how changes in land cover affect ecohydrological processes and how the ecohydrological changes could further affect the environment and public health. Therefore, landscape characteristics that can influence ecohydrological processes should be considered in watershed management, and landscape configuration at the subwatershed level should be considered in the SWAT model. The various interactions between class-level landscape metrics, water yield, and sediment yield are useful for providing guidelines on soil erosion prevention and sustainable hydrologic ecosystem services.