Climate Change Impact on Spatiotemporal Hotspots of Hydrologic Ecosystem Services : A Case Study of Chinan Catchment , Taiwan

Hydrologic ecosystem services are greatly affected by the changing climate. In this study, the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model was used to quantify hydrologic ecosystem services. Five general circulation models (GCMs) and two representative concentration pathways (RCPs) were selected to estimate hydrologic ecosystem services. The Local Indicators of Spatial Association (LISA) index was used to identify hydrologic ecosystem hotspots. The hotspots were used to evaluate the impact of climate change on the services. Results indicate that annual water yields vary from −17% to 8%, with significant intra-year fluctuation. Compared to baseline data, the CESM1-CAM5 predicts an increase of 45% in June, but HadGEM2-AO predicts a drop to only 12% in January. Sediment export results show a similar trend to water yield, with sediment export increasing significantly under RCP 8.5, and monthly sediment export increases concentrated from June and October. Nitrogen and phosphorous exports both show less significant changes but obvious intra-year variations. The CESM1-CAM5 predicts strong seasonal and spatial variation of the hydrologic ecosystem services. Our proposed approach successfully identifies annual and monthly hotspot spatial changes of hydrologic ecosystem services under climate change.


Introduction
Hydrological variation and water availability risks increase with climate changes, and yet water demand is increasing in most countries around the world [1].To deal with such risks, many countries have developed frameworks to achieve convergence between water supply and demand, and to ensure freshwater ecosystem services are sustained [1].Ecosystem services are directly supplied by ecosystem major components such as soil, water, vegetation and biota, or indirectly through their interactions [2].These linkages support practical implementations of water ecosystem services assessment and valuation methodologies [3].However, these details as well as additional ecosystem services have not yet been taken into serious consideration during river basin planning, and hydrologic ecosystem services management decision-making processes.Moreover, hydrologic ecosystem services refer to any and all of the water-related services provided to humans who rely on the ecosystem [4][5][6][7][8][9][10].The regional effects of hydrologic ecosystem services then, include benefitting both people downstream, and ecosystems throughout the watershed.Hydrologic ecosystem services can be divided into five categories: Diverted water supply, in situ water supply, water damage mitigation, spiritual and aesthetic, and supporting services [4,8].Improvement of water quality and quantity is just one type of hydrologic ecosystem service [11,12].Therefore, at the watershed scale, hydrological cycles connect to major ecosystem services such as water purification, water retention, and climate regulation [3,8].Water purification, represented by water quality, influences human well-being in various ways.For example, the effects of increased nutrient loads in water bodies, and the unevenly and spatiotemporally distributed water quality benefits different groups have access to [5].In terms of water quantity, since this hydrologic ecosystem service both directly and indirectly affects the supply and demand of water resources, it is crucial for human survival and development.Many studies have already demonstrated positive results when integrating ecosystem services into watershed management (e.g., the EU Water Framework Directive), and have recommended its inclusion as an essential decision indicator [3,6,13,14].
While the Millennium Ecosystem Assessment Board proposed an assessment of the effects of ecosystem change on human well-being [12], hydrologic ecosystem services should be defined and assessed for inclusion in decision-making processes [9,13].Further, assessments should be conducted primarily through mature-market or nonmarket activities in order to quantify the effect of hydrologic ecosystem services on human welfare [15].Hydrologic ecosystem services are sustained by maintaining particular ecosystem services or economic benefits [16,17].If a hydrologic ecosystem service is irreplaceable or is at risk of climate change-induced changes to the supply of hydrologic ecosystem services, then assessing hydrologic ecosystem services should be required and necessary for water-resource governance and decision-making processes.For example, hydrologic ecosystem services have been included in the decision-making process in the urban water-supply system [18]; water, food, energy supply [19]; and water-resource management [3].Recently, the most common concepts used in water governance are integrated water resources management (IWRM) and ecosystem services.Ecosystem services must constantly be reexamined to ensure their sustainability [20].Using ecosystem services as a foundation for developing an adaptive IWRM program [3,6,21,22] therefore enables evaluations of the balance between humans and the environment from a macroscopic perspective; and enables developing a multiple-decision model for water-resource management based on the concept of sustainability.This is an effective way to delineate spatiotemporal distributions of hydrologic ecosystem services for decision-making and water resources management.
Climate change influences the hydrological cycles and land cover that in turn influences surface albedo with feedback to the climate system [23][24][25][26][27].As the global environment continues to evolve, climate change is a reality that we all experience [28].In 2012, the Intergovernmental Panel on Climate Change (IPCC) published the Fifth Assessment Report (AR5) and confirmed that global warming will persist.Global climate change has threatened the stability of ecosystems and the value of the services they provide [29].Of the 24 ecosystem services, 60% have been confirmed to be in a state of continual degradation.This current state will continue to degrade unless appropriate measures are implemented [12].Studies have indicated that both climate change and human activities affect the value of ecosystem services [30].Further, climate change influences hydrological conditions and runoff [31].Changes in temperature and rainfall patterns reduce sediment retention and runoff capacity, causing negative effects on hydrologic ecosystem services, such as water availability and quality [32,33].The sustainable management of water resources becomes a more complex matter when combined with the uncertainty of climate change [34].It is therefore essential to identify spatiotemporal hotspots of hydrologic ecosystem services under climate change prior to watershed management decision-making.
Water-resource management in Taiwan is difficult because of the island's unique landform and rainfall patterns.Water resources are substantially affected by deteriorating climatic conditions.Older water-resource adaptive programs in which hydrologic ecosystem services had not been designed and fully considered may be insufficient under drastic global climate changes [35], particularly in areas such as Hualien, Taiwan [36] (along Taiwan's east coast), since it has more access to water than other areas in Taiwan.Hualien is the largest county by area in Taiwan, though its landform primarily consists of mountainous terrain.Since no artificial lakes or reservoirs have been constructed there, rainfall is the primary water source.Studies have reported the risks for a Hualien-area water shortage caused by climate change.Peng et al. [37], for example, simulated both the flow rate of each watershed in Hualien using generalized watershed loading functions [38], as well as the shortage index of each water-supply system using a system dynamics model [39].The purpose of this study, however, is to This study selected the Hualien-area water-supply system with the highest water deficiency, Chinan catchment, as the study area to explore changes in four hydrologic ecosystem services under a future climate scenario.The effect of climate change on the functions of hydrologic ecosystem services in the study area was assessed.Changes in hydrologic ecosystem services were demonstrated using spatial characteristics, which is a useful method for enabling decision makers to understand the situation and make informed decisions.To this end, we performed a regional spatial-correlation calculation to demonstrate hotspot spatial distributions of monthly hydrologic ecosystem service functions in order to inform watershed management decision-making.

