Attribution of Streamflow Variations in Southern Taiwan

Climate change and anthropogenic activity are the main factors impacting the hydrological environment. For sustainable water utilization, identifying the impact contribution of these two factors on the streamflow variations is an important topic in recent research. In this study, seven river basins in southern Taiwan were selected as the study area to evaluate the annual streamflow from 1980 to 2017. The decomposition and elasticity methods based on the Budyko hypothesis were applied to quantify the contribution of climate and anthropogenic factors to the streamflow variations. In addition, the normalized difference vegetation index (NDVI) was used to represent the actual situation of land cover and verify the parameters in the Budyko equation. The two quantitative methods consistently demonstrated that the streamflow variations from pre- to post-period occurred due to the climate factor. The elasticity coefficient of variables demonstrated that the streamflow change is more sensitive to precipitation and this influence reduces from pre- to post-period as the streamflow increase. In the NDVI variations, except for the Yanshui and the Linbain rivers, the Budyko equation parameters changed consistently with NDVI. The present study provides effective results on the contribution of streamflow variations in southern Taiwan to serve as a reference for future water management.


Introduction
The efficient management and protection of precious water resources are needed for the sustainable development of human society [1,2]. However, the global average surface air temperature has increased during the 20th century and the majority of assessments indicate that warming will begin in the future [3,4]. The hydrological cycle pattern demonstrated that dry regions get drier, whereas wet regions get wetter [5]. This change in the natural environment results in the uncertainty of controlling available water resources and, consequently, difficulties in managing water resources occur in many regions. The extreme phenomena, like droughts and floods caused by the anomalous changes in a precipitation pattern, occurred more frequently and dramatically around the world [6,7]. To meet the increasing water resources demand, an increasing number of scholars have proceeded to study this issue, focusing on evaluating water resource availability.
Primarily, among the hydrological variables, the related studies were mainly focused on the streamflow variations, as an intuitive amount reflecting both the climate variations and the patterns of water consumption [8,9]. Analyzing the impacts on the streamflow variations, the studies can reveal the problems in various hydrological systems and guide water resource management. These impacts can be divided into two types: Climate and anthropogenic factors [10]. The climate factor affects the hydrological amounts directly, changing the weather condition and causing the uneven distribution are more unpredictable for demand. According to the forecasting, under different scenarios with Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (AR4), the precipitation and streamflow will gradually decrease in the future [44]. Primarily, this tendency stems from the increase of temperature, water demand, rainfall intensity and non-rainfall days, along with the decline of a reservoir storage capacity. With the inherent disadvantage of climate conditions, southern Taiwan faces more serious water resources problems compared to other areas. Evidently, the change of the hydrological environment in southern Taiwan requires more attention to satisfy the future demand for water resources. Yeh and Tsao [45] have already evaluated the attributions of streamflow variations in southern Taiwan. However, just four river basins were selected as study areas in the previous publication, based on the available hydrological and meteorological datasets from gauging stations. Thus, the streamflow in the study can only represent the relatively small areas of each river basin above certain gauging stations for the southern district. To modify the disadvantage, this study applied a new dataset to present the variations of the main seven river basins and make the results more representative of southern Taiwan. Furthermore, the Budyko approach with the integration of decomposition and elasticity method has been applied in the present study to adopt the advantages and drawbacks of each method. Due to the inclusion of basic physical processes, these two approaches are more reasonable and reliable compared to other statistical methods for determining the streamflow change. The present study aims to evaluate the impacts on water resources in southern Taiwan as a quantifying contribution of climate and anthropogenic factors to the streamflow. The obtained results can serve as a reference for related research on analyzing the hydrological environment of river basins.

