Spatial conservation areas for water yield hydrological ecosystem services with their economic values effects under climate change: a case study of Teshio watershed located in northernmost of Japan

Ecosystem services (ESs) sustain and fulfill human well-being through provisioning, regulating, and cultural services that are formed based on various conditions, processes, and components of natural (such as forest and wetland) and artificial (such as farmland and paddy field) ecosystems (Wallace 2007). Among these ESs, hydrological provision of water yield is important in economic and agricultural development. It supplies energy through hydroelectric dams, water for drinking and agricultural irrigation. Hydrologic flow and cycle fluctuate greatly on the spatial scale, hence, it is often difficult to depict their characteristics on local and regional scales in heterogeneous environments (Nepal 2016). The heterogeneity of water yield is mainly influenced by meteorological factors such as precipitation, temperature, wind velocity, solar radiation, and others (Li et al. 2015). Therefore, assessing the effects of climate change on catchment water yield (including surface, lateral, and groundwater) is crucial for the future management of water supplies. However, undertaking these assessments on a regional scale with sufficient resolution, which is useful to water resources engineers and planners, is a challenge. The intensity and characteristics of the impact of climate change on water resources and freshwater ecosystems vary significantly from region to region. Therefore, the protection of water yield ESs is an important criterion in sustainable water resource management. This is predominant in cases when human water consumption is in competition with environmental water requirements. Furthermore, the impact of global warming can influence the availability and demand for water resources. This can have a further impact on the availability of water required to sustain ecological functions (Kirshen et al. 2005; Barron et al. 2012). Coupled with increasing water demand, it is likely to result in large increases in the number of people at risk of water scarcity. The increase in population, economic growth, and industrialization has resulted in higher water demands for household usage, irrigation, and hydroelectricity production. This has led to the excessive exploitation of water resources in many parts of the world (Fan & Shibata 2014). Water deficits have exacerbated in regions experiencing reduced rainfall and rising temperatures because of anthropogenic global warming (Zölch et al. 2018).

It is important for common conservation practice to protect from unreasonable exploitation of ESs. Currently, regular markets do not reflect ES values. Therefore, effective planning and management tools are required to regulate activities with potential negative impacts on ESs (Fan & Shibata 2014; Pan et al. 2017). In particular, it is widely recognized that different conservation planners need to account for economic costs and potential hydrological ES loss when designing an ecosystem reserve network. This has led to the widespread adoption of the systematic conservation planning approach. The conservation approach is a target-driven process that aims to identify networks of priority areas for ensuring the representation and long-term persistence of hydrological ES. Setting targets helps to increase transparency and measure progress. Also, it allows socio-economic data to be included in the planning process without influencing or endangering conservation goals. In the valuation of hydrological ES, cost information is used. Cost information is the amount of money made through each water supply in each planning unit (Gosselink et al. 1974; Naidoo et al. 2006). Much effort in systematic conservation planning of water resources has gone into devising measures and algorithms which are efficient at capturing the economic and ecological importance of different candidate areas. While cost is at the theoretical heart of the complementarity based analyses, little empirical attention has been paid to explicitly incorporate the cost into conservation planning. There are also few studies on integrating economic costs into systematic conservation planning of water resources across the watershed under climate changes. These studies did not spatially quantify the economic values of hydrological provision ESs and integrated these economic values into the systematic conservation planning of water resources to achieve different conservation targets (e.g. representing 10–90% of potential maximum water yield provision ESs) under climate changes. Furthermore, they also did not analyze how the economic value of hydrological provision ESs impacts the spatial hotspots for conservation areas for water resources, according to the relationship between water supply and demand under climate changes. It is inevitable to address this gap by utilizing the first estimate of economic costs for water resources conservation across the entire watershed under climate changes to examine the consequences of incorporating estimates of economic costs into conservation planning. Thus, this ecosystem-based reserve network can be designed so that it meet targets, while also minimizing impacts on stakeholders and increasing the likelihood of successful implementation.

ES-based management through economic values of ESs proved to be an exercise in system thinking and helped to guide our initial understanding of ESs panning (Costanza et al. 1997, 2017; Naidoo et al. 2006; Fan & Shibata 2014). From this, we gained perspective on system disturbances, drivers, alternate regimes, and cross-scale interactions. We were also better prepared to identify the conservation features which are in most need of protection. Spatial priority areas for ESs with replacement cost layers can be used for evaluating proposed conservation areas, for evaluating the loss arising from some high-quality area becoming unavailable for conservation (Cabeza & Moilanen 2006; Leathwick et al. 2008; Moilanen et al. 2008). In this study, we developed a systematic conservation framework of water resources through integrating the hydrology model into a systematic conservation model. This framework comprehensively considers the spatial heterogeneity of water yield affected by topography, meteorology, vegetation, and soil driving factors. In addition, it considers the spatial pattern of economic values of water yield, which is determined by locations of power plants, and spatial distributions of populations and agriculture. Therefore, this paper has two major goals: (1) to identify spatial conservation areas where conservation effects should be directed for hydrological ESs provision with different conservation target levels under climate changes across watershed scale; and (2) to contrast the impacts of the costs of water resources on spatial conservation areas for hydrological provision with different target levels under climate changes. In response to these objectives, this paper attempts to demonstrate that the economic values (cost) of hydrological provision ESs provided by ecosystems across the watershed should be involved in water resources planning and management. This will balance the contradiction between water conservation and exploitation under climate change.