Materials and Methods
To assess the effects of climate change on hydrologic ecosystem services in the Chinan catchment area, the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) was used to simulate the changes of the four hydrologic ecosystem service functions: Water yield, sediment export, nitrogen (N) export, and phosphorus (P) export.To identify key locations for water-resource management, the Local Indicators of Spatial Association (LISA) model was used to calculate clusters of the four hydrologic ecosystem services within the spatial distribution of the study area.These identified hotspots then, are where immediate attention is required.

Study Area
The Chinan water-supply system is a part of the Lau River Basin.The Lau River flows through Shoufeng and Sioulin Townships and is a tributary of the Hualien River.Located in eastern Taiwan, the Hualien River is a river managed by the central government and is one of the primary rivers in Hualien County.Its source is the Guangfu River, which originates from the Bazi Mountain.The main stream of the Hualien River is 57.28 km long, and the drainage system is 81 km long.The drainage area is 1507.09km 2 , and the mean discharge is 102 m 3 /s.The drainage system flows northeast along the Longitudinal Valley and enters the Pacific Ocean at Lingding at the northern tip of the Coastal Range.
Land use for the Chinan catchment area is displayed in Figure 1.The majority of the land is covered by forest.The Lau River flows eastwards through the catchment area, while the downstream region is for agricultural use (orchards).Some mining lands are located to the west of the catchment area, and some road-covered lands are scattered around its edges.shortage index of each water-supply system using a system dynamics model [39].The purpose of this study, however, is to explore the effects of poor weather conditions on the Hualien-area hydrologic ecosystem services in the context of climate change.This study selected the Hualien-area water-supply system with the highest water deficiency, Chinan catchment, as the study area to explore changes in four hydrologic ecosystem services under a future climate scenario.The effect of climate change on the functions of hydrologic ecosystem services in the study area was assessed.Changes in hydrologic ecosystem services were demonstrated using spatial characteristics, which is a useful method for enabling decision makers to understand the situation and make informed decisions.To this end, we performed a regional spatialcorrelation calculation to demonstrate hotspot spatial distributions of monthly hydrologic ecosystem service functions in order to inform watershed management decision-making.

Materials and Methods
To assess the effects of climate change on hydrologic ecosystem services in the Chinan catchment area, the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) was used to simulate the changes of the four hydrologic ecosystem service functions: Water yield, sediment export, nitrogen (N) export, and phosphorus (P) export.To identify key locations for water-resource management, the Local Indicators of Spatial Association (LISA) model was used to calculate clusters of the four hydrologic ecosystem services within the spatial distribution of the study area.These identified hotspots then, are where immediate attention is required.