Study Area
Southern Taiwan is a part of several administrative districts including Chiayi County, Chiayi City, Tainan City, Kaohsiung City, and Pingtung County with a total area of approximately 10,002 km 2 . The area lies within the tropical monsoon climate zone with annual mean precipitation of 2500 mm and temperature of 24 • C. The weather in wet seasons (from May to October) is influenced by typhoons and southwestern monsoons, resulting in 90% of the total precipitation [43]. This fact indicates that the weather of southern Taiwan is characterized by distinct wet and dry seasons. The disparity of the rainfall makes little water in the river during no-rainfall periods and the water becomes muddy due to a lot of sediment after concentrated rainfall.
Bazhang, Zengwen, Yanshui, Erren, Gaoping, Donggang, and Linbain River Basins were selected as the study areas. These basins are the major seven basins in southern Taiwan from the North to the South. The seven basins have a large terrain variation (average elevation from 34.5 to 1051.8 m) and the drainage areas also show big differences (from 339.2 to 3256.9 km 2 ). Most of the basins belong to sedimentary basins and the permeability of each region depends on different types of deposit. Based on Gaoping River Basin as the boundary, the geology structure of southern Taiwan can roughly be divided into the Chianan Plain and the Pingtung Plain (including Gaoping River Basin). In Chianan Plain, the northern region is the marine and continental deposit formed by sea level change in the past and the southern region shows high spatial heterogeneity due to tectonic movements. The upstream of Pingtung Plain is mainly composed of sandstone, shale, and their interbedding. The middle and downstream rivers are modern alluvial deposits with gravel and sand. The geographic distribution and more detailed information of river basins in the study area are shown in Figure 1 and Table 1.

Data Processing
The precipitation data were obtained from the Taiwan Climate Change Projection Information and Adaptation Knowledge Platform (TCCIP), which utilizes the geoinformation of reference stations and takes into account the climate signal of the target station to establish and complete the  The precipitation data were obtained from the Taiwan Climate Change Projection Information and Adaptation Knowledge Platform (TCCIP), which utilizes the geoinformation of reference stations and takes into account the climate signal of the target station to establish and complete the long-term gridded dataset in Taiwan. TCCIP provided a 5 km × 5 km daily gridded rainfall dataset across Taiwan via a natural neighbor interpolation scheme from 1969 stations. The evapotranspiration and potential evapotranspiration were acquired from the Global Land Evaporation Amsterdam Model (GLEAM). Through a set of algorithms, GLEAM provided datasets cataloging the evapotranspiration, bare-soil evaporation, interception loss, open-water evaporation, and sublimation on a 0.25 • latitude-longitude regular grid with a daily temporal resolution [46,47].
The streamflow data reflected in this study were calculated via the water balance equation [48]. Generally, the water balance equation is expressed as follows: where P, E, and Q are the precipitation, evapotranspiration, and streamflow, respectively, and ∆S represents the change in water storage. ∆S can be taken as zero in the long-term period (10 years or more) [49]. Before applying the data, all variables have been aggregated from daily to yearly values.

Time Series Analysis
The Mann-Kendall test [50,51] and Sen's slope estimator [52,53] have been performed in the preceding operation to understand the trend of each hydrological variable. It is worth noting, that the two change point tests were employed to determine the change points in the streamflow series. In the framework of attribution analysis, the result can provide a rational estimation to divide the data sequence into two unequal periods. The Pettit test [54] can evaluate the probability of the change points occurring at each time. The Mann-Kendall-Sneyers test [55] can approximate the beginning of a change in a time series based on the progressive and retrogressive statistical analysis by the graphical technique [56]. Within this test belonging to the null hypothesis, the intersection of forward and backward sequences inside a confidence interval indicates the beginning of a step change point [57]. After dividing the study period into two sub-periods, the change of streamflow was calculated as the difference of periods. In this study, we determine the total change of streamflow in a watershed as: where Q pre and Q post are the mean streamflow during the pre-and post-period, respectively, and the total change of it (∆Q) is assumed to be the combination of climate change (∆Q C ) and human activity (∆Q H ) [20]. The next sections describe two methods quantifying the impact of both factors on the streamflow.

The Budyko Framework
The relationship between the hydrological variables can be expressed through Equation (1). Budyko [34,35] proposed that the evaporation ratio (the ratio of actual evapotranspiration to precipitation) is a function of the aridity index (the ratio of potential evapotranspiration to precipitation) which can be expressed as: where P, E, and PET are the precipitation, evapotranspiration, and potential evapotranspiration, respectively. Two boundary conditions are attributed to Equation (3): When PET/P→0, E/P→0, and when PET/P→∞, E/P→1. Table 2 outlines the equations derived from other studies on the basis of this hypothesis. Each equation is suitable for various environmental situations. To find an appropriate equation for our study area we screened the best fitting Budyko-type equation as the process shown in Figure 2. Finally, we applied the equation using the following methods to quantify the hydrological response in terms of natural and anthropogenic factors.