This study was conducted at the Teshio River watershed, northern Hokkaido, in northernmost Japan, which is located at 44.33°N, 142.25°E (Figure 1(a)). This river is the fourth longest (256 km) in Japan, originating from Mount Teshio and flowing into the Sea of Japan. The agricultural land, population and three hydropower plants concentrate on the upper and middle watershed, implying that there are close relationships between the water yield and anthropogenic activity (the water yield is used for irrigation, resident daily usage and hydropower electricity production). The interaction between humans and water resources is more active in the upper and middle watershed than that in lower places. Therefore, we focused on the upper and middle watersheds in this study. The catchment area of the upstream and middle of the Teshio River watershed is 2,908 km². Average water flow is 138.6 m³ s^–1, with a maximum of 1,278 m³ s^–1 in April from snowmelt and a minimum of 19.8 m³ s^–1 in February during the mid-snowpack season at Bifuka monitoring station (Figure 1(a); Katsuyama et al. 2009). Approximately 78% of the catchment is covered by forest, categorized as cool-temperate mixed forest, including deciduous broadleaf and evergreen coniferous species with a dense understorey of Sasa dwarf bamboo (Ileva et al. 2009). Other land uses are mainly farmland and paddy fields, with area percentages of 13 and 4%, respectively. The remaining 5% land use was urban and water body. The farmland and paddy fields are mainly distributed on both sides of the main Teshio channel and on the southwestern watershed (Figure 1(b)).

The soil is dominated by brown forest soil (Cambisol; IUSS Working Group WRB 2006); others are gray lowland soil (Gleyic Fluvisols; IUSS Working Group WRB 2006), brown lowland soil (Haplic Fluvisols; IUSS Working Group WRB 2006), grey soil (Gleyic Fluvisols; IUSS Working Group WRB 2006), and peat soil (Histosols; IUSS Working Group WRB 2006) (Figure 1(c)). The low values of field water capacity occur on both sides of the main Teshio channel and on the southwestern watershed (Figure 1(d)). The slope ranges from 0 to 83.8% with an average of 14.5% (Figure 1(e)). The climate of the area is characterized as cool and relatively humid, with an annual average temperature of 5 °C and precipitation of approximately 1,000 mm. There is snowfall from November through March, contributing to almost 40% of the annual precipitation measured at Bifuka station. There are spatial differences of precipitation across the watershed (Figure 1(f)). Precipitation in mountainous areas at higher altitudes (southeastern, northern and western parts of the watershed) is substantial, especially in October. Water discharge is generally characterized by one main peak in spring during snowmelt, and a second smaller peak, usually between late summer and early autumn (Ileva et al. 2009; Katsuyama et al. 2009).

The overall analytical framework consists of developing climate change scenarios, modeling of water yield hydrological ES, addressing the spatial pattern of economic values of water yields, and simulating spatial conservation areas for water yield hydrological ESs with their economic value effects under climate change (Figure 2). We developed the multiple climate change scenarios generated by a general circulation model (GCM) for the study watershed. Following the development of change scenarios for climate, we simulated the hydrological ESs using a watershed hydrology model (Soil and Water Assessment Tools, SWAT) under climate change scenarios at the watershed scale. We simulated the spatial conservation areas for water yield hydrological ESs using a systematic conservation model (Marxan model). We then analyzed the spatial configuration of conservation areas for water yield hydrological ESs with their economic value effects under climate change through fitting a regression model, Jaccard's index and Spearman correlation.

Here, and are the current and future mean monthly precipitation, respectively. and are simulated mean monthly precipitation under the current and future scenario climate conditions, respectively. Climate changes estimated, based on the IPCC-DCC Fourth Assessment Report SRB1 scenario (IPCC 2007), for the study site are given in Table 1. According to the above procedures, we calculated average data for short-, mid- and long-term climate changes for the following model simulation. The generation method of daily climate data for that simulation is described in the next section.

SWAT is a hydrologic and water quality model developed by the United States Department of Agriculture Agricultural Research Service (Arnold & Allen 1996; Abbaspour et al. 2007; Schuol et al. 2008). SWAT was developed to predict the impact of land use management practices on water, sediment, and chemical yields in large complex watersheds with varying soils, land use and management conditions over long periods of time. In this study, SWAT was used to assess the spatial and temporal distribution of water yield with varying soils, land use and climate conditions within basins. The model is a continuous-time, spatially distributed simulator of the hydrologic cycle and pollutant transport at catchment scales, and runs on a daily time step. The watershed is divided into multiple sub-watersheds, which are then divided into units of unique soil, land-use and slope characteristics called hydrologic response units (HRUs). These HRUs are defined as homogeneous spatial units characterized by similar geomorphologic and hydrologic properties. The overall hydrologic balance is simulated for each HRU, including canopy interception of precipitation, partitioning of precipitation, snowmelt water, and irrigation water between surface runoff and infiltration, redistribution of water within the soil profile, evapotranspiration (ET), lateral subsurface flow from the soil profile, and return flow from shallow aquifers.

The major input data of running the SWAT model included climate data, a terrain map, soil properties and a land-use map. The climate data were obtained from the Japan Meteorological Agency (JMA). The SWAT model requires meteorological data such as daily precipitation, maximum and minimum air temperature, wind speed, relative humidity, and solar radiation. We used Sharpley and Williams's WXGEN weather generator within the SWAT to generate daily temperature and precipitation data for the target climate scenarios (Sharpley & Williams 1990). To obtain statistical weather parameters of stochastic processes and distributions in WXGEN (Williams & Griffiths 1995; Neitsch et al. 2011), we used observed climate data from 1997 to 2009, because continuous daily data for precipitation and temperature were available for this period.

The DEM was acquired from the Japanese Geographic Survey Institute (GSI; 50 × 50 m). The land-use map from 2006 and vector soil data were obtained from the National Land Information Office in the Ministry of Land, Infrastructure, Transport and Tourism (MLIT; 50 × 50 m).

