Impacts of Climate Change and Land Use / Cover Change on Streamﬂow in Beichuan River Basin in Qinghai Province, China

: Climate change (CC) and land use / cover change (LUCC) are the main drivers of streamﬂow change. In this study, the e ﬀ ects of CC and LUCC on streamﬂow regime as well as their spatial variability were examined by using the Distributed Hydrology Soil Vegetation Model (DHSVM) for the Beichuan River Basin in the northeast Tibetan Plateau. The results showed that CC increased annual and maximum streamﬂow in the upstream but decreased them in the downstream. CC also enhanced minimum streamﬂow in the whole river basin and advanced the occurrence of daily minimum streamﬂow. Temperature change exerted greater inﬂuence on streamﬂow regime than wind speed change did in most situations, but the impact of wind speed on streamﬂow reﬂected the characteristics of accumulative e ﬀ ects, which may require more attention in future, especially in large river basins. As for LUCC, cropland expansion and reservoir operation were the primary reasons for streamﬂow reduction. Cropland expansion contributed more to annual mean streamﬂow change, whereas reservoir operation greatly altered monthly streamﬂow pattern and extreme streamﬂow. Reservoir regulation also postponed the timing of minimum streamﬂow and extended durations of average, high, and low streamﬂow. Spatially, CC and LUCC played predominant roles in the upstream and the downstream, respectively.


Introduction
Streamflow is one of the key components of hydrologic cycle, an important water resource for maintaining sustainable social and economic development and an ecosystem with diverse life forms [1,2]. However, with the increase of atmospheric concentration of radiatively active gases [3], such as carbon dioxide, and continuous development of human society, global natural water balance has been affected [4]. The former, that is, climate change (CC), affects the hydrologic cycle, such as in precipitation and evaporation [5]. The latter primarily refers to land use/cover change (LUCC) caused by human activities, such as deforestation, urbanization, cropland expansion, and reservoir operation, which affect landscape morphological and physiological properties by altering albedo, leaf area index (LAI), surface roughness, soil infiltration rate, and interruption of stream channels, among others [6,7].
The diversity makes BRB a perfect study area for quantifying the heterogeneous effects of various factors on streamflow.
The area of the BRB is about 2800 km 2 , with the outlet at the Qiaotou station ( Figure 1). The elevation ranges from 2400 to 4800 m and drops off from the northwest to the southeast. Its climate is dominated by East Asian monsoon as well as plateau monsoon [39]. The BRB mean annual air temperature is around 2.8 • C, and the total annual precipitation is approximately 508 mm. Beichuan river originates from the south slope of Daban Mountain of the Qilian Mountain Range (Figure 1). It is a major branch of the Huangshui River, one large tributary of the Yellow River. Baoku River and Heilin River are the two tributaries of the Beichuan River. The BRB supplies 70% of the fresh water consumption (about 1.35 × 10 8 m 3 domestic and industrial use) to Xining City, the capital city of Qinghai Province, with a population of more than 2 million people, being the largest city on the Tibetan Plateau. There are about 10,000 people living in the BRB who herd cattle and grow wheat, canola, and potatoes. The basin also provides irrigation water to the downstream agriculture field, a main crop production region for the Qinghai Province.
Water 2020, 12, x FOR PEER REVIEW 3 of 26 The area of the BRB is about 2800 km 2 , with the outlet at the Qiaotou station ( Figure 1). The elevation ranges from 2400 to 4800 m and drops off from the northwest to the southeast. Its climate is dominated by East Asian monsoon as well as plateau monsoon [39]. The BRB mean annual air temperature is around 2.8 °C, and the total annual precipitation is approximately 508 mm. Beichuan river originates from the south slope of Daban Mountain of the Qilian Mountain Range (Figure 1). It is a major branch of the Huangshui River, one large tributary of the Yellow River. Baoku River and Heilin River are the two tributaries of the Beichuan River. The BRB supplies 70% of the fresh water consumption (about 1.35 × 10 8 m 3 domestic and industrial use) to Xining City, the capital city of Qinghai Province, with a population of more than 2 million people, being the largest city on the Tibetan Plateau. There are about 10,000 people living in the BRB who herd cattle and grow wheat, canola, and potatoes. The basin also provides irrigation water to the downstream agriculture field, a main crop production region for the Qinghai Province.

Model
The Distributed Hydrology Soil Vegetation Model (DHSVM) [40,41] was used in this study. The DHSVM is a fully distributed, physically-based hydrologic model, and has been identified as one of the most suitable models for hydrologic study in complex mountainous terrain among 30 hydrological models [42]. It is typically applied at high resolutions varying from 25 to 200 m at a subdaily time scale for watersheds of up to 10,000 km 2 . The DHSVM has been widely applied in both tropical [43][44][45] and temperate catchments [10,46] in many research fields, including forest-snow interactions [47], climate change [48], landscape change [27,49], human activities [27], urbanization [10,44], stream temperature [50], water quality simulations [51], and sediment transportation [52].
Using a grid, the model explicitly accounts for spatial variations of topography, meteorological elements, soil, and vegetation in a basin, and describes physical processes in energy and water balances. Built on digital elevation model (DEM)-derived grid, each cell with specified resolution has explicit elevation, soil and vegetation types, soil depth, and meteorological variables interpolated