Study Area
The Chinan water-supply system is a part of the Lau River Basin.The Lau River flows through Shoufeng and Sioulin Townships and is a tributary of the Hualien River.Located in eastern Taiwan, the Hualien River is a river managed by the central government and is one of the primary rivers in Hualien County.Its source is the Guangfu River, which originates from the Bazi Mountain.The main stream of the Hualien River is 57.28 km long, and the drainage system is 81 km long.The drainage area is 1507.09km 2 , and the mean discharge is 102 m 3 /s.The drainage system flows northeast along the Longitudinal Valley and enters the Pacific Ocean at Lingding at the northern tip of the Coastal Range.The primary use of water in the Hualien area is for agricultural purposes, accounting for 92% of water use, whereas domestic water accounts for 4.2%, and industrial water only accounts for 3.8%.The primary use of water in the Hualien area is for agricultural purposes, accounting for 92% of water use, whereas domestic water accounts for 4.2%, and industrial water only accounts for 3.8%.In total, 96.3% of agricultural water is used for irrigation.Regarding water supply, the principal source is surface water, accounting for 92% of the total water supply [41].If the surface water supply is insufficient, groundwater is allocated for supplementation.The water-supply structure in the Hualien area is depicted in Figure 2. Chinan is one of the 14 water-supply systems, supplying water for the northern jurisdiction.In total, 96.3% of agricultural water is used for irrigation.Regarding water supply, the principal source is surface water, accounting for 92% of the total water supply [41].If the surface water supply is insufficient, groundwater is allocated for supplementation.The water-supply structure in the Hualien area is depicted in Figure 2. Chinan is one of the 14 water-supply systems, supplying water for the northern jurisdiction.

Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST)
InVEST is an application software combined with a geographic information system designed by The Natural Capital Project [42].The model is used to investigate changes in human interest caused by ecosystem changes, employing conditions and processes that produce ecosystem service functions as production functions to quantify and evaluate the effects of changes in the level of ecosystem service output.The InVEST model can assess and quantify 14 final ecosystem services (directly contributing to human welfare), including carbon storage and sequestration, water yield, nutrient retention, and sediment retention [42].InVEST is a model tool for estimating and explicitly demonstrating the spatial information of ecosystem services [42,43].This study used the InVEST model (InVEST3.3.3) to quantify four hydrologic ecosystem services, including water yield, sediment export, N export, and P export.Rainfall, evapotranspiration, and infiltration were used as the InVEST model input data to calculate water yield.Sediment export for each analyzed cell refers to the total amount of sediment exported from each cell that reaches the waterbody.The primary factors affecting the amount of sediment export are rainfall amount, rainfall intensity, soil property, topography, vegetation, and land use.Increase in sediment export affects water quality and reservoir management [42].The sediment export is calculated from the Universal Soil Loss Equation (USLE) and is used to calculate the amount of soil deposited downstream along the streamline of a river [42].The yield of the nutrients, including N and P, are calculated according to the loading of nutrients of each land use type.The export of nutrients is then calculated by the water yield, nutrient yields, and digital elevation model (DEM) data.The InVEST user guide can be consulted for further details [42].

Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST)
InVEST is an application software combined with a geographic information system designed by The Natural Capital Project [42].The model is used to investigate changes in human interest caused by ecosystem changes, employing conditions and processes that produce ecosystem service functions as production functions to quantify and evaluate the effects of changes in the level of ecosystem service output.The InVEST model can assess and quantify 14 final ecosystem services (directly contributing to human welfare), including carbon storage and sequestration, water yield, nutrient retention, and sediment retention [42].InVEST is a model tool for estimating and explicitly demonstrating the spatial information of ecosystem services [42,43].This study used the InVEST model (InVEST3.3.3) to quantify four hydrologic ecosystem services, including water yield, sediment export, N export, and P export.Rainfall, evapotranspiration, and infiltration were used as the InVEST model input data to calculate water yield.Sediment export for each analyzed cell refers to the total amount of sediment exported from each cell that reaches the waterbody.The primary factors affecting the amount of sediment export are rainfall amount, rainfall intensity, soil property, topography, vegetation, and land use.Increase in sediment export affects water quality and reservoir management [42].The sediment export is calculated from the Universal Soil Loss Equation (USLE) and is used to calculate the amount of soil deposited downstream along the streamline of a river [42].The yield of the nutrients, including N and P, are calculated according to the loading of nutrients of each land use type.The export of nutrients is then calculated by the water yield, nutrient yields, and digital elevation model (DEM) data.The InVEST user guide can be consulted for further details [42].A DEM with a 40 × 40 m cell resolution from the Department of Land Administration, Ministry of the Interior was used.Data on land use was retrieved from the National Land Surveying and Mapping Center [44].Soil-depth data was obtained from the Soil and Water Conservation Bureau, Council of Agriculture, Executive Yuan.Plant available water capacity was calculated according to the InVEST user guide.The Chinan area soil erodibility (K) was informed by the Soil and Water Conservation Handbook [45], and the reference evapotranspiration rate was calculated using Hamon's Equation.The rainfall erosivity index (R) was calculated using the formula R = ap b , where p is the monthly rainfall; a and b are parameters which differ each month.The production rates of N and P exports surveyed by Wen and Chang [46] were used as InVEST model input data.

Local Indicators of Spatial Association (LISA)
LISA [47] was used to analyze the spatial correlation of each hydrologic ecosystem service, to identify the hotspots of hydrologic ecosystem service functions, and to present the distribution characteristics of hydrologic ecosystem service indicators.Hotspots were identified based on the values generated from the InVEST model.A weighted calculation method was then used on these hotspots to determine the optimal space for a key hydrologic ecosystem services protection area, for planning purposes.LISA provides a function for viewing outliers within the designated study area.LISA can display the degree to which the spatial unit affects the autocorrelation value of the entire research space with respect to two aspects: The aggregation point of the space in a region; and the degree of aggregation.In other words, the location of the cluster in the study area [47].For each spatial unit, LISA offers a statistical analysis of the value of spatial autocorrelation.That is, the similarity between each space and the adjacent region, and whether this value is greater than that of a farther region [48,49].These spatial autocorrelation values correspond to the local Moran I value areas (i.e., areas surrounded by features with similar values); and local Moran P values, which are substantially similar to its adjacent values [46].Spatial hotspots are thusly identified for each space in the study area.
By rejecting the null hypothesis that does not involve a spatial correlation, LISA determines the amount of spatial aggregation.If results pass a significance test and exhibit positive significance, then a positive spatial correlation is identified with attributes similar to that of the surrounding area.Observed values of the spatial unit and the surrounding area were either higher than average (HH), or lower than average (LL).A negative significance level indicated that the surrounding attributes were dissimilar, constituting a negative spatial correlation.Moreover, this negative significance suggests that the observed value of the spatial unit was higher than the average value, but the observed value of the surrounding area was lower than the average value (HL) or vice versa (LH).A non-significant result indicated no spatial correlation [50].Because the total area of the Chinan catchment is not large, the output of the hydrologic ecosystem service function exceeded the average grid.The area within 200 m of the grid must also exceed the average grid in order for the grid to be identified as a hotspot.

Climate Change Scenario and Models
The UN Intergovernmental Panel on Climate Change (IPCC) published the Fifth Assessment Report (AR5) in 2013 [51].In this report, future global climate projections have been simulated through a variety of scenarios, which are conceived of by increases of radiative forcing, that reflect trajectories of greenhouse gas increases in the atmosphere.Four representative concentration pathways (RCPs), namely RCP2.6,RCP4.5, RCP6.0, and RCP8.5, are labelled as such to represent corresponding radiative forcing values of +2.6, +4.5, +6.0, and +8.5 W/m 2 , respectively, in the year 2100, and are relative to pre-industrial revolution values.The projected future climate in 2100 is computed by atmosphere-ocean general circulation models (GCMs) [51]  World Meteorological Organization (WMO), the applicability for 41 GCMs have been evaluated.On a global scale, their results are generally in agreement.However, further evaluation is needed to test the effectiveness at the regional scale.
This study adopted the AR5 future climate scenarios simulated by the five general circulation models (GCMs) under two RCPs, 2.6 and 8.5 [51,52].Of these GCMs, three are suitable for eastern Taiwan:

InVEST Analysis Results for the Effect of Climate Change on Hydrologic Ecosystem Services
For average annual water yield, four hydrologic ecosystem service outputs increased the most in the RCP8.5 scenario when compared to the baseline.The average annual water yield increased by 8% (CESM1-CAM5), the sediment export increased by 11% (CESM1-CAM5), and the N and P exports increased by 1% (GISS-E2-R) and 3% (GISS-E2-R), respectively, compared with the baseline.Specifically, the water yield and sediment export in the RCP8.5 CESM1-CAM5 scenario increased the most, by 10% (Table 1).The annual trend for each scenario in RCP2.6 and RCP8.5 is the same.However, the water yield decreases during the dry season and increases during the flood season (Figure 3a,b).
The spatial distribution of water yield gradually decreases from west to east (Figure 4); January and September, for example.Since the primary parameter for calculating water yield is rainfall, the water yield and rainfall spatial distribution trends are the same, gradually decreasing from the higher western terrain to the lower eastern side.The spatial distribution of water yield gradually decreases from west to east (Figure 4); January and September, for example.Since the primary parameter for calculating water yield is rainfall, the water yield and rainfall spatial distribution trends are the same, gradually decreasing from the higher western terrain to the lower eastern side.Sediment export values in the Chinan catchment area under different climate change scenarios show that sediment export is substantially low from November to May, with the high sedimentexport season from June to October (Figure 5).The amount of sediment export is higher than the baseline, particularly in September.This trend is the same as that of the water yield.The spatial  The spatial distribution of water yield gradually decreases from west to east (Figure 4); January and September, for example.Since the primary parameter for calculating water yield is rainfall, the water yield and rainfall spatial distribution trends are the same, gradually decreasing from the higher western terrain to the lower eastern side.Sediment export values in the Chinan catchment area under different climate change scenarios show that sediment export is substantially low from November to May, with the high sedimentexport season from June to October (Figure 5).The amount of sediment export is higher than the baseline, particularly in September.This trend is the same as that of the water yield.The spatial Sediment export values in the Chinan catchment area under different climate change scenarios show that sediment export is substantially low from November to May, with the high sediment-export season from June to October (Figure 5).The amount of sediment export is higher than the baseline, particularly in September.This trend is the same as that of the water yield.The spatial distribution of sediment export is concentrated in the entire downstream farmland area, indicating that the sediment was washed from the upstream to the downstream area by rainfall and distributed across the agricultural land.
distribution of sediment export is concentrated in the entire downstream farmland area, indicating that the sediment was washed from the upstream to the downstream area by rainfall and distributed across the agricultural land.The amount of N export is relatively high during certain months in some of the models, for example, in December for the RCP8.5 Had GEM2-AO model; in April and October for RCP8.5 MIROC5; and in October for RCP2.6 CESM1-CAM5.The simulation results of other models were similar to the baseline.A relatively high amount of N export is observed in the baseline, with the same observations for most models occurring from April to June, and in October (Figure 6).The effect of climate change on regional N export is minute, probably because land use in the study area is simple and remains unchanged.The spatial distribution of N export is concentrated in mining areas and agricultural land in the downstream area.Nutrients are produced by surface water, point-source pollution, and agriculture [53], while the N export of agricultural land should be produced by chemical fertilizers used in the orchards [54].In the simulation, however, some N exports of the west side of the mining area might have been caused by the imported parameters of InVEST that used the values from InVEST's biophysical table.Moreover, the default value of N in the mining area is relatively high.The amount of N export is relatively high during certain months in some of the models, for example, in December for the RCP8.5 Had GEM2-AO model; in April and October for RCP8.5 MIROC5; and in October for RCP2.6 CESM1-CAM5.The simulation results of other models were similar to the baseline.A relatively high amount of N export is observed in the baseline, with the same observations for most models occurring from April to June, and in October (Figure 6).The effect of climate change on regional N export is minute, probably because land use in the study area is simple and remains unchanged.The spatial distribution of N export is concentrated in mining areas and agricultural land in the downstream area.Nutrients are produced by surface water, point-source pollution, and agriculture [53], while the N export of agricultural land should be produced by chemical fertilizers used in the orchards [54].In the simulation, however, some N exports of the west side of the mining area might have been caused by the imported parameters of InVEST that used the values from InVEST's biophysical table.Moreover, the default value of N in the mining area is relatively high.
distribution of sediment export is concentrated in the entire downstream farmland area, indicating that the sediment was washed from the upstream to the downstream area by rainfall and distributed across the agricultural land.The amount of N export is relatively high during certain months in some of the models, for example, in December for the RCP8.5 Had GEM2-AO model; in April and October for RCP8.5 MIROC5; and in October for RCP2.6 CESM1-CAM5.The simulation results of other models were similar to the baseline.A relatively high amount of N export is observed in the baseline, with the same observations for most models occurring from April to June, and in October (Figure 6).The effect of climate change on regional N export is minute, probably because land use in the study area is si Regarding changes in the P export hydrologic ecosystem service, the amount of average annual P export is relatively high in RCP8.5 GISS-E2-R and CCSM4.In other models, the amount of average annual P export is similar to the baseline.Most models and baseline values indicate that the P export is high from April to November (Figure 7).The spatial distribution of P export is only determined for agricultural land in the downstream area.The P export from agricultural land may be due to chemical fertilizers used in the orchards.
Water 2019, 11, x FOR PEER REVIEW 9 of 20 Regarding changes in the P export hydrologic ecosystem service, the amount of average annual P export is relatively high in RCP8.5 GISS-E2-R and CCSM4.In other models, the amount of average annual P export is similar to the baseline.Most models and baseline values indicate that the P export is high fro