SWAT was initially set up directly for streamflow during 2001–2002, 2003–2007 and 2008–2009 these three time periods were taken as warm-up, calibration and validation periods, respectively. The coefficient of determination (R²) and Nash–Sutcliffe index (NSI; Nash & Sutcliffe 1970) were used to evaluate model performance. Sensitivities are usually defined by computing partial derivatives of the objective functions with respect to the parameters. The objective functions are usually some lumped measures such as the total mass export, sum of square error between modelled and observed values or sum of absolute errors. A sensitivity index, S, can be calculated for a small change of one given parameter while the other input parameters are held constant. A sensitivity analysis was then performed using Latin hypercube (LH), a one factor-at-a-time method (OAT) (Arnold & Allen 1996; Abbaspour et al. 2007). LH could subdivide the distribution of each parameter into N strata with a probability of occurrence equal to 1/N. This approach results in N non-overlapping realizations and the model is run N times. The OAT consists of repetitions of calculating the sensitivity index S whereby the derivatives are calculated for each parameter by adding a small change to the parameter. The change in model outcome can then be unambiguously attributed to such a modification by means of an elementary effect, S_i. The LH-OAT performs LH sampling followed by OAT sampling. It starts with taking N LH sample points for N intervals, and then varying each LH sample point P times by changing each of the P parameters one at a time, as is done in the OAT design. The method operates by loops and each loop starts with a LH point. Around each LH point j, a partial effect (sensitivity index) S_i,j (sensitivity index) for each parameter is calculated, therefore a loop requires P + 1 runs. A final effect is calculated by averaging these partial effects of each loop for all LH (N loops). The method is very efficient as the N intervals in the LH method require a total of N*(P + 1) runs. The final effects can be ranked from the largest effect to the smallest effect.

The calibrated and validated SWAT model was used to simulate annual average water yield hydrological ESs under climate change scenarios. The water yield predicted by the SWAT is defined as the net water amount leaving the sub-basin that contributes to streamflow, which is directly obtained from SWAT results. Water yields on land include surface, lateral, and groundwater in each HRU.

The economic value for ESs can be classified into direct and indirect economic values (Costanza et al. 1997). In this study, we focused on the direct economic value of hydrological provision ESs as financial cost. Water in the upper and middle Teshio River watershed is mainly utilized for hydrologic power, residents, and paddy field irrigation. In this study, we addressed the spatial pattern of economic values of hydropower electricity production, irrigation for crop production, and resident water use (as payoff fees). The water demands for the hydropower electricity production, crop production, and resident water usage are annual average time scale. The explicit calculation processes and results are obtained from Fan & Shibata (2014) (Figure 3).

1. Economic values of hydropower electricity production. We assumed the unit price of electricity at 20.7 yen per Kw h, based on the current price in Japan from ‘International Energy Statistics’. There are three hydroelectric power plants (Pon-Teshio, Iwaonai and Niupu-gawa) in the study watershed. The total cost of electricity production at each plant was estimated using annual electric production. The economic value of water yield for electric power generation in each HRU of the upper watershed for each power plant was allocated according to the relative contribution of that yield in each HRU to its total yield, as described in Equation (4) (electricity economic value):

   

   (4)
2. Economic values of resident water use (as payoff fees). The unit price of resident daily water usage was assumed as 155 yen per cubic meter, based on the current price in the city of Sapporo. The economic value in each HRU was calculated using this water price and population in the HRU, as described by Equation (5) (resident economic value):

   

   (5)
3. Economic values of irrigation for crop production. The unit price of rice was assumed as 235 yen per kilogram, based on the current price from the Ministry of Agriculture, Forestry and Fisheries. Annual average rice and water yields in each HRU were obtained from the SWAT model. We aggregated the HRUs into a sub-watershed level, according to the location of irrigation stations for paddy fields. The predicted irrigation economic value (Equation (6)) in each aggregated sub-basin was assigned to each HRU according to its relative water yield in the sub-basin. The locations of resident daily water supply and irrigation station in this study watershed were obtained from the National Land Information Office in the MLIT:

   

   (6)

The economic value of water yield hydrological provision ES obviously had varying spatial patterns. Owing to the geographical locations of Iwaonai, Pon-teshio and Niupu-gawa hydropower stations, the larger economic values of water yield used in hydrologic electricity production occurred in the southeast and north of the study watershed (Figure 3(a)). Because the population was mainly distributed in the cities of Shibetsu and Nayoro, the greater values related to residential water usage appeared in the west and east (Figure 3(b)). Due to spatial distributions of paddy fields and irrigation stations, the larger values for paddy field irrigation were in the west, east and north (Figure 3(c)). Actually, the water yield used for hydropower generation was also used for the irrigation of paddy fields and residential water usage in terms of the spatial patterns of their economic values (Figure 3(a) and 3(c)). Especially, the water yield originated from the northern study watershed was not only used for hydropower electricity production but also for the irrigation of paddy fields. The spatially distributed pattern of total economic value of water yields was almost identical to those of irrigation (Figure 3(d)). The economic values of water yield hydrological provision ESs for hydropower electricity production, residential usage and paddy field irrigation were 0–656 million, 0–203 million, and 0.221–1,756 million yen, respectively. The total economic value of these three hydrological provision ESs ranged from 0.246 to 1,960 million yen.

Systematic conservation planning is a widely used approach for designing protected area systems and ecological networks (Ball & Possingham 2001; Harwood & Stokes 2003). The MARXAN model is one of the systematic conservation models, which generally involves dividing the planning region into a series of planning units and using optimization algorithm to select outcomes of these units that meet specified conservation targets whilst minimizing conservation costs (Airame et al. 2003; Cook & Auster 2005). The MARXAN model allows the stakeholders to consider many aspects of the conservation problem of ESs: the number and types of conservation ESs included in the analysis, a target for each individual conservation ES, how important it is to meet targets for these conservation ESs, the status of planning units and the cost of each site that could be involved in the reserve network (Possingham et al. 2000; Culver et al. 2004; White & Kiester 2008). MARXAN has been widely applied in the following aspects: prioritizing areas for land acquisition, critiquing an existing proposal reserve network, providing a means for diverse stakeholder groups to develop proposals that represent their own interests, investigating the scope and scale of possible designs for effective broad scale networks, and determining where to focus conservation efforts (Possingham et al. 2000; Chan et al. 2006; Egoh et al. 2011; Nackoney & Williams 2013). Therefore, we took the MARXAN model as an effective tool to construct priority conservation areas of water yield hydrological ESs under climate changes with different economic cost indicators of water resources.