Model
The Distributed Hydrology Soil Vegetation Model (DHSVM) [40,41] was used in this study. The DHSVM is a fully distributed, physically-based hydrologic model, and has been identified as one of the most suitable models for hydrologic study in complex mountainous terrain among 30 hydrological models [42]. It is typically applied at high resolutions varying from 25 to 200 m at a sub-daily time scale for watersheds of up to 10,000 km 2 . The DHSVM has been widely applied in both tropical [43][44][45] and temperate catchments [10,46] in many research fields, including forest-snow interactions [47], climate change [48], landscape change [27,49], human activities [27], urbanization [10,44], stream temperature [50], water quality simulations [51], and sediment transportation [52].
Using a grid, the model explicitly accounts for spatial variations of topography, meteorological elements, soil, and vegetation in a basin, and describes physical processes in energy and water balances. Built on digital elevation model (DEM)-derived grid, each cell with specified resolution has explicit elevation, soil and vegetation types, soil depth, and meteorological variables interpolated from nearby stations. All energy and water balance variables are calculated for each cell. Stream channels are also distributed on the grid on the basis of their actual existence and are connected on the basis of the elevation. Water and energy balances are computed using physical laws, for instance, potential evapotranspiration is calculated by the Penman-Monteith equation for both overstory and understory vegetation. Surface runoff and sub-surface flow are modeled cell-by-cell and first accumulated to the streams, being then routed downstream using the linear reservoir method on the basis of elevation. Vertical unsaturated moisture movement through the soil layers is calculated by using Darcy's law. Canopy snow interception and release is modeled using a one-layer mass-energy balance model, whereas snow accumulation and melt below the canopy are simulated using a two-layer mass-energy balance model [40]. Given rugged terrain and complicated physical processes in the study area, we hereby chose the DHSVM so that high spatial resolution could be fully implemented. Figure S1 shows conceptual representation and flow chart of the DHSVM.
The spatial resolution of the DHSVM in BRB was 100 m. This 100 m resolution was used because the original land use/cover data we obtained was in 100 m resolution, and 100 m resolution can sufficiently reflect terrain in this area. The 100 m resolution also balanced the physical representation of the basin and the efficiency of the performance of the computer. Stream network and soil depth were generated on the basis of the DEM for the BRB. The model was run at 3 hour time step, and the simulation results were aggregated to daily and monthly so that calibration and validation could be performed at daily and monthly steps, as there was only observed daily streamflow available. According to Cuo et al. [46], sensitive parameters of the DHSVM are lateral saturated hydraulic conductivity, exponential decrease in lateral hydraulic conductivity with depth, porosity, LAI, vegetation height, and vapor pressure deficit, which were calibrated in this study. The model warming-up period was 1965, and the analysis was conducted from 1 January 1966 to 31 December 2016 after model calibration and validation.
Model performance was evaluated using Nash-Sutcliffe efficiency (E ns ), the ratio of root mean square error and standard deviation of the observation (RSR), and relative error (E r ). E ns is a normalized statistic that compares the residual variance to observed data variance and is regarded as the best objective function to reflect the overall fit between simulation and observation (see Equation (1)) [53,54]. RSR is normalized root mean square error (RMSE; Equation (2)), and it considers both an error index statistic and a scaling factor, making it applicable to multifarious constituents [53]. E r describes the difference in the means between simulation and observation (Equation (3)).