Spatial Distribution of Hotspots Identified Using LISA
Water yield hotspots are calculated based on the average annual value.Results revealed that the regions with high water yield throughout the entire year (hotspots) were all located on the west side of the catchment area, specifically in the higher upstream region.The hotspot areas are slightly smaller than that of the baseline, both for annual and monthly calculations.Hotspot distribution for monthly water yield during the summer, moves downstream from time to time.Further, four out of five GCMs calculate almost all hotspots downstream in October (Figures S1-S6).
The average annual sediment export calculations show that the east side of the catchment area (the downstream area), and the agricultural land on both sides of the river, are areas with high sediment output throughout the entire year.The other small hotspot area is on the west side of the mining area.The spatial distribution might be attributable to the fact that while some of the sediment is carried and deposited in this area [55], the rest of the sediment is washed to the downstream area and scattered across the agricultural land on both sides of the river.There are no differences observed among different years regarding hotspot spatial distribution.
The average area of hotspots when calculated on a monthly basis, is the same as that of the baseline for the same period, and no difference was noted in the spatial distribution of hotspots for each month.
The LISA results revealed that regions with relatively high N export throughout the entire year included the mining area and forest area, which is located on the west side of the catchment area, and also the surrounding area of the downstream agricultural land on the east side of the catchment area.No difference was observed.
In the spatial distribution of hotspots for each year.The areas of hotspots in the climate-change scenarios calculated with the five models are almost the same as the baseline, and no difference was observed in the spatial distribution of hotspots for any month.

Spatial Distribution of Hotspots Identified Using LISA
Water yield hotspots are calculated based on the average annual value.Results revealed that the regions with high water yield throughout the entire year (hotspots) were all located on the west side of the catchment area, specifically in the higher upstream region.The hotspot areas are slightly smaller than that of the baseline, both for annual and monthly calculations.Hotspot distribution for monthly water yield during the summer, moves downstream from time to time.Further, four out of five GCMs calculate almost all hotspots downstream in October (Figures S1-S6).
The average annual sediment export calculations show that the east side of the catchment area (the downstream area), and the agricultural land on both sides of the river, are areas with high sediment output throughout the entire year.The other small hotspot area is on the west side of the mining area.The spatial distribution might be attributable to the fact that while some of the sediment is carried and deposited in this area [55], the rest of the sediment is washed to the downstream area and scattered across the agricultural land on both sides of the river.There are no differences observed among different years regarding hotspot spatial distribution.
The average area of hotspots when calculated on a monthly basis, is the same as that of the baseline for the same period, and no difference was noted in the spatial distribution of hotspots for each month.
The LISA results revealed that regions with relatively high N export throughout the entire year included the mining area and forest area, which is located on the west side of the catchment area, and also the surrounding area of the downstream agricultural land on the east side of the catchment area.No difference was observed.
In the spatial distribution of hotspots for each year.The areas of hotspots in the climate-change scenarios calculated with the five models are almost the same as the baseline, and no difference was observed in the spatial distribution of hotspots for any month.
The LISA results indicated that areas with relatively high P export throughout the entire year are those surrounding both the downstream agricultural land that is on the east side of the catchment Water 2019, 11, 867 10 of 20 area, and the agricultural land in the southeast.No difference was observed in the spatial distribution of hotspots for any year in the climate-change scenarios.The areas of P-export hotspots in the climate-change scenarios calculated with the five GCMs are almost the same as the baseline, and no difference was observed in the spatial distribution of hotspots for any month.

Comparison of Changes to Hydrologic Ecosystem Services within and among Years
Though climate change impacts are strongly evident, few actions have been taken because of the lingering distributional conflicts regarding the risks of shortages [1].Accordingly, this study used calculations of future changes in the water yield.Compared with the baseline, the average annual water yield in the RCP2.6 MIROC5 model increases by 4%, while in the GISS-E2-R model it decreases by 9%.In the RCP8.5 CESM1-CAM5 model, the average annual water yield increases by 8%, while in the HadGEM2-AO model it decreases by 17%.This study infers that if climate change occurs in accordance with the projected RCP2.6 scenario, then the Chinan area water yield will not change substantially (Figure 8a-d).If climate change occurs as per the RCP8.5 scenario, however, greenhouse gas emissions are presumably high, and the annual increase and decrease ratio could vary as much as two-fold.
change scenarios calculated with the five GCMs are almost the same as the baseline, and no difference was observed in the spatial distribution of hotspots for any month.