Water yield hydrological ES used in the MARXAN model was obtained from each HRU in the SWAT model, and all calculation in the MARXAN was based on planning unit hexagons (in this study, the area size is 100 ha). The basic planning unit hexagons can produce more efficient results, defined as the extent to which priority areas minimize costs whilst also meeting conservation targets. The detailed process is as follows: (i) the spatial distributions of simulated water yield hydrological ESs of HRUs (shape file) were established based on simulated results of the SWAT model; (ii) the spatial distributions of simulated water yield ESs of HRUs (shape file) were changed into grid-based data using convention tools in the ArcGIS (grid file); (iii) the grid-based data sets of water yield ESs were then changed into planning unit hexagons data by spatial analysis tools (shape file); (iv) the planning unit hexagons data sets were taken as input data of MARXAN; (v) the MARXAN calculated optimal priority conservation areas and selection times for water yield ESs through the whole of the watershed.

Moreover, the MARXAN can consider the distribution of ESs, connectivity responses, variable emphases on the highest local quality of ESs, and cost efficiency for ESs conservation planning. Especially, the cost can effectively be used to reflect corridors that connect protected areas with one another, or biodiversity features with protected areas, by lowering the cost value of the planning units within these identified corridors. So, in the spatial analysis through Marxan, the costs, i.e. constraints, are one of the main drivers of the spatial outcome. The objective function that Marxan uses has the constraint that the cost of a selected area is the linear combination of costs of all planning units within the very same selected area. The cost of the selected area can be of any number of measures including non-monetary and monetary values. In this study, the area of planning unit, economic value of electric power production, economic value of resident water usage, economic value of irrigation for rice production, and total economic value of three water utilizations of each planning unit were taken as a measure of cost, respectively.

To understand the trade-offs of spatial conservation priority areas for water yield ESs with different economic values, we compare the similarity of protected areas and correlation of their selection times among water yield ESs with different economic values. This study measured the pair-wise similarity of best runs from each analysis for water yield ESs with different economic values using the Jaccard's index to investigate similarities in their spatial pattern (Shriner et al. 2006). The index is calculated as A/(A+B+C), where A represents planning units present in both outcomes, and B and C represent planning units present in only one of the respective outcomes. To compare the selection frequency patterns between the different outcomes of water yield ESs under different economic values, which is one aspect of the spatial patterns of the outcomes identified by MARXAN, this study also used the Spearman correlation associated with planning units selection times from water yield ESs to verify trade-offs between economic values (Warman et al. 2004).

The whole Teshio watershed was divided into 61 subwatersheds and 648 HRUs (Figures 4 and 5). Ten parameters were selected as the most sensitive by the LH sensitivity analysis (Table 2). The optimal values of each parameter obtained through calibration are also listed in Table 2. In this study, we chose the sum of squares error between the modelled and observed values as an objective function. The final effect sensitivity index value is an indicator of sensitivity to each parameter (a higher value indicates greater sensitivity). The parameter of the initial runoff curve (CN2) was most sensitive, followed by snow-related ones (SFTMP, SMTMP, SMFMX, SMFMN, and TIMP). Since the study site is in northern Hokkaido, the land surface is covered with snowpack for almost half a year, from November through March. This suggests that the winter climate increases sensitivity to snow-related parameters such as SFTMP, SMTMP, SMFMX, SMFMN, and TIMP. The distribution of the measured average annual water quantity and quality in the calibrated and validated periods was close to the simulated values, and the two distributions had a high degree of correlation, that is, R² and NSE were greater than 0.85 (Figure 6). Hence, the simulated results were credible and favorable, which indicated that the SWAT model had outstanding applicability in the Teshio watershed. The calibrated and validated SWAT model can be used to simulate water yield ESs under climate change scenarios.

The annual average precipitation increased from 958.7 mm under baseline climate conditions to 1,099 mm with long-term climate change (Figure 7). There was an obvious reduction of annual average snowfall owing to the increase of temperature in winter (Table 1). The annual average water yield across the watershed increased from 581.8 to 603.3, 623.3, and 671.2 mm for baseline climate conditions, short-term, mid-term, and long-term climate changes, respectively. The actual ET also showed an increasing trend under a warmer climate. All climate change scenarios resulted in greater monthly precipitation than that of the baseline climate conditions, especially in September, but May and October had different tendencies (Figure 8(a)). In October, all climate change scenarios reduced monthly precipitation over that of the baseline climate conditions. The maximum reduction in snowfall was from October through April for the long-term climate change scenario, in which the annual average temperature increased by 3.88° (Figure 8(b)). Snowfall decreased in November, December, and March for all climate change scenarios. Snowfall in October and April decreased to zero in the mid- and long-term climate change scenarios, and changes of snowfall in January and February were not apparent. It meant that the snow-cover-free period increased from 5 months (May–September) under baseline climate conditions to 7 months (April–October) under future climate change scenarios. Monthly average surface runoff increased under three climate change scenarios in January, February, March, August, September, November, and December, with maximum values in March (Figure 8(c)). Potential evapotranspiration (PET) and ET increased under climate change scenarios, especially during the snowmelt (March and April) and summer seasons (June and July) (Figure 8(d)). Despite changes in the seasonal distribution of ET, the change in annual average ET was relatively small. The maximum value of water yield appeared in April under the baseline climate conditions, short- and mid-term climate changes (Figure 8(e)). Under the long-term climate change, the water yield peaked in March, earlier than in other scenarios. The spatial patterns of water yield hydrological provision ESs under climate change scenarios were mostly the same. However, there is a subtle difference between the magnitudes of water yields. We took the spatial distribution of water yield under the long-term climate change scenario as an example to demonstrate its spatial characteristics and constructed a priority conservation area for water resources. This is because the largest values of water yields occurred in the long-term climate change scenario (Figure 9). The annual average water yield ranged from 437.33 to 833.65 mm H₂O. The largest values appeared in the northern, southeastern, and southwestern areas of the watershed, for example, the largest value (833.65 mm) occurred in the northern part.