Decomposition of the Budyko type curve
Following the movement of data points on the Budyko curves, Wang and Hejazi [15] proposed a method, determining the proportions of streamflow variations influenced by climate and anthropogenic factors. The decomposition method assumes that if a watershed changes without the direct anthropogenic impact, the data point will still move along a single curve due to climate change [65]. Thus, the anthropogenic contribution can be expressed with:

•
Decomposition of the Budyko type curve Following the movement of data points on the Budyko curves, Wang and Hejazi [15] proposed a method, determining the proportions of streamflow variations influenced by climate and anthropogenic factors. The decomposition method assumes that if a watershed changes without the direct anthropogenic impact, the data point will still move along a single curve due to climate change [65]. Thus, the anthropogenic contribution can be expressed with: where P post , E post /P post are the precipitation and evaporative ratio during the post-period, respectively, and E post /P post is the evaporative ratio calculated using the same aridity index (PET post /P post ) in the pre-period. Furthermore, the contribution from climate change (∆Q C ) can be calculated as Equation (2). •

Climate Elasticity Method
The concept of elasticity was proposed by Schaake [38] to evaluate the sensitivity of streamflow to the climate variables. The elasticity coefficient can be defined as the proportional change of streamflow divided by the proportional change of impact factors [66][67][68]: where x is the factor affecting the streamflow. The magnitude of elasticity correlates with the percent change of streamflow as the impact factor increased by 1%. In this study, we determine the change generated by precipitation and potential evapotranspiration as the impacts of climate on the streamflow Water 2020, 12, 2465 7 of 21 (∆Q C ) and the change caused by watershed characteristics as the anthropogenic impacts on streamflow (∆Q H ). After calculating the elasticity of each variable, the impact contribution can be expressed as follows: where Q, P, and PET are the mean annual streamflow, precipitation, and potential evapotranspiration in the whole study period, respectively, and n is a parameter introduced in the parametric Budyko-type equations.

Time Series of Hydrological Data
It follows from the hydrological data collected in this study that the annual mean precipitation in southern Taiwan is about 2200 mm. Owing to the increase of water entering (precipitation), the time series revealed that the streamflow exhibits an increasing trend in seven river basins. Figure 3 shows the annual streamflow trend and the variation of precipitation in each river basin from 1980 to 2017, whereas Table 3 provides the statistical information on each series. The change points test demonstrated that the change points in the Bazhang, Zengwen, Yanshui, Erren, Gaoping, Donggang, and Linbain River basins occurred in 2000, 2003, 2004, 1996, 2001, 1996, and 2004, respectively. To examine the existence of change points in the streamflow sequence, all years have been divided into three categories according to the streamflow amount. Specifically, we defined the largest 30% as the wet years, the smallest 30% as dry ones, and the rest as the normal years. Figure 4 demonstrates the distribution of years. It follows from Figure 4 that transitions occurred with a greater magnitude before and after the change points in each river basin. The variations of the streamflow between two periods in each river basin are given in Table 4. The following attribution analysis is focused on these values.

Fit the Best Budyko-Type Equations
Following the processes illustrated in Figure 2, we used the equations given in Table 2 to simulate the streamflow and compare the obtained simulations with several criteria. The performance of each equation is shown in Table 5. The estimation pattern of the observation and each equation including correlation, root-mean-square deviation (RMSD), and amplitude of its deviations, simultaneously were quantified through the Taylor diagram [69] shown in Figure 5. As a single parameter represents the watershed characteristics, the single-parameter Budyko-type equations exhibited a better simulation ability than non-parameter ones. Overall, the equation proposed by Wang and Tang [36] exhibited the best performance among the eight Budyko-type equations with the R 2 and MAE of 0.992 and 22.9 mm, respectively. The elasticity coefficient will be derived from the Wang-Tang equation in the next paragraph.