Comparison of Changes to Hydrologic Ecosystem Services within and among Years
Though climate change impacts are strongly evident, few actions have been taken because of the lingering distributional conflicts regarding the risks of shortages [1].Accordingly, this study used calculations of future changes in the water yield.Compared with the baseline, the average annual water yield in the RCP2.6 MIROC5 model increases by 4%, while in the GISS-E2-R model it decreases by 9%.In the RCP8.5 CESM1-CAM5 model, the average annual water yield increases by 8%, while in the HadGEM2-AO model it decreases by 17%.This study infers that if climate change occurs in accordance with the projected RCP2.6 scenario, then the Chinan area water yield will not change substantially (Figure 8a-d).If climate change occurs as per the RCP8.5 scenario, however, greenhouse gas emissions are presumably high, and the annual increase and decrease ratio could vary as much as two-fold.
According to the results of the five models for the RCP8.5 scenario (Figure 8e-h), the annual water yield of the CCSM4 model was equivalent to that of the baseline, and the results of the CESM1-CAM5, GISS-E2-R, and MIROC5 models all revealed increased trends.However, the annual water yield of the HadGEM2-AO model decreased by 17%.Some studies have proposed that HadGEM2-AO might underestimate water yield in the Eastern Pacific region (subtropical zones) if used to simulate past climates for comparison with observation data [56].In this study, simulation results also revealed a similar model characteristic in that the average monthly water yield exhibits the lowest values.Therefore, increases in the annual water yield hydrologic ecosystem service might be a possible future change (Tables S1-S4).According to the results of the five models for the RCP8.5 scenario (Figure 8e-h), the annual water yield of the CCSM4 model was equivalent to that of the baseline, and the results of the CESM1-CAM5, GISS-E2-R, and MIROC5 models all revealed increased trends.However, the annual water yield of the HadGEM2-AO model decreased by 17%.Some studies have proposed that HadGEM2-AO might underestimate water yield in the Eastern Pacific region (subtropical zones) if used to simulate past climates for comparison with observation data [56].In this study, simulation results also revealed a similar model characteristic in that the average monthly water yield exhibits the lowest values.Therefore, increases in the annual water yield hydrologic ecosystem service might be a possible future change (Tables S1-S4).
Most other studies have only simulated annual changes in water yield [6][7][8]22,57], which does not reflect seasonal and inter-annual variations [58].This study used monthly data to conduct more detailed simulations.The results indicate that in the RCP2.6 scenario, the total water yield of the high-flow period calculated by all the models did not change considerably when compared with the baseline.However, the total water yield of the low-flow period was significantly lower than that of the baseline.In the RCP8.5 scenario, the increase in the high-flow period and decrease in the low-flow period are greater than in the RCP2.6 scenario, which is the same as average annual water yield results.Therefore, the simulation results for total annual water yield all increased in both the RCP2.6 and RCP8.5 scenarios, and the increases all occurred in the summer.The water yield results were expected to be low in the winter, while the simulated water yield results were even lower than that of the baselines.The fact that water yield increase occurs only in the summer will certainly bring attention to the difficulty of storing water, and complicate water-resource allocation.Since Chinan has long been known as the most water-deficient catchment area in Hualien [59], the predicted inter-annual variation of hydrologic ecosystem service is disconcerting.
Under the RCP2.6 scenario, water yield in the five models decreased mostly in February and May and increased mostly in July.Therefore, the variation of both water yield and the baseline during these months exhibited the greatest change.In the RCP8.5 scenario, water yield decreased mostly in January and April, and increased mostly in July and October, which reveals an intensified change in the high-flow and low-flow periods (Figure 9a1,a2).Specifically, the gap between the high-flow and low-flow periods is substantial (Figure 9b1,b2).
Regarding the comparison of annual and average monthly water yield, the variation in monthly average increases in all of the models.The highest variation in water yield is 8%, observed in RCP8.5 CESM1-CAM5.The variation in the average monthly water yield increases by 45% (June) and decreases by 60% (April).The average annual water yield of RCP8.5 HadGEM2-AO decreases by 17%.The variation of average monthly water yield increases by 19% (August) and decreases by 88% (January), indicating that using the average monthly value to simulate the amount of water yield produces a larger variation in the future.The aforementioned results also prompted another discussion.The substantial decrease in water yield during the low-flow period is a warning for water-resource management [60].The simulation results of RCP8.5 HadGEM2-AO demonstrate the lowest average annual water yield, and the average monthly water yield notably reduces in the low-flow period.The HadGEM2-AO simulation results for the Eastern Pacific region seem to underestimate the amount of rainfall [56].There are also indications of the possibility of drastic climate change that we should seriously consider.Whether or not the current adaptation design of water-conservation engineering is sufficient for coping with water depletion in the future low-flow period, is also a topic worth consideration.
The general trend and variation pattern in the difference between average annual sediment export and the baseline is similar to that of water yield as demonstrated above [61,62].The average annual sediment export increases by 11% in the RCP8.5 CESM1-CAM5 model and decreases by 26% in the RCP8.5 HadGEM2-AO model.The difference between the simulation results of N and P export and the baseline does not exceed 1%.Climate change has relatively little effect on these two hydrologic ecosystem services.
The average monthly sediment export also exhibits a sharp increase.Compared with the baseline, the average sediment export of the RCP8.5 CESM1-CAM5 model increases by 11% annually and by 70% in June.The average sediment export of the RCP8.5 HadGEM2-AO model decreases by 26% annually and by 79% in January.Furthermore, the increase in sediment export primarily occurs from June to October, while four out of five models' sediment export simulation results show that a number of monthly yields increase by more than 50%.During this period of flood season, rainfall is relatively high and the sediment export hydrologic ecosystem service surges, which could potentially increase the erosion of gullies and river banks.Excessive soil loss might also result in reduced nutrient strength and water-storage capacity [42].