There are no clear pre-existing conservation targets for water supply ES in the Teshio watershed. Therefore, we tested the different target levels in spatial prioritizing water yield ES (Egoh et al. 2010, 2011; Fan & Shibata 2016). We simulated the forest land and agricultural land priority areas for water yield ES with different maximum total amount target levels (from 10 to 90% of 671.2 mm which is the water yield under long-term climate change) under climate changes and different cost layers (Figures 10–14). Both the forest and agricultural land priority conservation areas for water yields under all climate change scenarios increased as the conservation targets increased. However, the variation trends of forest land priority conservation areas for water yields under all climate change scenarios were more regular than those of agricultural land priority conservation areas. This is because forest land is the dominant vegetation in Teshio watershed and it determines the main hydrology and ecological processes. The forest and agricultural land spatial priority conservation areas for water yields were used to satisfy conservation target levels' demands. Under all the different economic cost indicators of water resources, there was a decreasing trend from baseline climate conditions to long-term climate change scenarios. Hence, the marked curves of agricultural and forest land spatial priority conservation areas for annual average water yield under baseline climate conditions were always above the climate change scenario curves. The marked curves of agricultural land spatial priority conservation areas for annual average water yield without cost layer were downward convex under all climate change scenarios. Particularly, there were some typical tipping points that occurred in the agricultural land priority areas for annual average water yield with the cost layer calculated from hydropower electricity production. The marked curves of agricultural land spatial priority conservation areas for annual average water yield with cost layer calculated from hydropower electricity production were firstly downward convex until the conservation target was up to 40% of 671.2 mm. It then changed into upwardly convex until the conservation target was up to 60% of 671.2 mm. However, all marked curves of agricultural land spatial priority conservation areas for annual average water yield with cost layers were upwardly convex. These resulted from economic cost layers calculating the resident water usage, irrigation for rice production, and the summation of the three economic values of water resources. The shape and varied trends of these curves were almost similar. Under all climate change scenarios, all marked curves of forest land spatial priority conservation areas for annual average water yield were linearly related to the conservation targets with and without cost layers. For the cost from area, the marked curves of forest land spatial priority conservation areas for annual average water yield began to separate from each other when the conservation target was up to 40% of 671.2 mm. For the three different cost layers calculated from water resources, all marked curves of forest land spatial priority conservation areas for annual average water yield began to separate from each other when the conservation target was up to 20% of 671.2 mm.

We compared the results using fitting regression analysis (Figures 10–14; Tables 3 and 4) (model standard errors of all fitting regression equations were less than 0.1, all the fitting coefficients fell into 95% confidence intervals). This was carried out to further understand the quantitative relationships between the spatial priority conservation areas for annual average water yield and conservation target levels under climate changes with and without cost layers. The linear correlated relationships between forest land spatial priority conservation areas for annual average water yield and conservation targets were stronger than those between agricultural land spatial priority conservation areas. The linear correlated relationships between agricultural land spatial priority conservation areas for annual average water yield with cost from irrigation for rice production and total cost and conservation targets decreased from baseline climate conditions to long-term climate change (Table 3). The linear correlated relationships between forest land spatial priority conservation areas for annual average water yield with cost from resident water usage and cost from irrigation for rice production and conservation targets decreased from baseline climate conditions to long-term climate change (Table 4). The agricultural land spatial priority conservation areas for annual average water yield were more sensitive to conservation targets compared to forest land spatial priority conservation areas.