Sensitivity Analysis with the Elasticity Coefficient
The results presented in the last section demonstrated that the Wang-Tang Budyko-type equation is applicable for southern Taiwan and the elasticity coefficients can be subsequently derived. In a long-time period, the mean annual streamflow can be derived from Wang and Tang [36] and Equation (1) as follows: where Q, P, and PET are the streamflow, precipitation, and potential evapotranspiration, respectively, and n is the parameter in the Budyko-type equation. Applying Equation (5), the elasticity coefficients of three hydrological variables ( P ε ) can be expressed as:

Sensitivity Analysis with the Elasticity Coefficient
The results presented in the last section demonstrated that the Wang-Tang Budyko-type equation is applicable for southern Taiwan and the elasticity coefficients can be subsequently derived. In a long-time period, the mean annual streamflow can be derived from Wang and Tang [36] and Equation (1) as follows: where Q, P, and PET are the streamflow, precipitation, and potential evapotranspiration, respectively, and n is the parameter in the Budyko-type equation. Applying Equation (5), the elasticity coefficients of three hydrological variables (ε P ) can be expressed as: Table 6 summarizes the elasticity coefficients calculated via Equations (9)- (14). The overall impacts of three hydrological variables in every river basin suggest that the precipitation causes positive effects whereas the other factors affect the streamflow negatively. Take the annual streamflow variations in the Bazhang River basin as an example. When the precipitation, potential evapotranspiration, and parameter n increase by 1%, the streamflow changes by 1.98%, −0.98%, and −0.49%, respectively. The greatest, intermediate, and the smallest absolute values of the coefficients correspond to the precipitation, potential evapotranspiration, and the parameter n, respectively. Consequently, the change in streamflow in southern Taiwan was more sensitive to the precipitation compared to the other two variables. It is worth noting that the variation of elasticities indicates that the influence of these three variables reduces with the increasing trend of streamflow. The magnitude of changes depended on the relative variations between streamflow and variables in each river basin.

Attribution Analysis
After understanding the change characters and sensitivity of hydrological variables analyzed in previous paragraphs, the climate elasticity and decomposition methods were applied to quantify the contribution of climate and anthropogenic factors impacting the streamflow. Figure 6 shows variations of Budyko-type curves from pre-period to post-period for each river basin fitted by the Wang-Tang equation. As seen from Figure 6, all parameters decrease in the studied river basins except for the Erren and Linbain River Basins. This result indicates that the watershed characteristics change Water 2020, 12, 2465 13 of 21 differently than others and the anthropogenic factor declines the streamflow in these two basins. The quantitative results are shown in Figure 7. The results obtained by two methods used in this study consistently demonstrated that the climate is the major influential factor in each river basin. Specifically, the streamflow variations given in Table 4 are primarily originated from the climate factor. The results presented in this section illustrate the relative impacts of climate and anthropogenic factors on the streamflow variations for reference to water management.

Watershed Parameters Versus NDVI Variations
The n parameter in the Budyko fitting curves, can represent the watershed characteristics, such as vegetation, soil properties, and topography to illustrate the environment condition in the river basins. According to the previous studies, the magnitude of parameter n increases as the vegetation cover becomes abundant [40,70]. Thus, to verify the condition of each river basin mentioned in the previous section, we applied the actual vegetation changes (e.g., NDVI) reflecting the observed

Watershed Parameters Versus NDVI Variations
The n parameter in the Budyko fitting curves, can represent the watershed characteristics, such as vegetation, soil properties, and topography to illustrate the environment condition in the river basins. According to the previous studies, the magnitude of parameter n increases as the vegetation cover becomes abundant [40,70]. Thus, to verify the condition of each river basin mentioned in the previous section, we applied the actual vegetation changes (e.g., NDVI) reflecting the observed variations of n in the Wang-Tang equation. The NDVI (Normalized Difference Vegetation Index) in each river basin is computed in ArcGIS from the satellite images of Landsat5 and Landsat8 in United States Geological Survey (USGS) EarthExplorer. Figure 8 shows the variations of normalized difference vegetation index (NDVI) and the differences from pre-to post-period are presented in Table 7. Compared to the change of parameter n, except for the Yanshui and Linbain River Basins, the NDVI of other river basins manifested the consistent variations from pre-to post-period (see Table 8). Namely, the n in these river basins can clearly display the specific physical mechanisms to describe the actual environmental changes.