Suitability of Various Circulation Models for Small Areas
Selecting an appropriate climate model is the primary challenge of conducting research on the effects of climate change.A total of 41 GCMs were adopted in the IPCC Fifth Assessment Report [51].Because the horizontal grid resolution exported from a GCM is approximately 150-300 km, downscaling was necessary to obtain a higher resolution.Otherwise, using boundary conditions of the regional climate model would be necessary to obtain information with higher resolution [63].Therefore, uncertainty included in a number of parameters used for climate model inherently provides uncertainty while assessing climate change.Poor results are often observed in studies that use climate-model simulations.Numerous GCMs have been found to possess insufficient features to simulate regional climate characteristics [64].For example, mountains are a problem that has not been solved in all models [65] such that some studies have not considered the effects of local terrains [66] on water yield.These insufficiencies increase the number of discrepancies between simulation and observation results.For example, by comparing the simulation results of MIROC5 with those of other GCMs in this study, we find that water yield from May to September was relatively high.In particular, the values of water yield in June (RCP2.6),July (RCP8.5),and August (RCP8.5)were the highest among all models, while those in May (RCP2.6 and 8.5), June (RCP8.5),July (RCP2.6),and September (RCP8.5)were second highest among all models, which implies the occurrence of heavy rainfall during the summer.Model performance of the summer monsoon and the annual circulations of temperature and rainfall in Southeast Asia were determined to be satisfactory, however, an assessment of the summer monsoon and summer monsoon indices has revealed bias [67].This is also consistent with the result observed in this study.And yet, inconsistencies with other models can also been found.For example, a study that assessed the similarities among circulation models using climate-model genealogy revealed that CESM1-CAM5 and CCSM4 exhibited high similarity [52].However, in the present study the simulation results of these two models were obviously different with respect to water yield changes.Furthermore, CESM1-CAM5 demonstrate a greater response than CCSM4 to all future climate parameters, including higher temperature, greater precipitation, and changes in sea-level pressure [68].In this study, the annual water yield of RCP2.6 CESM1-CAM5 was lower than that of CCSM4; in the RCP8.5 scenario, the annual water yield of CESM1-CAM5 was higher than that of CCSM4, and the average monthly water yield was not substantially higher than that of CCSM4.As a whole, some model characteristics are not necessarily applicable to small areas.
Taiwan is geographically located at the intersection of the continental high-pressure zone and subtropical anticyclone, which is also at the lower edge of the monsoon gyre during the summer.Its rainfall is affected by complicated mechanisms [66].Therefore, the selection of a suitable circulation model has long met challenges in the research of climate change in Taiwan [63].One study suggests that the prediction of precipitation should be a priority function when selecting a GCM for Taiwan [52], and another asserts that the results of most models simulating the total annual rainfall of summer monsoons in East Asia are satisfactory [64].A comprehensive assessment of various model performances (e.g., the similarity of the climate zone, climate characteristics, and precipitation timing to those of Taiwan) might aid in the establishment of a more rigorous selection indicator.
As mentioned in the methodology section, five GCM models (i.e., CESM1-CAM5, GISS-E2-R, CCSM4, HadGEM2-AO, and MIROC5) were selected based on the applicability rankings in different climate zones of Taiwan.Among them, HadGEM2-AO has been verified to perform better in terms of observation data of weather stations in Taiwan [52].CESM1-CAM5 ranks at the top when only the weather stations in eastern Taiwan were taken into account.In addition, its performance for mountain areas and offshore islands is better.If we take the average from our five simulations as our consideration basis, the total annual water yield will increase and the difference between annual high-flow and low-flow periods will become more noticeable.Four out of five applied GCMs are able to predict both of the situations above, while HadGEM2-AO can only reveal the latter.Among the four, the CESM1-CAM5 predicted the two climate characteristics mentioned above.This finding is consistent with the conclusion of a previous study [52] stating that CESM1-CAM5 is the most suitable GCM for application in eastern Taiwan.

LISA Results Helped Identify the Spatial-Distribution Characteristics of Hotspots
This study used the hydrologic ecosystem service functions of the InVEST simulation as LISA input data to identify the hotspots of hydrologic ecosystem service functions through a weighted calculation.The analysis was conducted to determine the optimal space for a key hydrologic ecosystem services protection area [55].Cells identified as hotspots indicate that those regions require attention to the decline in hydrologic ecosystem service functions in a climate-change scenario, and that pre-adjustments should also be made.Maps and diagrams transform the scientific results of assessing hydrologic ecosystem services into information that can be understood by policy makers, stakeholders, and the public [9,69].This is a very important step for incorporating hydrologic ecosystem services into decision-making procedures.
Compared to the baseline, the hotspots computed by LISA using simulated data from five GCMs and two RCP pathways demonstrate very insignificant change.As previously presented [61], our results show that the same geographic relation identified as hotspots are generally concentrated in the high elevation areas.However, the monthly hotspot distributions show a significant change in the summer time, and pronounced in October (Figure 10), as distributions move to downstream areas.This phenomenon reminds us that even if the annual hotspot remains the same as the baseline, summer water yield may still substantially change its geographical distribution.Summer marks the rainy season in Taiwan, and the major water supply season as well.Such a hydrologic ecosystem service change draws our attention since water intake locations may have to be reconsidered.However, spatial heterogeneity of land use and distributed precipitation rates throughout areas with various hydrological properties may cause additional hydrologic ecosystem services as a result of these hydrological processes [7,8].In this study, the areas with the highest hotspot concentrations have higher fragmental levels due to hydrological conditions under RCP8.5 when compared to those under RCP2.6, with the except of the CCSM4 and MIROC5 GCM models.In various climate-change scenarios, the average annual sediment export hotspot areas obtained from the ten scenarios (i.e., five models under two RCPs) are the same as that of the baseline (Figure S7).The hydrologic ecosystem service function of sediment export simulated using InVEST revealed that climate change increased sediment export from June to October.The spatial distribution of hotspots calculated by LISA exhibited no difference for any month, though the amount of sediment export increased during the flood season [60,70].Therefore, water-resource managers are advised to pay closer attention to the hotspot areas identified by LISA during this time period [33,71].
The N and P hotspots located in agricultural land in the study area distributed mainly along the riverine areas [26,54,57].For the spatial distribution of N export, the average hotspot area calculated with the ten models is similar to that of the baseline.Regarding the spatial distribution of average monthly N export hotspots, some models revealed that the hotspot area is smaller from April to In various climate-change scenarios, the average annual sediment export hotspot areas obtained from the ten scenarios (i.e., five models under two RCPs) are the same as that of the baseline (Figure S7).The hydrologic ecosystem service function of sediment export simulated using InVEST revealed that climate change increased sediment export from June to October.The spatial distribution of hotspots calculated by LISA exhibited no difference for any month, though the amount of sediment export increased during the flood season [60,70].Therefore, water-resource managers are advised to pay closer attention to the hotspot areas identified by LISA during this time period [33,71].
The N and P hotspots located in agricultural land in the study area distributed mainly along the riverine areas [26,54,57].For the spatial distribution of N export, the average hotspot area calculated with the ten models is similar to that of the baseline.Regarding the spatial distribution of average monthly N export hotspots, some models revealed that the hotspot area is smaller from April to October (wet season) and larger in November to March (dry season) of the next year (Figure S8).Increasing temperatures could reduce soil erosion due to cultivating period [26,72].This indicates that the area of relatively high N export increases during the low-flow period.On the other hand, because the area is mostly located in the upstream zone, the effects of N export on water quality and cost of water cleaning should be heeded [73].The average annual P export is almost the same as that of the baseline, and no difference is observed in the area or spatial distribution of hotspots for any month (Figure S9).These areas coincide with farmland; thus the P export can be attributed to fertilization [26,55,72].However, this study also finds that a threshold for nutrient loads may be reached [7,74].
Policy-making is essential in mitigating the effects of climate change to ensure suitable adaptation under variate uncertainty levels [7,8,22,26].This study evaluated and identified hydrologic ecosystem service hotspots under RCP2.6 and 8.5 scenarios and concluded that an effective adaptation strategy is one that is suited to deal with climate change on hydrologic ecosystem services in the study area [22,26].Our analysis shows that by integrating hydrologic ecosystem services with hotspots analysis, this approach provides relevant spatiotemporal information on ecosystems and hydrologic ecosystem services at appropriate scales, which is suitable for decision-making under climate change.In this study, we used equal weighting to identify hotspots of hydrologic ecosystem services to select strategies under climate change RCP2.6 and 8.5.It is essential to study the tradeoffs between ecosystem services, but a suitable systematic planning framework [6] (such as the approach used in our study) offers methods to effectively identify valuable synergies for hydrologic ecosystem services conservation and management.However, further studies should consider trade-offs among hydrologic ecosystem services under all possible scenarios.Moreover, to balance ecological integrity, it is suggested that social utility and urban development should be taken into account [22] in hydrologic ecosystem services management under climate change at the watershed scale.
In this study, the combination of various models, including InVEST and LISA, was used to quantitatively identify hydrologic ecosystem service hotspots (i.e., water yields, sediment export, Nitrogen and Phosphorous exports) under various climate change scenarios.Unlike other relevant studies [75,76] on climate changes, our approach with combines InVEST and LISA has the advantage of identifying spatial autocorrelations among neighboring cells for ecosystem services valuation [77], and for the planning and management of spatial adaptation of water resources under climate change.Annual and monthly hydrologic ecosystem service hotspots changes can be analyzed to show intra-annual fluctuations and annual variation of the services to enhance hydrologic ecosystem services management and planning for climate change impacts.