We simulated the forest land and agricultural land spatial priority conservation areas for water yield with different maximum total levels (from 10 to 90% of 671.2 mm) under the long-term climate change scenario. In this research, we found that there were obvious tipping points in the forest and agricultural land spatial priority conservation areas for annual average water yield with a 20 and 40% maximum total amount target level (20 and 40% of 671.2 mm). Therefore, we chose the 20 and 40% of the maximum total amount values of water yield ES as the target levels for the spatial priority conservation planning to explore the target on outcomes. These determined 20 and 40% targets could reflect spatial characteristics of priority conservation, that is, especially when the aim is to conserve hydrological provision ES across the watershed involving more forest land. The selected typical target levels could also be used as the conservation thresholds for avoiding water resources conservation and water exploitation, and ensure highly efficient water utilization in the study watershed. The spatial priority conservation identified for water yield ES with 20% of the maximum total amount target level with different cost layers under the long-term climate change scenario ranged from 15.83 to 20.46% of the landscape (Figure 15(a)). The spatial priority conservation areas for water yield with area cost were mainly distributed in the north, southwest, and southeast of the study site. These priority conservation areas were related to the spatial pattern of water yield (Figure 8). The districts with higher water yield were preferentially involved into priority areas for meeting conservation targets. Spatial priority conservation areas for water yield adding individual and total costs largely excluded the places where water yield was used for three different water utilizations including electric power production, resident water usage and irrigation for rice production. The spatial priority conservation areas for annual average water yield with the cost calculated from hydropower electricity production was mainly concentrated in the northeastern regions. The number of conservation patches from spatial priority conservation area for water yield was approximately one and the spatial structure was more compact compared to other costs. For the cost calculated from resident water usage, the spatial priority conservation areas for water yield focused on the northern and riverine areas. The common spatial priority conservation areas for water yield with cost from area and cost from resident water usage appeared in the north of the study site. The spatial pattern of priority conservation areas for water yield with cost calculated from irrigation for rice production was similar to that of priority conservation areas for water yield with total cost. This indicated that the cost from irrigation for rice production significantly impacted on total cost and spatial priority conservation of water yield on an annual basis. Moreover, the number of conservation patches was less and the spatial structure was more compact compared to total cost. Some of the priority conservation areas were the same priority conservation areas among three different economic values of water resources. Nonetheless, the selection times (priority ranking) of conservation areas were different under long-term climate change with different cost layers (Figure 15(b)). For the north of the study site, the selection times with cost calculated from resident water usage were greater than those with cost from area. Furthermore, the selection times with the cost calculated from irrigation for rice production were always greater than those with total cost. The spatial priority conservation areas identified for water yield ES with 40% of the maximum total amount target level with different cost layers under long climate change scenario ranged from 38.10 to 39.58% of the landscape (Figure 16(a)). The number of conservation patches increased and turned were more fragmented compared with 20% of the maximum total amount target level. There were more southwestern and southeastern regions involved into spatial priority conservation areas for water yield with area cost. The spatial priority conservation areas for water yield with cost calculated from hydropower electricity production extended from the northeast to the northwest of the study watershed. There were more riverine places involved in the spatial priority conservation areas for water yield with costs from resident water usage, irrigation for rice production and total cost. The larger selection times focused on the southwestern areas with all different cost layers calculated from the three water utilization types. Particularly, the selection times with cost calculated from hydropower electricity production were higher than those with other costs (Figure 16(b)). In a word, spatial priority conservation areas identified for water yield ES with three individual and total costs under the same target levels (20 and 40% of the maximum total amount) were different from those for water yield with area cost. The spatial structures of spatial priority conservation areas for water yield with the cost from hydropower electricity production were more compact than those for water yield with other cost categories. The spatial priority conservation areas for water yield with area cost were closely related to the spatial distribution of water yield ES. The places with higher water yield amounts would be previously involved in a reserve network for water resources.

The pair-wise Jaccard's indexes of spatial priority conservation areas of water yield ES with different cost layers under long-term climate change are shown in Table 5. The Jaccard's indexes of water yield with cost from resident water usage, cost from irrigation for rice production and total cost were high (approximately 0.4). Contrarily, the Jaccard's indexes of water yield with cost from area, cost from hydropower electricity production, cost from irrigation for rice production, and cost from resident water usage were low (approximately 0.04). The Jaccard's index increased as the target level increased from 20 to 40%. The pair-wise correlations of selection times of spatial priority conservation areas of water yield ES with different cost layers under long-term climate change are shown in Table 6. The correlations of selection times of spatial priority conservation areas for water yield with cost from area, cost from hydropower electricity production, cost from resident water usage and total cost were negative. The magnitude of negative correlation became larger as the target level increased from 20 to 40%. It is suggested that there were tradeoffs between spatial priority conservation areas for water yield with cost from area, cost from hydropower electricity production, cost from resident water usage and total cost because of the lower Jaccard's index and negative correlation. The Jaccard's index provides an indication of the amount of flexibility in number and location of water yield with individual and total costs calculated from water utilization, which could measure the ecological, economic and climate effectiveness of protected area networks in water yield ES systematic conservation planning. The correlation of selection times that is widely used as surrogates to measure tradeoff existed in spatial priority conservation areas for water yield with different cost layers from an economic and ecological standpoint in the reserve network designs. The use of Jaccard's index and correlation of selection times could be used to understand the tradeoffs among spatial priority conservation areas for water yield with individual and total costs from the three water utilization types. This will help to promote the construction of highly compact priority areas of watershed conservation systems to satisfy ecological and economic targets.

The high water yield appearing in the southwest of the study watershed was mainly from baseflow, because of its lowest temperature and thickest snowpack. This is typical in snow-dominated cool regions (Marshall & Randhir 2008). There were no crops (which are originally covered by farmland and paddy fields except the winter season, Figure 1(b)) growing in winter which resulted in low ET (Brown et al. 2005). Consequently, this allowed snowmelt water at the bottom of the snowpack to readily penetrate into the soil. Dominant soils were brown lowland and gray lowland at the lower elevation (Figure 1(c)), which easily saturate because of low field capacity water contents and flat topography (Figure 1(d) and 1(e)). The low field capacity water content resulted in low water retention due to low the clay content of the brown lowland and gray lowland (Pachepsky et al. 2001). Moreover, the gentle topography could result in low soil water availability and shorten the lag between when water exits the soil profile and enters the shallow aquifer owing to the short depth to the water table and low water retention (Tani 1997; Soulis & Valiaantzas 2013). Therefore, it suggested that these characteristics, such as land cover (no crops in winter), soil (low field capacity water content) and topography (flat), created the spatial pattern of water yield near the river channel in the southwestern area (Figure 9). The high water yield in the northern region mainly originated from snowmelt, contributing to peak flow during the snowmelt period (March and April). Snowmelt peaked in April, and snowfall contributed almost 50% of the total precipitation. Additionally, the water yield in the west was higher than in other areas, and is associated with pasture and farmland, brown lowland and peat soils, lower temperature, ET, and greater precipitation (including rainfall and snowfall). It also further suggested that greater precipitation in winter (Wassamu, Shibetsu and Bifuka) generated higher water yields from snowmelt in April (Figure 1(a)). Another high value of water yield occurred in the southeastern part because of increased precipitation and its topographical location. This region is also located in the original water source of the Teshio watershed. The slope is steeper in this area than in other places, suggesting that it enhanced the rapid response of water yield to storm rainfall in autumn. The steeper slope may also increase soil water availability on the lower slope, causing rapid water saturation of the soil during storms. That is, the spatial redistribution of surface runoff due to topography resulted in greater soil water availability on the lower slope, which contributed to a spatial water yield gradient on the slope (Tsukamoto & Ohta 1988; Tani 1997; Matson et al. 2002; Lin et al. 2007; Qiu et al. 2017).