Basic Assumptions in This Study
The evaluation presented in this study can easily and efficiently present the attribution of streamflow variations. Nevertheless, some applications or assumptions in this study must be emphasized. First, the water balance equation we used to calculate the streamflow was based on the assumption that the water storage change can be neglected in the long-term period. This approach is not applicable to some specific areas with certain geological conditions. However, this application was also found in previous research for evaluating spatial and temporal differences of each hydrological variable [48]. The technique can not only ensure the hydrological variables in the long-term period satisfy the assumptions according to Budyko framework, but also can conquer the limitation of the gauging stations. Second, the environment of the pre-period was accepted as natural without any influence. If various impacting factors are needed to be considered, other methods have to be applied. Third, the impacts of climate and anthropogenic factors were defined through the variations of climate variables and the parameter n in the Budyko equation, respectively. To evaluate the consolidated impacts of these factors, more efficient methods are needed to be employed in further discussion. Delmonte Oliveira et al. [19] proposed a new hydrological approach in the decomposition method to separate the climate impacts from the variations of parameter n. This method can provide more precise quantitative results compared with the traditional ones.

Using the NDVI Variations to Verify Watershed Parameters
The utilization of the normalized difference vegetation index (NDVI) to confirm the change of parameter n is a practicable operation in previous studies [71,72]. Theoretically, the same pattern of change will be observed in the parameter n and the NDVI from the pre-to post-period [70]. However, in the Yanshui and Linbain River Basins, the changes of n differed from NDVI variations. According to the NDVI differences of basins given in Table 8, the change of NDVI magnitude lay between −0.1 and 0.1, indicating that the NDVI changed less between two periods. It follows from the evaluation mentioned above that, most likely, the inconsistent changes that occurred in the Yanshui and Linbain River Basins stem from the insignificant change of NDVI values. Furthermore, the parameter n can also reflect the soil, topography, or other watershed characteristics [73]. The inconsistent changes in these two river basins can probably be influenced by other factors as well. Thus, other related indexes are encouraged to substitute the NDVI for various comparisons in the future.

Compared with the Previous Publication
The results of this study were shown to be similar with a previous publication [45], that streamflow variations are mainly due to the climate factor. However, this research extended the study area to almost the whole of southern Taiwan and screened the best fitting Budyko-type equation with given processes. These techniques can not only make the results more representative of southern Taiwan, but also improve the simulation of quantitative results. Furthermore, the usage of the two change point tests and the comparison of each elasticity from pre-to post-period can also assist to present a better discussion. Through the improvement and novelty, as mentioned, this study can provide attribution analysis of the streamflow variations in southern Taiwan more clearly than the previous publication.

Conclusions
The hydrological data of southern Taiwan from 1970 to 2017 was collected in this study to evaluate the impact contribution of climate and anthropogenic factors on the streamflow variations. The time series analysis revealed the increasing trend of streamflow in each basin. Based on two statistical methods, the change points of streamflow series in the Bazhang, Zengwen, Yanshui, Erren, Gaoping, Donggang, and Linbain River basins occurred in 2000, 2003, 2004, 1996, 2001, 1996, and 2004, respectively. The simulation processes demonstrated that the Wang-Tang equation is appropriate for southern Taiwan. The elasticity of each variable shows that the streamflow variations are more sensitive to the change of precipitation followed by potential evapotranspiration and the parameter n. Furthermore, the absolute value of elasticity coefficients decreased from pre-to post-period and exhibited that the influence of the three variables reduced as the trend of streamflow series increased. As for the attribution analysis, the consistent results derived by the climate elasticity and decomposition methods revealed that climate is the main influential factor in each river basin. The n parameter in the Budyko fitting curves decreased from pre-to post-period in most of the river basins except for the Erren and Linbain River Basins. Compared to the change of NDVI, except for the Yanshui and Linbain River Basins, the parameter n of other river basins exhibited the consistent variations from pre-to post-period. Consequently, the parameter n in the Budyko equation can clearly display the watershed characteristics to describe the actual environmental changes.