Conclusions
The impact of climate change on hydrologic ecosystem services should not just be considered annually, but also monthly to show intra-annual fluctuations and annual variation of the services.Under five climate-change scenarios proposed in AR5, it is apparent that until 2035, the annual average value of water yield and sediment export in the study area are both more abundant in wet seasons but more deficit in dry seasons compared to the base-line period from 1996 to 2009.Moreover, the intra-annual fluctuations have increased significantly due to climate change, thus that considering the maximum extent of the impact that hydrologic ecosystem services may face is necessary.However, in general, the annual total water volume will increase while intra-annual fluctuations between the wet and dry seasons will become even more significant through 2035.Accordingly, decision makers should adopt water conservation plans that secure water supply for future water demand.
Hydrologic ecosystem service values are spatially correlated.Unlike other relevant studies, our LISA analysis of the four hydrologic ecosystem services in the study catchment's hotspots shows that through 2035, the annual average water-yield for hotspots generally remained the same as the base-line period, though water-yield hotspots were relocated in downstream areas during summer.Taking this information into account, sustainable planning of the water intake point setting is necessary, as well as planning to increase the summer water intake in lower reach hotspot concentrated areas.Five GCMs were adopted in this study and their performance within the small study area was compared.Results indicate that CESM1-CAM5 is particularly able to predict two climatic characteristics in Taiwan, demonstrating its suitability.As the case in this study, common GCMs are generated for the purpose of large-scale applications.When applied to a small area, however, simply using statistical downscaling methods may result in low predictability due to unknown local effects.Therefore, choosing models suitable for specific areas is necessary to better undertake comparisons among GCMs.This study's method is successfully applied in the study area and can be extended to other regions in the future.Our proposed approach yields information under climate change that can be considered in water resources and watershed planning and management.
of poor weather conditions on the Hualien-area hydrologic ecosystem services in the context of climate change.

Figure 1 .
Figure 1.The study area.(a) Topography and river systems in the Hualien area [40]; (b) land use of the Chinan area.

Figure 1 .
Figure 1.The study area.(a) Topography and river systems in the Hualien area [40]; (b) land use of the Chinan area.

Water 2019 ,
11, x FOR PEER REVIEW 4 of 20
study, the baseline refers to 1996-2009, and the climate change scenario refers to 2035.Historical weather data used in this study was collected from the Central Weather Bureau and Water Resources Agency and was then downscaled to a 40 × 40 m cell resolution with the Kriging method.

Figure 4 .
Figure 4. Spatial distribution of water yield in January and September.

Figure 4 .
Figure 4. Spatial distribution of water yield in January and September.

Figure 4 .
Figure 4. Spatial distribution of water yield in January and September.

Figure 8 .
Figure 8. Four hydrologic ecosystem services, water yield (a,e), sediment export (b,f), N export (c,g) and P export (d,h) under RCP2.6 (a-d) and RCP8.5 (e-h) scenarios.Note: The six columns along the x-axis are the baseline and five general circulation models (GCMs), while the y-axis is the annual average output from InVEST.

Water 2019 ,
11, x FOR PEER REVIEW 15 of 20 have higher fragmental leve

Figure 10 .
Figure 10.Hotspot spatial distribution of water yield in October.

Figure 10 .
Figure 10.Hotspot spatial distribution of water yield in October.
using parameters associated with RCPs above.Through the efforts of Coupled Model Inter-comparison Project Phase 5 (CMIP5), sponsored by

Table 1 .
Hydrologic ecosystem services calculated for different climate change scenarios.