The study region is classified as conforming to the ‘temperate zone’, showing a prominent seasonal cycle (Yasunaka & Hanawa 2006; Fan & Shibata 2014). The future climate projections are mostly consistent and the overall increase in mean annual precipitation is the likely consequence of changing climate (Dibike & Coulibaly 2007; IPCC 2007). Meteorological variables such as precipitation and temperature are found to be highly sensitive in predictions of the effect of climate change on water resources. Consequently, investigating the historical and varied trends of such variables would help to reveal the effects of climate change on water resources (Table 1). In this basin, water yield is somewhat more sensitive to changes in precipitation than to changes in temperature, as found by Fan & Shibata (2015). In theory, changes in water yield are reported to be more sensitive to changes in temperature than to changes in precipitation (Singh & Woolhiser 2002; Bouraoui et al. 2004). The simulated results also indicate that the variation in water yield would be larger than the proportional variation in precipitation. Furthermore, the hydrology model had some difficulty in modeling the seasonal water balance components under different climate change scenarios. These problems can be attributed to uncertainty in the input climate data, especially in the spatial and temporal coverage of precipitation and temperature. The weather stations used in this study site were distributed outside the mountainous area, with higher altitudes. Hence, the stations may have significantly different weather conditions from other areas. The mountainous area may enhance the impact of weather on water balance components under different climate change scenarios. Thus, more accurate weather data for mountainous areas may be necessary to obtain more accurate water balance components. Despite some variations in seasonal water balance components, the annual simulation of water yields are quite similar under different climate changes, especially between the observed (validation) period and baseline climate conditions. So, the total water yield results simulated from the hydrology model under different climate change scenarios can be used as input data to following systematic conservation of water yield ESs. Moreover, alternative climate change scenarios should be applied and compared in future related studies.

Water yield (water supply) is a major ES provided by the ecosystems in the watershed of the Teshio River (Ileva et al. 2009; Fan & Shibata 2016). Human preferential utilization for water resources revealed the relationship between the different stakeholders. Comparing the spatial priority conservation areas for water yields with three different cost layers, synergy and tradeoff relationships were observed (Tables 5 and 6). In this study, incorporating cost into hydrological ES prioritization can increase the cost-effectiveness of resulting priority sets and alter the priority attached to conserving different kinds of ecosystems (Marshall & Randhir 2008; Jeppesen et al. 2009; Somura et al. 2009; Chiang et al. 2012; Fan & Shibata 2015; Fan et al. 2017). To demonstrate the potential impact of incorporating cost in water yield ES prioritization exercises, we examined the cost-effectiveness of different economic values derived from water resources utilization for ecoregions across the Teshio watershed. We illustrated changes in cost-effectiveness by considering what could be achieved with diverse conservation targets. Here, the key hydrological provision ES-based management strategy was the most cost-effective strategy with a high conservation value and meeting of conservation targets. The cost-effectiveness means the minimum cost should maximize the conservation of ES (Moilanen et al. 2009; Fan & Shibata 2014; Odgaard et al. 2017). Water yields in some sites which have been already utilized for human activities were assigned higher costs and should be prioritized for conservation. This can guarantee the current provision of hydrological ESs across the watershed. However, unused water yields by humans have lower costs and can be regarded as the potential water resource for future agricultural expansion and urban development. Therefore, the cost efficiency in the systematic conservation model could balance the tradeoff between the protection and water resource development of hydrological provision ESs. We indicated that the selection times (ranking) map for the annual average water yield adding the cost layer differed from the ranking map for annual average yield without adding the cost layer (Figures 15(b) and 16(b)). The area with high ranking in the selection times map of water yield without the cost layer shifted to low rankings in the selection times map of annual average water yield with the cost layer. Similarly, the low ranking areas in the selection times map of water yield turned to high rankings in the selection times map of annual average water yield with the cost layer based on the tradeoff between water yield and its economic value (Naidoo et al. 2006; Arponen et al. 2010).

The three different economic cost indicators are all based on information of biophysical land suitability for water utilization (Fan & Shibata 2014). However, each cost indicator uses different assumptions associated with the size and location of hydro-power plant, resident and paddy field, and relative electricity, tap-water and rice prices. The economic costs vary by several orders of magnitude, such as depending on local economic development, electricity production potential of hydro-power plant, population density and agricultural land distribution. Accordingly, improving cost-effectiveness in watershed conservation prioritization requires spatially explicit and detailed information on conservation costs. The most direct way to consider cost to all groups of stakeholders is to use the sum of the respective costs (Rodríguez et al. 2006). This way, the overall cost can be minimized, but differential effects among different stakeholder groups remain implicit. In contrast, including each cost indicator separately in comparative conservation prioritization analyses makes spatially and socially relevant tradeoffs between stakeholder groups explicit. This approach does not lead to a single solution that includes all costs, but the single-cost scenarios can additionally be compared to a scenario that combines several cost indicators. Besides showing efficiency conservation gains, incorporating economic costs into conservation planning can be used to demonstrate the tradeoffs between obtaining conservation targets and the cost necessary to obtain it (Naidoo et al. 2006; Fan & Shibata 2016).

The construction of spatial conservation areas for hydrological provision ESs is typically fraught with uncertainty and imprecision. However, if people hope that the constructed reserve network can provide a reliable basis for water resource protection and future water utilization and be received by the public, it should be more objective and accurate. In this study, we developed an integrated approach to build spatial priority areas for ESs with and without their economic value effects based upon the hydrology model and systematic conservation model. The combination of the hydrology model (considering relationships between hydrological cycle and vegetation, meteorology, topography, and soil) and systematic conservation model (considering connectivity and cost of ESs reserve network) provides quantitative and spatially explicit conservation areas for water yield ESs across a watershed. In order to link the water yield-based reserve network and conservation targets, we also quantified forest and agricultural land priority areas for water yield with and without cost layers under different maximum total amount target levels. Our proposed approach, based on hydrology cycle processes and economic values of water yield, is a step towards more objective and accurate constructed reserve networks for water resources.