Data
There are four hydrological stations and one meteorological station in the BRB ( Figure 1). Datong (DT) had daily temperature maxima and minima, precipitation, and wind speed operated by Qinghai Meteorological Bureau. Daily streamflow and precipitation were collected from four hydrological stations: Niuchang (NC) above Heiquan dam, Xiamen (XM) on Baoku River, Heilin (HL) on Heilin River, and Qiaotou (QT) on Beichuan River, all managed by the Qinghai Bureau of Hydrology and Water Resources Survey. The station operation periods were different (Table 1). XM was removed in 2001 when Heiquan reservoir was put into use, and NC was built in 2002.
In light of the reality that no meteorological observations were available in the altitude higher than 2459 m in the BRB, we installed an automatic weather station (AWS) in September 2013, situated at 3463 m to measure meteorological elements at hourly time steps. Due to the short periods, the meteorological data observed at AWS could not be used as model forcing. Instead, we used daily mean temperature at this station and DT to compute temperature lapse rate in the BRB for the overlapping period in 2014 and 2017, obtaining −0.81 • C/100 m as the BRB temperature lapse rate, which was then used to interpolate maximum and minimum temperature across the basin. The average of daily maximum and minimum temperature was used as daily mean temperature. We assumed that every station had the same daily wind speed as DT. The missing precipitation was filled by applying the linear relationship for monthly precipitation between the two stations in the overlapping period.
The method of Nijssen et al. [55] was adopted to calculate relative humidity, downward shortwave, and longwave radiation on the basis of daily precipitation and temperature. An iterative algorithm between dewpoint temperature and radiation was used to calculate both variables involving factors of temperature range, minimum temperature, and annual precipitation [56,57]. The estimated dewpoint temperature was used to calculate actual vapor pressure, and with saturation vapor pressure, relative humidity was calculated. Spline interpolation was used to disaggregate daily to hourly temperature and relative humidity. Daily precipitation was evenly distributed to hourly time steps. Daily wind speed was regarded as the same as hourly interval wind. Hourly downward shortwave and longwave radiation were calculated from disaggregated hourly temperature and humidity by taking various factors into consideration, such as location, local time, transmissivity, emissivity, and cloudiness. Hourly forcings were then aggregated or averaged to obtain three-hourly forcing to drive the model. These station data were gridded to 100 m spatial resolution using the Cressman interpolation approach, which calculated weighted average values from multiple stations on the basis of the distance as a radius from the stations.
The DEM was extracted from Shuttle Radar Topographic Mission (SRTM) data (http://srtm. csi.cgiar.org/srtmdata/). The BRB soil class map was downloaded from Harmonized World Soil Database (HWSD) (http://webarchive.iiasa.ac.at/Research/LUC/External-World-soildatabase/HWSD_ Data/), which showed that loam and sandy loam were the main soil types occupying nearly 90% of the basin, being consistent to our laboratory measurement of field-collected soil samples. Land use/cover data for three periods: 1980, 2000, and 2015, were derived from Resource and Environment Data Cloud Platform (http://www.resdc.cn/Default.aspx). Maps of 1980 and 2015 were the earliest and latest available for the BRB, respectively, and the 2000 map was the year before the reservoir operation. Because 1980 and 1990 maps were very similar, we did not use the 1990 map. As there were no existing land cover products before 1980, we created land cover of 1960 on the basis of the land cover map of 1980. According to the older generation who have lived in this area since the 1950s, cropland appeared in the late 1960s. We then constructed land cover for 1960 with assumptions that (1) cropland area in 1980 transformed into grassland, as it was the dominant land cover type of the whole basin, and (2) other land cover types were not changed during 1960-1980.

Statistical Analysis
Investigating long-term trends of hydroclimatic variables improve the understanding of hydroclimatic conditions and changes [58]. In this study, Sen's slope estimator, Mann-Kendall (MK) test [59,60], and moving t-test technique were employed to examine the trends and change-points in various hydroclimatic variables in the BRB. Details of these techniques are provided in the Supplementary Materials section. Table 1 shows mean annual, annual trend, and variability of observed air temperature, wind speed, precipitation, and streamflow obtained from daily time series at each station. Figure 2 shows mean monthly and monthly trend of these variables. The mean annual air temperature at DT station was 4.84 • C (Table 1), high in July and low in January (Figure 2a). Both annual and monthly temperature increased significantly, ranging from 0.056 to 0.118 • C year −1 for the monthly and 0.08 • C year −1 for the annual over 1966-2016. The coefficient of variation (CV yr ) of annual temperature was 26.49%. The average annual wind speed was 1.59 m/s and decreased significantly at −0.008 m s −1 per year during the period. The CV yr of annual wind speed was 22.56%. Wind speed showed a strong seasonal pattern as well, high in spring and low in summer ( Figure 2b). Wind speed significantly decreased in January-May and November-December, and the rates ranged from −0.018 to −0.011, increasing insignificantly in July-September.

Observation
Due to monsoon influence, precipitation was concentrated in May-September, accounting for 80% of annual total, and peaked in August for all stations (Figure 2c-g). Precipitation did not show a significant change, except in February and June at XM and February at DT. Table 1 shows that annual precipitation and elevation had a generally positive relationship in BRB, with the exception of HL. Although HL was located at the second highest place, lower than NC, HL measured the highest annual precipitation (600.88 mm) among all stations, perhaps due to it having the closest distance to Qinghai Lake, the largest inland lake in China. Abundant water vapor from the lake and orographic lifting could favor precipitation in winter and spring. CV yr of annual precipitation was between 11.19% and 14.98%.    The monthly streamflow pattern (Figure 2h-k) followed that of temperature and precipitation in general, with its peak 1 month lag behind that of precipitation at NC, HL, and QT. There were a few significant monthly streamflow trends at NC, HL, and QT, which were different from monthly precipitation change. At NC, HL, and XM, annual streamflow trends followed those of annual precipitation without any statistical significance (Table 1). However, at QT, annual streamflow decreased while annual precipitation increased slightly due to reservoir operation ( Table 1). CVyr of annual streamflow in the BRB was less than 27%, which is relatively small compared to other rivers in China. Figure S2 displays the change point detection of annual data. Figure S2a demonstrates that an obvious upward trend of temperature started in 1993 when the UF (Test statistic, details are provided in the Supplementary Materials section) curve exceeded the 0.05 confidence interval. However, because the intersection of the curves was located outside the confidence interval, the MK method failed to detect the exact year of the abrupt change point. We then used the moving t-test technique to detect again. Figure S2b shows that abrupt changes appeared in 1982 and 1984 for n = 15 and n = 10, respectively. We chose 1983 as the change point, and the annual temperature showed upward trends of 0.012 °C year −1 and 0.113 °C year −1 before and after 1983 on the basis of Sen's slope test. For wind speed, UF and UB (Test statistic, details are provided in the Supplementary Materials section) curves intersected in 1980 ( Figure S2c), an abrupt change year. After 1993, the UF curve was constantly lower than the confidence interval, indicating wind speed consistently and significantly decreasing at p < 0.05. The MK change point test identified no mutation point for precipitation and streamflow for all stations (not shown). In the following analysis, we used temperature and wind speed change to represent climate change during 1966-2016 in the BRB. The monthly streamflow pattern (Figure 2h-k) followed that of temperature and precipitation in general, with its peak 1 month lag behind that of precipitation at NC, HL, and QT. There were a few significant monthly streamflow trends at NC, HL, and QT, which were different from monthly precipitation change. At NC, HL, and XM, annual streamflow trends followed those of annual precipitation without any statistical significance (Table 1). However, at QT, annual streamflow decreased while annual precipitation increased slightly due to reservoir operation (Table 1). CV yr of annual streamflow in the BRB was less than 27%, which is relatively small compared to other rivers in China. Figure S2 displays the change point detection of annual data. Figure S2a demonstrates that an obvious upward trend of temperature started in 1993 when the UF (Test statistic, details are provided in the Supplementary Materials section) curve exceeded the 0.05 confidence interval. However, because the intersection of the curves was located outside the confidence interval, the MK method failed to detect the exact year of the abrupt change point. We then used the moving t-test technique to detect again. Figure S2b shows that abrupt changes appeared in 1982 and 1984 for n = 15 and n = 10, respectively. We chose 1983 as the change point, and the annual temperature showed upward trends of 0.012 • C year −1 and 0.113 • C year −1 before and after 1983 on the basis of Sen's slope test. For wind speed, UF and UB (Test statistic, details are provided in the Supplementary Materials section) curves intersected in 1980 ( Figure S2c), an abrupt change year. After 1993, the UF curve was constantly lower than the confidence interval, indicating wind speed consistently and significantly decreasing at p < 0.05. The MK change point test identified no mutation point for precipitation and streamflow for all stations (not shown). In the following analysis, we used temperature and wind speed change to represent climate change during 1966-2016 in the BRB.

LUCC
As shown in Figure 3 and Table 2, the main landscape types in BRB were grassland, closed shrubland, and cropland, together accounting for more than 80% of the total area throughout the period. Grassland took up about 50% and was distributed across the whole river basin. More than 20% of the basin was occupied by closed shrubland, concentrated mostly in the upstream and midstream areas. The third largest land cover type was cropland, which made up 11.64%-12.75% of the downstream and lowland areas of the basin in 1980, 2000, and 2015. There was no cropland in 1960.
Water 2020, 12, x FOR PEER REVIEW 8 of 26

LUCC
As shown in Figure 3 and Table 2, the main landscape types in BRB were grassland, closed shrubland, and cropland, together accounting for more than 80% of the total area throughout the period. Grassland took up about 50% and was distributed across the whole river basin. More than 20% of the basin was occupied by closed shrubland, concentrated mostly in the upstream and midstream areas. The third largest land cover type was cropland, which made up 11.64%-12.75% of the downstream and lowland areas of the basin in 1980, 2000, and 2015. There was no cropland in 1960.  Table S1 illustrates that from 1980 to 2000, the most significant features of LUCC were urban and built area expansion and interconversion between grassland and other land cover types. Urban and built area increased by 1.57 times and reached to 41.65 km 2 , of which 37% (15.44 km 2 ) was from cropland. Grassland interconverted mostly with closed shrubland, cropland, and bare ground. As a result, grassland increased by 45.52 km 2 , whereas closed shrubland, cropland, wooded grassland, woodland, and bare ground decreased moderately. Cropland decreased due to its conversions to grassland (38.03 km 2 ), closed shrubland (5.76 km 2 ), and woodland (5.06 km 2 ). Water increased by 19.41% (0.64 km 2 ) from cropland and closed shrubland transformation. The transmission matrixes of land use in Tables S1 and S2 illustrate the interconversion between different landscape types from 1980 to 2000 and from 2000 to 2015. From 1960 to 1980, mere grassland transformed to cropland (358.52 km 2 ). Table S1 illustrates that from 1980 to 2000, the most significant features of LUCC were urban and built area expansion and interconversion between grassland and other land cover types. Urban and built area increased by 1.57 times and reached to 41.65 km 2 , of which 37% (15.44 km 2 ) was from cropland. Grassland interconverted mostly with closed shrubland, cropland, and bare ground. As a result, grassland increased by 45.52 km 2 , whereas closed shrubland, cropland, wooded grassland, woodland, and bare ground decreased moderately. Cropland decreased due to its conversions to grassland (38.03 km 2 ), closed shrubland (5.76 km 2 ), and woodland (5.06 km 2 ). Water increased by 19.41% (0.64 km 2 ) from cropland and closed shrubland transformation.  (Table S2) due to the nature reserve that was approved for establishment by the Qinghai People's Government in 2005, being promoted to a national nature reserve by the State Council in 2013, as well as other conservation policies. The greatest change was from bare ground to grassland, which was only 6.61 km 2 . Urban and built area only increased by 3.18 km 2 during this period, indicating limited human activity after 2000. Moreover, cropland decreased slightly, owing to the increment of urban/built area and water area. Water area expanded by two times and was increased from 3.92 km 2 in 2000 to 8.11 km 2 in 2015 due to the operation of Heiquan reservoir since 2001.

Calibration and Validation Results
Model calibration was conducted for 1971-1980 at XM station. HL (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988), XM (1991XM ( -2000, and NC (2002-2016) were used for model validation ( Figure 4 and Table 3). The E ns for monthly streamflow was 0.76, and monthly RSR and E r were low, which indicated that the DHSVM-simulated streamflow matched the "very good" observation during calibration periods (Table 3). Peakflows were overestimated in 1979 (Figure 4a). To further investigate model performance of extreme streamflow, Table 4 compares the standard deviation and relative error of peakflows and low flows selected from the top 10 maximum flow and bottom 10 minimum flow in each year during the calibration period. It can be seen that the model could capture the basic features of peakflows, but there was large relative error (48.64% due to low observed value as denominator) for low flow simulation, although the absolute difference (only at 0.57 m 3 s −1 ) was reasonable.   For model validation, monthly , RSR, and showed "good" or "very good" agreements between simulation and observation, wherein NC station had the highest monthly at up to 0.8  For model validation, monthly E ns , RSR, and E r showed "good" or "very good" agreements between simulation and observation, wherein NC station had the highest monthly E ns at up to 0.8 (Table 3). In contrast, the simulation at HL station was inferior to XM and NC stations; this may have been caused by its low annual streamflow, which was only about 33% and 25% of NC and XM annual streamflow, respectively. Simulated peakflows could always meet the "satisfactory" standard, especially upstream at HL and NC, whereas low flows showed "dissatisfactory" status, owing to the small value (Table 4). In general, the above analysis indicates that the DHSVM model could reproduce the flow pattern fairly well in the BRB.
After calibration and validation, model uncertainty resulting from the parameter uncertainty was evaluated. Sensitive parameters including lateral conductivity, exponential decrease, porosity, field capacity, leaf area index, and vapor pressure deficit were increased (or decreased) by 10% at the same time, and two extra simulations were conducted for the ±10% changes. The streamflow of the two simulations were also shown in Figure 4. The figure shows that model parameter uncertainty yielded very small uncertainty ranges in streamflow, as uncertainty ranges of annual and monthly streamflow of all stations in BRB were within ±5% and ±12%, respectively.
To further evaluate the performance of the DHSVM, we also compared the annual and monthly simulated and observed streamflow regime indices in the longer overlapping periods at NC and HL (Table 5 and Figure S3), where human activity is not intense. Table 5 compares the mean annual; annual variability; and the trends of annual streamflow, annual 1-day peakflow (FX1day), and 1-day lowflow (FN1day). The relative errors of mean annual streamflow were only −8.21% and −10.77% at NC and HL, respectively. The simulated annual variability could also reflect the observed variability, especially at HL. At both stations, the directions of trends of annual mean and extreme streamflow were the same for the simulation and observation, although the quantities were slightly different.  Figure S3 shows the simulated and observed means and trends of monthly streamflow. In general, simulated streamflow could reflect the intra-annual variability to some extent, but with low spring flow and early occurrence of maximum monthly flow, as shown by the lines. The bar graph shows that the directions of most of the monthly trends of simulation and observation were the same at the two stations, and all six opposite trends (January, February, April, and May for NC, July and September for HL) in the simulation and observation were insignificant. Overall, the DHSVM could reflect real hydrological conditions and changes reasonably well in the upper BRB with limited human activity. Next, we investigated CC and LUCC impacts on streamflow regime and the spatial variation of the impacts across the basin using DHSVM throughout 1966-2016, which is longer than the observation periods at most stations except QT.

Scenarios and Streamflow Regime Metrics
As annual precipitation at all stations did not show statistical significant change over the period (Table 1), we did not examine long-term precipitation change effects on streamflow in this study. To evaluate the impacts of temperature and wind speed change on streamflow regime, we compared historical climate with the trend removed climate conditions. Figure S4 shows  LUCC refers to cropland expansion, grassland change, urbanization, reservoir regulation, and so on. This study adopted the scenario method, which has been widely utilized [6,10,17,61] to investigate CC and LUCC impacts and to segregate their impacts. Six scenarios were used: (1) S1, a baseline represented by climate in 1966 and land cover in 1960; (2) S2, represented by historical wind speed condition, temperature in 1966, and land cover in 1960. The comparison between S2 and S1 showed the impact of wind speed change; (3) S3 was represented by historical temperature condition, wind speed in 1966, and land cover in 1960. The comparison between S3 and S1 denoted the impact of temperature change; (4) S4 was represented by historical climate conditions and land cover in 1960. The comparison between S4 and S1 showed the influence of CC in general; (5) S5 was represented by a series of simulations in the sequence using historical climate conditions bracketing land cover conditions in corresponding periods, that is, simulations with land cover in 1960 with 1966-1970 climate, 1980 land cover with 1971-1990 climate, 2000 land cover with 1991-2010 climate, and 2015 land cover with that of 2011-2016. The comparison between S5 and S4 revealed the impacts of LUCC without reservoir impacts; (6) S5 R was the observed streamflow at QT station. The comparison of S5 R and S4 represented the LUCC, including reservoir operation influences on streamflow at QT with the rationale that observed streamflow was affected by the combined impacts of CC and LUCC.
The full streamflow regime represented by magnitude, frequency, variability, timing, and duration [35] was examined. This study selected 18 streamflow metrics on the basis of the Indicators of Hydrologic Alteration (IHA) [62,63]. A detailed description of these metrics is provided in Table 6. The rigorously calibrated and validated DHSVM model was then employed to compare the differences of flow regimes among all the scenarios.   (Figure 5c). QT had decreased flow in the warm season, except June, and increased flow in autumn and winter. The increased winter flow at the four stations could be due to cryosphere degradation. The combined impact of wind speed and temperature generally resembled that of temperature impact, except in June and July at QT, indicating temperature was the predominant impact factor for changing monthly streamflow magnitude (Figure 5d). The impacts of LUCC on monthly streamflow amplified from the upstream to the downstream areas, largely reflecting the reality in the basin (Figure 5e). As cropland expansion was the main feature of land cover change taking place in the middle to lower reaches in 1980, NC and HL streamflow was barely affected, whereas streamflow at XM dropped slightly during May-September (reduced by 0.23 m 3 s −1 ). Streamflow at QT trended similarly with XM, but with 10-fold reduction from May to October without considering the reservoir (left y-axis). However, with the reservoir, the change was mixed, with increases and decreases across months (Figure 5e, right y-axis). Streamflow reduced by 10.65 m 3 s −1 on average in June-September, and increased by 4.37 m 3 s −1 on average in March and April, resulting from flood control in the rainy season and water release for spring irrigation. Under the combined CC and LUCC (Figure 5f), monthly streamflow at NC, HL, and XM increased in most months, similar to the changes caused by CC impacts, indicating the dominance of CC impacts. Although monthly streamflow without reservoir influence at QT had similar variability caused by CC impacts (Figure 5d), the magnitude of monthly streamflow was mainly negative, like that by LUCC impacts (Figure 5e), with the greatest decline (4.19 m 3 s −1 ) in August. With reservoir regulation, October-May streamflow increased, whereas June-September streamflow decreased dramatically, especially in August, which dropped markedly by 14.93 m 3 s −1 . Figure 6 demonstrates that extreme streamflow (Q max1 and Q max3 , hereafter Q max with reference to left y-axis; Q min1 and Q min3 , hereafter Q min with reference to right y-axis) increased with basin area (Figure 6a) and responded diversely to CC and LUCC. Wind speed decrease resulted in higher Q max , which was enhanced in the downstream area, but had a smaller negative change in Q min (Figure 6b). With greater temperature, Q max increased at NC, HL, and XM, but decreased at QT, whereas Q min consistently increased in the whole basin (Figure 6c). Figure 6b-d also shows that temperature was the dominant impact factor on Q max and Q min in the whole river basin. The impacts of LUCC with and without the reservoir were different (Figure 6e). LUCC without the reservoir displayed negative effects on all Q max and Q min , whereas LUCC with the reservoir had positive effect on Q min but greater negative influence on Q max . By comparing Figure 6d-f, it can be found that CC impacts on Q max and Q min were dominant in the upstream (NC and HL) and the midstream (XM) areas. As for the downstream area (QT), LUCC was the principal factor for extreme streamflow change.
The variability of extreme streamflow (CV max1 and CV max3 , hereafter CV max ; CV min1 and CV min3 , hereafter CV min ) are depicted in Figure 7. Average values of CV max and CV min in the BRB in baseline scenario were around 55% and 12%, respectively. In general, the two indices showed an inverse relationship with basin area, with CV max decreasing and CV min increasing with area ( Figure 7a). CV max decreased with the wind speed change at all stations, whereas CV min decreased at NC and HL, but increased at XM and QT (Figure 7b). Temperature increasing led to declining CV max but rising CV min at each station (Figure 7c). Comparing Figure 7b-d shows that CV max and CV min were dominated by temperature increase, except CV max at QT. LUCC, especially LUCC with the reservoir, exerted positive effects on CV max1 and CV min (Figure 7e). The great positive effects on variability could make extreme streamflow less predictable. The impacts of the combined CC and LUCC (Figure 7f), especially including the reservoir, increased the variability of CV min at all stations. However, CC and LUCC effects had negative impacts on CV max at most stations. Figures 6 and 7 show inverse change in Q max and CV max , but similar change in Q min and CV min caused by CC and LUCC. This is because streamflow magnitude changes can arouse standard deviation change of data series, and different magnitudes of extreme flow can result in various trends of variability.

Figure 5.
Monthly streamflow of baseline scenario S1 (a) and its change in scenarios of wind speed change (b); temperature change (c); wind speed and temperature change, that is, climate change (d); land use/cover change (e); and both climate change and land use/cover change (f). QTR represents Qiaotou streamflow change caused by land use/cover change, which covered reservoir regulation (the same below for QTR). Figure 6 demonstrates that extreme streamflow (Qmax1 and Qmax3, hereafter Qmax with reference to left y-axis; Qmin1 and Qmin3, hereafter Qmin with reference to right y-axis) increased with basin area (Figure 6a) and responded diversely to CC and LUCC. Wind speed decrease resulted in higher Qmax, which was enhanced in the downstream area, but had a smaller negative change in Qmin (Figure 6b). With greater temperature, Qmax increased at NC, HL, and XM, but decreased at QT, whereas Qmin consistently increased in the whole basin (Figure 6c). Figure 6b-d also shows that temperature was the dominant impact factor on Qmax and Qmin in the whole river basin. The impacts of LUCC with and without the reservoir were different (Figure 6e). LUCC without the reservoir displayed negative effects on all Qmax and Qmin, whereas LUCC with the reservoir had positive effect on Qmin but greater negative influence on Qmax. By comparing Figure 6d-f, it can be found that CC impacts on Qmax and Qmin were dominant in the upstream (NC and HL) and the midstream (XM) areas. As for the downstream area (QT), LUCC was the principal factor for extreme streamflow change. Monthly streamflow of baseline scenario S1 (a) and its change in scenarios of wind speed change (b); temperature change (c); wind speed and temperature change, that is, climate change (d); land use/cover change (e); and both climate change and land use/cover change (f). QT R represents Qiaotou streamflow change caused by land use/cover change, which covered reservoir regulation (the same below for QT R ). The variability of extreme streamflow (CVmax1 and CVmax3, hereafter CVmax; CVmin1 and CVmin3, hereafter CVmin) are depicted in Figure 7. Average values of CVmax and CVmin in the BRB in baseline scenario were around 55% and 12%, respectively. In general, the two indices showed an inverse relationship with basin area, with CVmax decreasing and CVmin increasing with area ( Figure 7a). CVmax decreased with the wind speed change at all stations, whereas CVmin decreased at NC and HL, but increased at XM and QT (Figure 7b). Temperature increasing led to declining CVmax but rising CVmin at each station (Figure 7c). Comparing Figure 7b-d shows that CVmax and CVmin were dominated by temperature increase, except CVmax at QT. LUCC, especially LUCC with the reservoir, exerted positive effects on CVmax1 and CVmin (Figure 7e). The great positive effects on variability could make extreme streamflow less predictable. The impacts of the combined CC and LUCC (Figure 7f), especially including the reservoir, increased the variability of CVmin at all stations. However, CC and LUCC effects had negative impacts on CVmax at most stations. Figures 6 and 7 show inverse change in Qmax and CVmax, but similar change in Qmin and CVmin caused by CC and LUCC. This is because streamflow magnitude changes can arouse standard deviation change of data series, and different magnitudes of extreme flow can result in various trends of variability.

Frequency
Non-exceedance probability of annual maximum peakflows for each station are shown in Figure  8. Peakflows with 95% non-exceedance probability for baseline (S1) corresponded to 78. 35, 29.47, 124.26, and 323.34 m 3 s −1 at NC, HL, XM, and QT, respectively. CC and LUCC caused little change in NC, HL, and XM peakflow frequency, with characteristics the same as CC and LUCC impacts on their magnitude of extreme flow. Peakflow frequency changed moderately at QT affected by CC and LUCC. The 95% non-exceedance probability of peakflow for S1, S4, and S5R corresponded to 323.34, 285.11, and 190.2 m 3 s −1 , respectively.

Frequency
Non-exceedance probability of annual maximum peakflows for each station are shown in Figure 8. Peakflows with 95% non-exceedance probability for baseline (S1) corresponded to 78. 35, 29.47, 124.26, and 323.34 m 3 s −1 at NC, HL, XM, and QT, respectively. CC and LUCC caused little change in NC, HL, and XM peakflow frequency, with characteristics the same as CC and LUCC impacts on their magnitude of extreme flow. Peakflow frequency changed moderately at QT affected by CC and LUCC. The 95% non-exceedance probability of peakflow for S1, S4, and S5 R corresponded to 323.34, 285.11, and 190.2 m 3 s −1 , respectively.

Timing
Three timing indices were selected-the occurrences of centroid of annual mean daily streamflow (JDCT) [63], annual maximum peakflow (JDmax), and annual minimum streamflow (JDmin), all represented by Julian day. In baseline scenario, JDCT, JDmax, and JDmin appeared on the 214th, 218th, and 88th days on average for all stations, respectively, during the course of a year on average, which were in early August for JDmax and JDCT and in late March for JDmin (Figure 9a). Figure 9b shows that wind speed had little influence on JDCT, JDmax, and JDmin, with only a 1 day difference for most stations. Temperature rising, on the one hand, postponed JDCT, but advanced JDmin and resulted in mixed change in JDmax (Figure 9c). Figure 9d reveals that JDmin appeared earlier due to CC, similar to

Timing
Three timing indices were selected-the occurrences of centroid of annual mean daily streamflow (JD CT ) [63], annual maximum peakflow (JD max ), and annual minimum streamflow (JD min ), all represented by Julian day. In baseline scenario, JD CT , JD max , and JD min appeared on the 214th, 218th, and 88th days on average for all stations, respectively, during the course of a year on average, which were in early August for JD max and JD CT and in late March for JD min (Figure 9a). Figure 9b shows that wind speed had little influence on JD CT , JD max , and JD min , with only a 1 day difference for most stations. Temperature rising, on the one hand, postponed JD CT , but advanced JD min and resulted in mixed change in JD max (Figure 9c). Figure 9d reveals that JD min appeared earlier due to CC, similar to temperature increase effect, whereas JD CT and JD max had various changes at each station. Figure 9e shows that when there was no reservoir, LUCC changed a little in JD CT . Moreover, there were 8 days in advance for JD max and JD CT after the reservoir operation (solid triangle of QT R in Figure 9e, values on the right y-axis). Additionally, under the reservoir operation, JD min was delayed for 62 days compared with baseline, which was in late May. On the basis of Figure 9d-f, CC played a leading role in changing timing indices for NC, HL, and XM. However, CC and LUCC showed equal impacts on QT timing indices when there was no reservoir (left y-axis), and LUCC was the dominant factor when reservoir began to operate (right y-axis). Overall, timing indices only changed greatly under reservoir operation, whereas in other situations the changes in timing was within 7 days.
Water 2020, 12, x FOR PEER REVIEW 18 of 26 temperature increase effect, whereas JDCT and JDmax had various changes at each station. Figure 9e shows that when there was no reservoir, LUCC changed a little in JDCT. Moreover, there were 8 days in advance for JDmax and JDCT after the reservoir operation (solid triangle of QTR in Figure 9e, values on the right y-axis). Additionally, under the reservoir operation, JDmin was delayed for 62 days compared with baseline, which was in late May. On the basis of Figure 9d-f, CC played a leading role in changing timing indices for NC, HL, and XM. However, CC and LUCC showed equal impacts on QT timing indices when there was no reservoir (left y-axis), and LUCC was the dominant factor when reservoir began to operate (right y-axis). Overall, timing indices only changed greatly under reservoir operation, whereas in other situations the changes in timing was within 7 days.

Duration
Three streamflow duration indices encompassing number of days with daily streamflow higher than annual mean streamflow (DAve), higher than high streamflow (DHigh), or lower than low streamflow (DLow) were chosen in this study. At baseline, the average durations were 128-132 days, 51-57 days, and 2-7days for DAve, DHigh, and DLow, respectively (Figure 10a). Compared to baseline, wind speed change only lengthened DAve within 2 days, but exerted no influence on DHigh and DLow (Figure 10b). Temperature change lengthened the duration of DHigh by 1-3 days, but cut down DLow

Duration
Three streamflow duration indices encompassing number of days with daily streamflow higher than annual mean streamflow (D Ave ), higher than high streamflow (D High ), or lower than low streamflow (D Low ) were chosen in this study. At baseline, the average durations were 128-132 days, 51-57 days, and 2-7days for D Ave , D High , and D Low , respectively (Figure 10a). Compared to baseline, wind speed change only lengthened D Ave within 2 days, but exerted no influence on D High and D Low (Figure 10b). Temperature change lengthened the duration of D High by 1-3 days, but cut down D Low by about 0-4 days, and also exerted mixed influences on D Ave at different stations (Figure 10c). Combined CC effect was dominated by temperature change effect (Figure 10d). As for LUCC, it had almost no influence when there was no reservoir (Figure 10e). However, once the reservoir was put into operation, D Ave , D High , and D Low were all extended, especially D Low , which was prolonged by 15 days. Figure 10f suggests that the combined impacts of CC and LUCC were generally similar to CC impacts when there was no reservoir, and reservoir operation on top of CC and LUCC still had the greatest impacts on the duration. by about 0-4 days, and also exerted mixed influences on DAve at different stations (Figure 10c). Combined CC effect was dominated by temperature change effect ( Figure 10d). As for LUCC, it had almost no influence when there was no reservoir (Figure 10e). However, once the reservoir was put into operation, DAve, DHigh, and DLow were all extended, especially DLow, which was prolonged by 15 days. Figure 10f suggests that the combined impacts of CC and LUCC were generally similar to CC impacts when there was no reservoir, and reservoir operation on top of CC and LUCC still had the greatest impacts on the duration. Figure 10. Duration of streamflow higher than annual mean streamflow (DAve), higher than high streamflow (DHigh), or lower than low streamflow (DLow) in baseline scenario S1 (a), and the impacts of wind speed (b), temperature (c), climate change (d), land use/cover change (e), and both climate change and land use/cover change (f).

Deconvolve the Impacts of CC and LUCC
The coefficients of variation in annual streamflow (CVyr) for all scenarios fluctuated slightly around 24% (Table 7). Temperature played a more important role in changing CVyr relative to wind speed. Comparing the impacts of CC and LUCC on CVyr, it is clear that CC impact was dominant at NC and HL in the upstream area, whereas LUCC with or without reservoir was dominant at XM and Figure 10. Duration of streamflow higher than annual mean streamflow (D Ave ), higher than high streamflow (D High ), or lower than low streamflow (D Low ) in baseline scenario S1 (a), and the impacts of wind speed (b), temperature (c), climate change (d), land use/cover change (e), and both climate change and land use/cover change (f).

Deconvolve the Impacts of CC and LUCC
The coefficients of variation in annual streamflow (CV yr ) for all scenarios fluctuated slightly around 24% (Table 7). Temperature played a more important role in changing CV yr relative to wind speed. Comparing the impacts of CC and LUCC on CV yr , it is clear that CC impact was dominant at NC and HL in the upstream area, whereas LUCC with or without reservoir was dominant at XM and QT in the downstream area. In general, combined CC and LUCC moderately increased CV yr at NC, XM, and QT, but decreased CV yr at HL.  Table 8 demonstrates the deconvolved effects of multiple impact factors. During the period of 1966-2016, most changes caused by CC and LUCC on annual mean streamflow passed the 0.01 significance test. On the basis of the first two columns of WS and T, it is clear that temperature affected mean annual streamflow more than wind speed, except for QT. In addition, the third column of CC was not the linear addition of temperature and wind speed impacts. Streamflow at different stations responded differently to CC and LUCC. NC, HL, and XM stations, which are in the upstream and midstream areas, were mostly affected by CC, which accounted for 100%, 81.25%, and 88.89%, respectively (Table 8), indicating that streamflow change at NC was totally controlled by CC because of the protection from national natural reserve. Streamflow at QT, the outlet of the BRB, was decreased by 1.96 m 3 s −1 (p < 0.01) due to CC and LUCC, accounting for 10.06% of the mean annual streamflow of S1 at QT. LUCC was the predominant factor, causing 84.51% reduction in mean annual streamflow, whereas CC accounted for only 15.49%. In general, upstream and midstream flow were more controlled by CC, whereas downstream flow was mainly predominated by LUCC.

Discussion
Generally, the DHSVM provided good agreement between simulation and observation in the overlapping periods (Table 5 and Figure S3). However, some discrepancies appeared in monthly simulated and observed streamflow at NC and HL ( Figure S3). For example, in the line graph, maximum monthly streamflow was in September for the observation but in July for the simulation, and spring streamflow was underestimated. Both of them were caused by inaccurate estimation of cryospheric variables such as snow water equivalent (SWE) and frozen soil, which will be investigated in the future.
Depending on the elevation and location of the intensity of human activity, streamflow regime changes were different in the BRB. Decreased wind speed restrained the annual evapotranspiration by −0.24 mm in the BRB on average, resulting in increases in the means of maximum, monthly, and annual streamflow (Figures 5b and 6b; Table 8), leading to decreases in streamflow variability ( Figure 7b and Table 7). Although annual wind speed decreased significantly, monthly wind speed trends in July-September increased slightly (Figure 2b), which resulted in decreased streamflow in July-September, shuffling the timing of the maximum daily flow by usually advancing its occurrence; thus, JD max appeared earlier (Figure 9b). In the meantime, D Ave lengthened consistently in this river basin (Figure 10b), which happened in April-June when great increase in streamflow occurred.
Temperature increase impacts on evaporation could be explained by the following three ways. According to the Clausius-Clapeyron equation, saturation water vapor increases with air temperature, and thus vapor pressure deficit between land surface and air increases. The rate of evaporation is dependent on the amount of heat transferred, and warmer dry air may supply more energy to an evaporating surface. Less energy is required to evaporate warm water than cool water, as latent heat of vaporization (Lv) becomes smaller with higher temperatures, for instance, Lv at 20 • C and 35 • C are 2.44 and 2.40 MJ kg −1 , respectively. For the same energy input, more water will be evaporated at the higher temperature. In the BRB, temperature increase intensified the mean annual evapotranspiration by 51.57 mm on average in the BRB, and a similar amount in the upper basins.
Temperature increase also reduced mean daily SWE by 20.86 mm on average in the entire BRB, and by 32.01 mm and 26.47 mm in the upstream basins above NC and HL, respectively, on the basis of the model scheme. Reduction in SWE and increase in evapotranspiration due to temperature increase interrupted the original partition among water balance terms over the simulation period across the BRB, hence changing the streamflow regime indices as shown above. In the upstream (NC and HL) and midstream (XM), snow melt increase was greater than evapotranspiration increase over the year on average, which resulted in greater streamflow magnitude at most monthly and annual time steps (Figures 5c and 6c; Table 8). In the lower elevation, SWE increase could not balance the evaporation increase, leading to decreased maximum and annual flow in lower stream (Figure 6c and Table 8). With temperature increase, starting date of snow melt was brought forward, making lower JD min , whereas JD CT was postponed due to the delaying of starting freeze date, resulting in greater streamflow in September and October (Figure 9c).
For CC impacts, streamflow regime variables, such as magnitude (annual, monthly, and extreme streamflow), variability (annual, extreme flow), timing (minimum flow), and duration were more influenced by temperature than by wind speed at most stations. Notwithstanding, effects of wind speed change on mean annual and monthly streamflow, Q max , CV max , and CV yr increased from upstream to downstream, which indicated significant spatially accumulative effects that wind speed change exerted on streamflow. By contrast, temperature change had completely opposite impacts on some streamflow regime indices, especially for annual streamflow and CV max , which indicated that temperature is more important at high elevation. The accumulative effects resulted in greater wind speed influence than temperature effect on mean annual streamflow and CV max at the downstream QT (Table 8 and Figure 7b). It is thus inferred that accumulative effects of wind speed change on streamflow may be more severe in large river basins, requiring more studies to evaluate this phenomenon, as such studies are lacking. The analysis also showed that combined climate change impacts were not the linear addition of temperature and wind speed impacts (Tables 7 and 8).
Since the 1960s and 2001, the downstream BRB has experienced cropland expansion and Heiquan reservoir operation, respectively. The former increased evapotranspiration by 15.29 mm in the BRB, and thus led to the reduction of the flow magnitude (Figures 5e and 6e; Table 8) and peakflow (Figure 8), which resulted in lower means in maximum and annual flow and increased their variability (Figure 7e and Table 7). Reservoir operation caused re-distribution of monthly streamflow, decreasing maximum streamflow and its variability during the flooding season (Figures 5e and 6e), but increasing them for minimum streamflow (Figures 6e and 7e). This change narrowed the gap between maximum and minimum flow throughout the year, making intra-annual streamflow pattern flatter and prolonging durations (Figure 10e). The operation could also make JD min occur on any day in a year, not only in March, which could postpone JD min considerably. On the other hand, with the noteworthy increase of streamflow in January-April after reservoir operation, JD CT appeared earlier (Figure 9e). At QT mean, annual streamflow declined by 6.78% due to LUCC without considering reservoir operation. The Heiquan reservoir also caused it to fall by 5.54% (12.32% minus 6.78%; Table 8). It seems that the effects of cropland expansion on mean annual streamflow were slightly larger than that of reservoir regulation. However, Figure 5f shows that reservoir regulation played a stronger role in intra-annual streamflow re-distribution.
As there was no streamflow observation right below the reservoir after 2001 when the reservoir started operation and reservoir information was recorded, the reservoir effects cannot be accounted for in this study. When examining LUCC (including the reservoir) impacts on river flow, we used a simplified approach, as mentioned above. Reservoir regulation rules, reservoir storage, and open water evaporation were not considered, as they would introduce errors and bring uncertainty into the results. In a future study, we plan to incorporate a reservoir module by considering reservoir storage and regulations with the DHSVM to further investigate reservoir impacts on streamflow in detail.
Model uncertainties result from input data, model structure, and parameter values. In this study, we only examined parameter-related uncertainty. As we have relatively high confidence in our model input data, including soil and vegetation types and climate forcing, the uncertainty from input data was not analyzed. As for model structure uncertainty, it was beyond the scope of the current study, but certainly can be examined in the future by using an ensemble of hydrological models.

Conclusions
The BRB experienced prominent upward trends in temperature and downward trends in wind speed at p < 0.01 from 1966 to 2016, with abrupt changes appearing in the 1980s, whereas annual precipitation and streamflow exhibited non-significant trends.
In general, the impacts of temperature change dominated over those of wind speed change, especially at high elevation, but accumulative effects of wind speed change on streamflow should not be ignored. As for magnitude, CC enhanced annual and maximum flow in the upstream and midstream areas but reduced those in the downstream area, and also increased minimum flow at all stations. For variability, annual and minimum flow showed upward trends, wheras maximum flow decreased. The frequency of peakflow under CC was changed little. For timing, CC brought the early occurrence of daily minimum streamflow. The duration of maximum streamflow was always extended.
Since the 1960s, cropland expansion, interconversion between grassland and other land cover types, and the Heiquan reservoir operation were the major characteristics of BRB LUCC. Reservoir regulation was the predominant factor altering the streamflow regime in the downstream area. It significantly increased the streamflow variability and extended the duration of all indices while decreasing the maximum flow and flood frequency. Additionally, there was earlier timing for maximum and temporal centroid streamflow, but later timing for minimum streamflow. Cropland expansion and reservoir regulation were the two major land cover type changes, with the former exerting greater impacts on annual streamflow, and the later primarily affecting intra-annual streamflow distribution.
Spatially, upstream and midstream flow in high elevation were more sensitive to CC, whereas downstream flow in low elevation was mainly predominated by LUCC.
The results of this study expand our understanding of the wind speed change impacts on hydrological processes, and also provide evidence to support mitigation schemes to better manage the water resources synergistically in this river basin.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/12/4/1198/s1, Figure S1: Conceptual representation (a) and flow chart (b) of the DHSVM. Figure S2: Change point detection of annual data. Mann-Kendall test on annual temperature at DT (T, (a)), moving t-test technique applied to annual temperature at DT (T, (b)), Mann-Kendall test on annual wind speed at DT (WS, (c)). Figure Table S1: Transition matrix of land use in the Beichuan River Basin from 1980 to 2000 (km 2 ). Table S2