The successful performance of the water yield ESs reserve network was largely due to the innovative ecosystem management systems and mechanisms. In order to improve the actual performance of regional water yield ESs prioritization conservation efforts under climate change scenarios, an adaptive management paradigm needs to be established to integrate the government motivated ‘top-down’ approach (systematic conservation model) and the local stakeholder motivated ‘bottom-up’ approach (economic costs of ESs), with balanced considerations of the dynamics and sustainability requirements of the targeted ecological-socioeconomic coupled systems (Yang & Li 2017). The systematic conservation model is widely used in reversing environmental degradation and can contribute to the improvement of ESs and adaptability to climate variability. An adaptive management approach allows for flexibility in human and financial resource allocation, an expansion of knowledge on the dynamic socioeconomic-ecological coupled systems, and efficacy of management operations. The experience of change in key water yield ESs and their spatial prioritization conservation under climate change scenarios in the Teshio River watershed exemplified the positive effects of environmental policies and the necessity of adopting an adaptive management approach.

In this study, the economic costs derived from the static water supply in the Teshio River watershed did not consider population growth, industrial change and expansion of agricultural land in the future. The dynamic changes in water demands under climate changes are currently not included in spatial prioritization conservation of water resources. An increase in population means greater demand for water, particularly in developing countries, and it is becoming increasingly concentrated in large cities. This would have two implications, on the one hand water use is different in an urban environment than in a rural environment. On the other hand, the increasing population concentration of demand means greater pressure on water resources in specific areas. Industrial development (hydro-power energy production) increases the demand for water, but industrial restructuring may reduce it. As water is seen as more of an indirect economic good, it will be used more efficiently in the industry. The growth in irrigated areas will lead to more usage by agriculture, but this may be offset to a certain extent by improvements in irrigation efficiency. Therefore, future research incorporating feedback of dynamic changes in water demands under climate changes would be necessary to refine spatial prioritization conservation areas for water resources. Such feedback impacts on the spatial distribution of economic costs calculated from water resources which further influences spatial prioritization conservation areas for water yield ESs. In addition, increasing water demands for ecological protection should put additional constraints on water resource use in future research.

This study provided an analytical framework that integrated a catchment-scale hydrology model to a systematic conservation model for constructing spatial optimal priority conservation areas of water yield ES under future climate change scenarios. We analyzed the impact of climate changes on spatial priority conservation areas for water yield ES under 20 and 40% conservation targets in the Teshio River watershed in northernmost Japan. In northernmost Japan, the warming climate in the future will result in significant impacts on water resources as well as the ecosystems dependent on them. Projections to short-, mid-, and long-term effects of future climate scenarios on water resources in northernmost Japan were used to assess the associated risk to the sustainability of water yield ES and water-dependent ecosystems due to changes in water quantity. It is projected that climate change is likely to cause hydrological characteristics and conditions in the watershed, similarly to other regions with a temperate zone climate type (showing prominent seasonal cycle) worldwide. However, the impact is not uniform across the watershed and depends on both the intensity of the changes to water resources and the sensitivity of the water-dependent ecosystems to such changes. The spatial structures of priority conservation areas for water yield with and without cost layers under a 20% conservation target are more compact than those with a 40% target under all climate changes, because the increased precipitation and temperature, spatial characteristics of land use, soil categorization, and topographical driving factor alter the distribution and total amount of water yield ES.

Furthermore, there were interactive relationships (synergies and tradeoffs) among spatial priority conservation areas for water yield with three different cost layers. We investigated how watershed-scale spatial conservation prioritization for water yield is affected by the different indicators of economic values from anthropogenic water utilization. Three spatial conservation prioritizations with three corresponding cost indicators (costs calculated from hydro-power electricity production, resident water-usage, and irrigation for rice production, respectively) were compared, each using one of three economic cost indicators to prioritize planning units according to conservation cost-effectiveness. We found that the different economic cost indicators were at best only weakly to moderately correlated in space when the conservation targets increased from 20 and 40% of total amount water yield hydrological provision ES. The spatial modeling and assessment of water yield ES were integrated into a systematic conservation model and may facilitate a better understanding of the climate change challenges we face and support the formulation of solutions to optimally address those challenges. The systematic conservation of water yield ES has achieved remarkable success in engaging with stakeholders (farmers and environmental managers). The stakeholders can easily understand the conservation debates on spatial priority conservation of water yield ES under climate change. For water yield ES conservation planning in the watershed under future climate change scenarios, it is necessary to couple these models for improved and sustainable management of the watershed ecosystem. The framework used in this project is suitable for watershed-scale estimates to assess climate change impacts on water-dependent ecosystems and water yields, finer temporal-scale modelling and analysis would be required. The reported results provide guidance on where additional water management effort is needed in the study watershed.

This study was partly supported by the National Natural Science Foundation of China (No. 41601088), Key Research and Development Projects of Sichuan Science and Technology Plan (No. 2019YFS0057), Development Program of Science and Technology of Jilin Province (No. 20160520061JH), and a scholarship from the China Scholarship Council and doctor research fund and was supported by Southwest University of Science and Technology (No. 18LZX619; No. 17LZX690). The work was conducted under the Environment Research and Technology Development Fund (S-15) of the Ministry of the Environment, Japan and Integrated Research Program for Advancing Climate Models (TOUGOU program) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. We are also grateful to the editors and reviewers for fundamental improvement of this manuscript.