Abstract

Water resources prioritization conservation planners are increasingly becoming aware of the economic value of water supply ecosystem services (ESs) under climate changes. Here we assessed how the water yield ES framework is implemented in the current spatial prioritization conservation of the water resources under climate change across the Teshio River watershed. We applied the systematic conservation model to optimize the area for water resources which satisfied the protection targets with and without considering economic values of the water yield provision service. The model indicated that the areas of spatial optimal ES protection for water yield with considering economic values were totally different from those without considering economic values of water resources. The optimal priority conservation areas were concentrated in southwestern, southeastern, and some northern areas of this watershed. These places could guarantee water resources sustainability from both environmental protection and socio-economic development standpoints. Moreover, the spatial priority conservation areas for water yield with economic value from hydro-power electricity production were traded off against the areas for water yield with economic values from resident water-use and irrigation for rice. Therefore, the systematic conservation planning of water yield with economic values under climate changes may provide a useful argument to promote the conservation of water resources.

INTRODUCTION

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; (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.

STUDY SITE AND METHODS

Study site

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 m3 s–1, with a maximum of 1,278 m3 s–1 in April from snowmelt and a minimum of 19.8 m3 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)).

Figure 1

Spatial patterns in: (a) Upstream of Teshio River watershed, (b) land use, (c) soil type (RLS (Rocky land), BFS (Brown forest soil), RGS (Regoslos soil), GLL (Gray low land soil), BLL (Brown low land soil), GLS (Gley soil), PDS (Podzol soil), PTS (Peat soil)), (d) field water capacity in soil (mm H2O), (e) slope (%), and (f) precipitation (mm H2O year–1).

Figure 1

Spatial patterns in: (a) Upstream of Teshio River watershed, (b) land use, (c) soil type (RLS (Rocky land), BFS (Brown forest soil), RGS (Regoslos soil), GLL (Gray low land soil), BLL (Brown low land soil), GLS (Gley soil), PDS (Podzol soil), PTS (Peat soil)), (d) field water capacity in soil (mm H2O), (e) slope (%), and (f) precipitation (mm H2O year–1).

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).

Overview of study approach

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.

Figure 2

Flowchart of methods used in this study.

Figure 2

Flowchart of methods used in this study.

Climate change scenarios

Observed climate datasets of precipitation and temperature for 34 years (1976–2009) were used as the baseline climate scenario. We used a GCM prediction model to simulate the future climate during 2010–2039 (hereafter called short-term change), 2040–2069 (mid-term change) and 2070–2099 (long-term change). First, the future climate change scenarios were derived from the Interdisciplinary Research on Climate (MIROC) model (MIROC3.2-HI) under the Special Report on emissions scenarios B1 (SRB1) scenario of Intergovernmental Panel on Climate Change (IPCC 2007). All modeled data for MIROC3.2-HI was obtained from the Data Distribution Centre of IPCC. Since spatial resolutions of GCMs are too coarse to represent local climate characteristics in our watershed, simple downscaling between the baseline and climate scenario of the nearest GCM grid was applied directly. Future change of temperature in the study area was assumed equal to the difference between temperatures simulated using GCMs for future and current conditions at a weather station in the watershed. Thus, future change of temperature was estimated as follows (Tung et al. 2005): 
formula
(1)
Here, and are current and future mean monthly temperatures (°C), respectively; and are simulated mean monthly temperature under the current and future scenario climate conditions, respectively. Future change of precipitation was assumed to be the ratio of precipitation in the future condition to precipitation in the current condition (Tung et al. 2005): 
formula
(2)

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.

Table 1

Changes of temperature and ratios of precipitation predicted by MIROC3.2-HI models based on SRB1 scenarios

Month Station
 
MIROC3.2-HI (44.27°E, 142.88°N)
 
T (°C) P (cm) ΔT (°C)
 
Ratio (cm cm–1)
 
Baseline
 
Short-term change Mid-term change Long-term change Short-term change Mid-term change Long-term change 
–9.20 0.87 1.70 3.30 4.20 1.01 1.03 1.17 
–8.90 0.62 2.00 3.20 4.00 1.11 1.02 1.20 
–3.50 0.61 1.80 3.50 4.20 1.10 1.29 1.24 
3.70 0.50 2.50 4.40 4.90 1.10 1.09 1.16 
10.30 0.59 1.60 2.60 3.30 0.89 0.90 0.98 
15.10 0.60 1.60 2.60 3.30 1.13 1.13 1.23 
19.00 1.08 1.50 3.00 4.20 0.98 1.02 1.27 
20.20 1.26 1.90 2.80 3.50 1.06 1.21 1.18 
14.90 1.35 1.60 2.50 3.30 1.10 1.23 1.39 
10 8.20 1.36 1.80 3.00 3.70 0.93 0.92 0.82 
11 1.20 1.44 1.80 3.00 4.10 1.10 1.05 1.05 
12 –5.30 1.20 1.70 2.70 3.90 1.11 1.17 1.17 
Average 5.50 0.96 1.79 3.05 3.88 1.05 1.09 1.16 
Month Station
 
MIROC3.2-HI (44.27°E, 142.88°N)
 
T (°C) P (cm) ΔT (°C)
 
Ratio (cm cm–1)
 
Baseline
 
Short-term change Mid-term change Long-term change Short-term change Mid-term change Long-term change 
–9.20 0.87 1.70 3.30 4.20 1.01 1.03 1.17 
–8.90 0.62 2.00 3.20 4.00 1.11 1.02 1.20 
–3.50 0.61 1.80 3.50 4.20 1.10 1.29 1.24 
3.70 0.50 2.50 4.40 4.90 1.10 1.09 1.16 
10.30 0.59 1.60 2.60 3.30 0.89 0.90 0.98 
15.10 0.60 1.60 2.60 3.30 1.13 1.13 1.23 
19.00 1.08 1.50 3.00 4.20 0.98 1.02 1.27 
20.20 1.26 1.90 2.80 3.50 1.06 1.21 1.18 
14.90 1.35 1.60 2.50 3.30 1.10 1.23 1.39 
10 8.20 1.36 1.80 3.00 3.70 0.93 0.92 0.82 
11 1.20 1.44 1.80 3.00 4.10 1.10 1.05 1.05 
12 –5.30 1.20 1.70 2.70 3.90 1.11 1.17 1.17 
Average 5.50 0.96 1.79 3.05 3.88 1.05 1.09 1.16 

SWAT model and water yield ecosystem service simulation

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 (R2) 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, Si. 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) Si,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.

Economic value of hydrological provision ecosystem services

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).

Figure 3

Spatial pattern of economic value of: (a) electric power production, (b) resident water-use, (c) irrigation for rice production and (d) summarizing previous three water resources utilization according to water yield in HRU of the SWAT model in the Teshio River watershed (unit: million yen). See details of water yield hydrological ecosystem service economic valuation procedure in Fan & Shibata's (2014) related study.

Figure 3

Spatial pattern of economic value of: (a) electric power production, (b) resident water-use, (c) irrigation for rice production and (d) summarizing previous three water resources utilization according to water yield in HRU of the SWAT model in the Teshio River watershed (unit: million yen). See details of water yield hydrological ecosystem service economic valuation procedure in Fan & Shibata's (2014) related study.

The total economic values of hydrological provision ESs included three parts, and detailed information is shown in Equation (3): 
formula
(3)
The calculation of economic values of hydropower electricity production, irrigation for crop production, and resident water use are as follows:
  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): 
    formula
    (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): 
    formula
    (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: 
    formula
    (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.

MARXAN model

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.

The objective function in MARXAN is designed so that the lower the value the better, MARXAN's objective function calculates the total cost associated with a set of planning units as: 
formula
(7)
where PUs refers to planning units, SPF to conservation feature penalty factor, and BLM to boundary length modifier. The SPF weighs the cost associated with failure to meet the representation target of each conservation feature specified, this feature is hydrological ESs in the present analysis. The penalty component of the MARXAN objective function is the penalty to a reserve system for not adequately representing conservation features. The threshold penalty represents a cost for exceeding some maximum desired size of solutions. Especially, the parameters BLM and SPF play an important role in spatial configuration of systematic conservation for ESs. The BLM is a weight used to control the spatial aggregation of planning units, and boundary is the length boundary surrounding the reserve system. By including a boundary length term in the objective function we can control the level of fragmentation in the reserve system (BLM = 1 in this study). An SPF determines the importance of meeting targets, so higher penalties can be set for not meeting targets for the most important features to increase the likelihood of the target being met, or penalties can be set high for all features which are a requirement for conservation targets (SPF values for water yield hydrological ESs in this study is 4. This study used the MARXAN model (Ball & Possingham 2001) to perform simulated annealing, and a Monte Carlo procedure was used to minimize multivariate functions associated with a physical or biological system (100 runs with 1,000,000 iterations).

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.

Fitting regression models of agricultural and forest land spatial priority conservation areas for water yield hydrological ESs and conservation target levels

So as to quantify the varied trends of agricultural and forest land spatial priority conservation areas for annual average water yield as changes of conservation targets, the fitting regression analysis was used to examine their relationships based on their marked curves (Figures 1014). All the programming codes for fitting regression analysis were implemented in WinBUGS (Gelman & Hill 2006; Fan & Shibata 2016). The fitting regression model of spatial priority conservation areas for annual average water yield and conservation target levels can be written as: 
formula
(8)
where i is the number of sample (, the number of n is equivalent to 9 in this study). is the conservation target level of annual average water yield including 10–90% of the maximum total amount values of water yield hydrological ES, and is total spatial priority conservation areas for annual average water yield corresponding to the conservation target level of annual average water yield.

Comparison of spatial patterns of conservation areas for water yield hydrological ESs

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).

RESULTS

Simulated results of water yield hydrological provision ecosystem services under climate change scenarios

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, R2 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.

Table 2

The most sensitive parameters and final calibrated values of parameters

Parameter Description Sensitivity index Minimum Maximum Optimal value 
CN2 Initial SCS runoff curve number 0.146 35 98 73 
ALPH_BF Baseflow alpha factor (days) 0.075 0.05 
SFTMP Snowfall temperature (°C) 0.072 −5 
SMTMP Snowmelt base temperature (°C) 0.065 −5 
SMFMX Maximum melt rate for snow during years (mm °C−1 d−10.061 10 
SMFMN Minimum melt rate for snow during years (mm °C−1 d−10.027 10 0.1 
TIMP Snowpack temperature lag factor 0.025 0.9 
ESCO Soil evaporaion compensation factor 0.019 0.75 
SURLAG Surface runoff lag coefficient 0.014 0.05 24 
GW_DELAG Groundwater delay (days) 0.013 500 20 
Parameter Description Sensitivity index Minimum Maximum Optimal value 
CN2 Initial SCS runoff curve number 0.146 35 98 73 
ALPH_BF Baseflow alpha factor (days) 0.075 0.05 
SFTMP Snowfall temperature (°C) 0.072 −5 
SMTMP Snowmelt base temperature (°C) 0.065 −5 
SMFMX Maximum melt rate for snow during years (mm °C−1 d−10.061 10 
SMFMN Minimum melt rate for snow during years (mm °C−1 d−10.027 10 0.1 
TIMP Snowpack temperature lag factor 0.025 0.9 
ESCO Soil evaporaion compensation factor 0.019 0.75 
SURLAG Surface runoff lag coefficient 0.014 0.05 24 
GW_DELAG Groundwater delay (days) 0.013 500 20 
Figure 4

Spatial distribution of subwatersheds in Teshio River catchment obtained from the SWAT model.

Figure 4

Spatial distribution of subwatersheds in Teshio River catchment obtained from the SWAT model.

Figure 5

Spatial distribution of hydrologic response units (HRUs) in Teshio River catchment obtained from the SWAT model (including four local enlarged drawings).

Figure 5

Spatial distribution of hydrologic response units (HRUs) in Teshio River catchment obtained from the SWAT model (including four local enlarged drawings).

Figure 6

Temporal fluctuation in observed (blue diamond with solid line) and simulated (red square with dash line) of stream flow from 2003 to 2009 at five monitoring stations (Oku-Shibetsu, Shibetsu, Makunbetsu, Nayoro oohashi and Bifuka bashi) in the Teshio river watershed. Calibration and validation of streamflow in the Teshio watershed, coefficient of determination (R2) and Nash–Sutcliffle efficiency (NES) were greater than 0.85. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2019.009.

Figure 6

Temporal fluctuation in observed (blue diamond with solid line) and simulated (red square with dash line) of stream flow from 2003 to 2009 at five monitoring stations (Oku-Shibetsu, Shibetsu, Makunbetsu, Nayoro oohashi and Bifuka bashi) in the Teshio river watershed. Calibration and validation of streamflow in the Teshio watershed, coefficient of determination (R2) and Nash–Sutcliffle efficiency (NES) were greater than 0.85. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2019.009.

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 five months (May–September) under baseline climate conditions to seven 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 H2O. 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.

Figure 7

Annual average precipitation, snowfall, water yield and evapotranspiration under climate changes. B, S, M and L in the scenarios mean climate changes in baseline (observed for 1976–2009), short-term prediction (2010–2039), mid-term prediction (2040–2069) and long-term prediction (2070–2099), respectively.

Figure 7

Annual average precipitation, snowfall, water yield and evapotranspiration under climate changes. B, S, M and L in the scenarios mean climate changes in baseline (observed for 1976–2009), short-term prediction (2010–2039), mid-term prediction (2040–2069) and long-term prediction (2070–2099), respectively.

Figure 8

Water-balance components under climate changes.

Figure 8

Water-balance components under climate changes.

Figure 9

Spatial pattern of water yield in HRU of the SWAT model under long-term climate change scenario across the Teshio River watershed (unit: mm H2O).

Figure 9

Spatial pattern of water yield in HRU of the SWAT model under long-term climate change scenario across the Teshio River watershed (unit: mm H2O).

Figure 10

(a) Agricultural land priority areas for annual average water yield without cost layer under climate change scenarios; (b) Forest land priority areas for annual average water yield without cost layer under climate change scenarios.

Figure 10

(a) Agricultural land priority areas for annual average water yield without cost layer under climate change scenarios; (b) Forest land priority areas for annual average water yield without cost layer under climate change scenarios.

Figure 11

(a) Agricultural land priority areas for annual average water yield with cost layer calculated from electric power production under climate change scenarios; (b) Forest land priority areas for annual average water yield with cost layer calculated from electric power production under climate change scenarios.

Figure 11

(a) Agricultural land priority areas for annual average water yield with cost layer calculated from electric power production under climate change scenarios; (b) Forest land priority areas for annual average water yield with cost layer calculated from electric power production under climate change scenarios.

Figure 12

(a) Agricultural land priority areas for annual average water yield with cost layer calculated from resident water-use under climate change scenarios; (b) Forest land priority areas for annual average water yield with cost layer calculated from resident water-use under climate change scenarios.

Figure 12

(a) Agricultural land priority areas for annual average water yield with cost layer calculated from resident water-use under climate change scenarios; (b) Forest land priority areas for annual average water yield with cost layer calculated from resident water-use under climate change scenarios.

Figure 13

(a) Agricultural land priority areas for annual average water yield with cost layer calculated from irrigation for rice production under climate change scenarios; (b) Forest land priority areas for annual average water yield with cost layer calculated from irrigation for rice production under climate change scenarios.

Figure 13

(a) Agricultural land priority areas for annual average water yield with cost layer calculated from irrigation for rice production under climate change scenarios; (b) Forest land priority areas for annual average water yield with cost layer calculated from irrigation for rice production under climate change scenarios.

Figure 14

(a) Agricultural land priority areas for annual average water yield with cost layer calculated from summing electric power production, resident water-use and irrigation for rice production under climate change scenarios; (b) Forest land priority areas for annual average water yield with cost layer calculated from summing electric power production, resident water-use and irrigation for rice production under climate change scenarios.

Figure 14

(a) Agricultural land priority areas for annual average water yield with cost layer calculated from summing electric power production, resident water-use and irrigation for rice production under climate change scenarios; (b) Forest land priority areas for annual average water yield with cost layer calculated from summing electric power production, resident water-use and irrigation for rice production under climate change scenarios.

Agricultural land and forest land spatial conservation prioritization areas for water yield ecosystem service under long-term climate change

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 1014). 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 1014; 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.

Table 3

Fitting regression equations of agricultural land priority areas for annual average water yield and conservation targets with different cost layers under climate change scenarios

Cost layers Climate change scenarios
 
Baseline climate conditions Short-term climate change Mid-term climate change Long-term climate change 
Cost from area y = 0.02323 + 0.9523x y = 0.0351 + 0.9268x y = 0.01828 + 0.965x y = 0.01887 + 0.9607x 
Cost from electric power production y = 0.02097 + 0.958x y = 0.02089 + 0.9579x y = 0.01748 + 0.9654x y = 0.01603 + 0.9684x 
Cost from resident water-use y = 0.05572 + 0.8892x y = 0.05459 + 0.8929x y = 0.04834 + 0.9028x y = 0.04588 + 0.9091x 
Cost from irrigation for rice production y = 0.05085 + 0.8987x y = 0.05303 + 0.8932x y = 0.05628 + 0.89x y = 0.05428 + 0.8891x 
Total cost y = 0.04256 + 0.9162x y = 0.03851 + 0.9239x y = 0.04228 + 0.914x y = 0.04325 + 0.9114x 
Cost layers Climate change scenarios
 
Baseline climate conditions Short-term climate change Mid-term climate change Long-term climate change 
Cost from area y = 0.02323 + 0.9523x y = 0.0351 + 0.9268x y = 0.01828 + 0.965x y = 0.01887 + 0.9607x 
Cost from electric power production y = 0.02097 + 0.958x y = 0.02089 + 0.9579x y = 0.01748 + 0.9654x y = 0.01603 + 0.9684x 
Cost from resident water-use y = 0.05572 + 0.8892x y = 0.05459 + 0.8929x y = 0.04834 + 0.9028x y = 0.04588 + 0.9091x 
Cost from irrigation for rice production y = 0.05085 + 0.8987x y = 0.05303 + 0.8932x y = 0.05628 + 0.89x y = 0.05428 + 0.8891x 
Total cost y = 0.04256 + 0.9162x y = 0.03851 + 0.9239x y = 0.04228 + 0.914x y = 0.04325 + 0.9114x 
Table 4

Fitting regression equations of forest land priority areas for annual average water yield and conservation targets with different cost layers under climate change scenarios

Cost layers Climate change scenarios
 
Baseline climate conditions Short-term climate change Mid-term climate change Long-term climate change 
Cost from area y = 6.977E–4 + 0.9987x y = 4.807E–4 + 0.9991x y = 4.24E–4 + 0.9993x y = 5.667E–4 + 0.999x 
Cost from electric power production y = 0.002092 + 0.9962x y = 0.001509 + 0.997x y = 0.002028 + 0.996x y = 0.001142 + 0.9981x 
Cost from resident water-use y = 0.003688 + 0.9933x y = 0.004386 + 0.991x y = 0.005691 + 0.9884x y = 0.007617 + 0.9852x 
Cost from irrigation for rice production y = 9.721E–4 + 0.9976x y = 0.002334 + 0.9954x y = 0.002037 + 0.9958x y = 0.002113 + 0.9958x 
Total cost y = 0.002286 + 0.9957x y = 0.004601 + 0.9916x y = 0.001845 + 0.9968x y = 0.001486 + 0.9969x 
Cost layers Climate change scenarios
 
Baseline climate conditions Short-term climate change Mid-term climate change Long-term climate change 
Cost from area y = 6.977E–4 + 0.9987x y = 4.807E–4 + 0.9991x y = 4.24E–4 + 0.9993x y = 5.667E–4 + 0.999x 
Cost from electric power production y = 0.002092 + 0.9962x y = 0.001509 + 0.997x y = 0.002028 + 0.996x y = 0.001142 + 0.9981x 
Cost from resident water-use y = 0.003688 + 0.9933x y = 0.004386 + 0.991x y = 0.005691 + 0.9884x y = 0.007617 + 0.9852x 
Cost from irrigation for rice production y = 9.721E–4 + 0.9976x y = 0.002334 + 0.9954x y = 0.002037 + 0.9958x y = 0.002113 + 0.9958x 
Total cost y = 0.002286 + 0.9957x y = 0.004601 + 0.9916x y = 0.001845 + 0.9968x y = 0.001486 + 0.9969x 

Spatial conservation prioritization areas for water yield ecosystem service under long-term climate change

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.

Figure 15

(a) Spatial priority areas for water yield ecosystem service; (b) Selection times (ranking) with 20% of the maximum total amount target level under long-term climate change.

Figure 15

(a) Spatial priority areas for water yield ecosystem service; (b) Selection times (ranking) with 20% of the maximum total amount target level under long-term climate change.

Figure 16

(a) Spatial priority areas for water yield ecosystem service; (b) Selection times (ranking) with 40% of the maximum total amount target level under long-term climate change.

Figure 16

(a) Spatial priority areas for water yield ecosystem service; (b) Selection times (ranking) with 40% of the maximum total amount target level under long-term climate change.

Spatial correlated relationships between conservation prioritization areas for water yield ecosystem service with different cost layers under long-term climate change

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.

Table 5

Pair-wise Jaccard's indexes of spatial priority conservation areas of water yield ecosystem service with different cost layers under long-term climate change

Scenario Cost from area Cost from electric power production Cost from resident water-use Cost from irrigation for rice production Total cost 
Cost from area 1.000 0.128 0.153 0.0492 0.0482 
Cost from electric power production 0.176 1.000 0.0417 0.0584 0.124 
Cost from resident water-use 0.171 0.181 1.000 0.379 0.383 
Cost from irrigation for rice production 0.187 0.272 0.598 1.000 0.622 
Total cost 0.189 0.366 0.577 0.745 1.000 
Scenario Cost from area Cost from electric power production Cost from resident water-use Cost from irrigation for rice production Total cost 
Cost from area 1.000 0.128 0.153 0.0492 0.0482 
Cost from electric power production 0.176 1.000 0.0417 0.0584 0.124 
Cost from resident water-use 0.171 0.181 1.000 0.379 0.383 
Cost from irrigation for rice production 0.187 0.272 0.598 1.000 0.622 
Total cost 0.189 0.366 0.577 0.745 1.000 

Notes: Pair-wise Jaccard's index of optimal priority conservation areas of water yields with 20% of the maximum total amount target level under long-term climate change and pair-wise Jaccard's index of optimal priority conservation areas of water yields with 40% of the maximum total amount target level under long-term climate change.

Table 6

Pair-wise correlations of selection times of spatial priority conservation areas of water yield ecosystem service with different cost layers under long-term climate change

Scenario Cost from area Cost from electric power production Cost from resident water-use Cost from irrigation for rice production Total cost 
Cost from area 1.000 –0.0466 0.0235 –0.0874 –0.143 
Cost from electric power production –0.111 1.000 –0.171 –0.104 0.0399 
Cost from resident water-use –0.00900 –0.154 1.000 0.518 0.542 
Cost from irrigation for rice production –0.0282 0.0618 0.718 1.000 0.889 
Total cost –0.0253 0.205 0.687 0.941 1.000 
Scenario Cost from area Cost from electric power production Cost from resident water-use Cost from irrigation for rice production Total cost 
Cost from area 1.000 –0.0466 0.0235 –0.0874 –0.143 
Cost from electric power production –0.111 1.000 –0.171 –0.104 0.0399 
Cost from resident water-use –0.00900 –0.154 1.000 0.518 0.542 
Cost from irrigation for rice production –0.0282 0.0618 0.718 1.000 0.889 
Total cost –0.0253 0.205 0.687 0.941 1.000 

Notes: Pair-wise Jaccard's index of optimal priority conservation areas of water yields with 20% of the maximum total amount target level under long-term climate change and pair-wise Jaccard's index of optimal priority conservation areas of water yields with 40% of the maximum total amount target level under long-term climate change.

DISCUSSION

Impacts of climate changes on water yield ecosystem services

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.

Impacts of economic values of water resources on spatial prioritization conservation areas for water yield ecosystem services

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. 2009a; 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).

Implications for spatial prioritization conservation areas for water yield ecosystem services and future work

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.

CONCLUSIONS

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.

ACKNOWLEDGEMENTS

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.

REFERENCES

REFERENCES
Abbaspour
K. C.
,
Yang
J.
,
Maximov
I.
,
Siber
R.
,
Mieleitner
J.
,
Zobrist
J.
&
Srinivasan
R.
2007
Modelling hydrology and water quality in the pre-alpine/alpine Thur watershed using SWAT
.
J. Hydrol.
333
(
2–4
),
413
430
.
Airame
S.
,
Dugan
J. E.
,
Lafferty
K. D.
,
Leslie
H.
,
McArdle
D. A.
&
Warner
R. R.
2003
Applying ecological criteria to marine reserve design: a case study from the California Channel Islands
.
Ecol. Appl.
13
,
S170
S184
.
Arnold
J. G.
&
Allen
P. M.
1996
Estimating hydrologic budgets for three Illinois watersheds
.
J. Hydrol.
176
,
57
77
.
Arponen
A.
,
Cabeza
M.
,
Eklund
J.
,
Kujala
H.
&
Lehtomäki
J.
2010
Costs of integrating economics and conservation planning
.
Conserv. Biol.
24
(
5
),
1198
1204
.
Ball
I. R.
&
Possingham
H. P.
2001
MARXAN – A Reserve System Selection Tool
.
The Ecology Center, University of Queensland
,
Australia
.
Barron
O.
,
Silberstein
R.
,
Ali
R.
,
Donohue
R.
,
McFarlane
D. J.
,
Davies
P.
,
Hodgson
G.
,
Smart
N.
&
Donn
M.
2012
Climate change effects on water-dependent ecosystems in south-western Australia
.
J. Hydrol.
434–435
,
95
109
.
Bouraoui
F.
,
Grizzetti
B.
,
Grandlund
K.
,
Rekolainen
S.
&
Bidodlio
G.
2004
Impact of climate change on the water cycle and nutrient losses in a Finnish catchment
.
Clim. Change
66
,
109
126
.
Brown
A. E.
,
Zhang
L.
,
McMahon
T. A.
,
Western
A. W.
&
Vertessy
R. A.
2005
A review of paired catchment studies for determining changes in water yield resulting from alterations in vegetation
.
J. Hydrol.
310
,
28
61
.
Chan
K. M. A.
,
Shaw
M. R.
,
Cameron
D. R.
,
Underwoood
E. C.
&
Daily
G. C.
2006
Conservation planning for ecosystem services
.
PLoS Biol.
4
(
11
),
e379
,
2138–2152
.
Chiang
L. C.
,
Chaubey
I.
,
Hong
N. M.
,
Lin
Y. P.
&
Huang
T.
2012
Implementation of BMP strategies for adaptation to climate change and land use change in a pasture-dominated watershed
.
Int. J. Environ. Res. Public Health
9
,
3654
3684
.
Costanza
R.
,
d'Arge
R.
,
de Groot
R.
,
Farber
S.
,
Grasso
M.
,
Hannon
B.
,
Limburg
K.
,
Naeem
S.
,
O'Neill
R. V.
,
Paruelo
J.
,
Raskin
R. G.
,
Sutton
P.
&
van den Belt
M.
1997
The value of the world's ecosystem services and natural capital
.
Nature
387
,
253
259
.
Costanza
R.
,
de Groot
R.
,
Braat
L.
,
Kubiszewski
I.
,
Fioramonti
L.
,
Sutton
P.
,
Farber
S.
&
Grasso
M.
2017
Twenty years of ecosystem services: how far have we come and how far do we still need to go?
Ecosyst. Serv.
28
,
1
16
.
Culver
D. C.
,
Christman
M. C.
,
Sket
B.
&
Trontelj
P.
2004
Sampling adequacy in an extreme environment: species richness patterns in Slovenian caves
.
Biol. Conserv.
112
,
191
126
.
Egoh
B. N.
,
Reyers
B.
,
Carwardine
J.
,
Bode
M.
,
O'Farrell
P.
,
Wilson
K. A.
,
Possingham
H. P.
,
Rouget
M.
,
De Lange
W.
,
Richardson
D. M.
&
Cowling
R. M.
2010
Safeguarding biodiversity and ecosystem services in the Little Karoo, South Africa
.
Conserv. Biol.
24
(
4
),
1021
1030
.
Egoh
B.
,
Reyers
B.
,
Rouget
M.
&
Richardson
D. M.
2011
Identifying priority areas for ecosystem service management in South African grasslands
.
J. Environ. Manage.
92
,
1642
1650
.
Gelman
A.
&
Hill
J.
2006
Data Analysis Using Regression and Multilevel/Hierarchical Models
.
Cambridge University Press
,
Cambridge
,
UK
.
Gosselink
J. G.
,
Odum
E. P.
&
Pope
R. M.
1974
The Value of the Tidal Marsh
.
Center for Wetlands Resources, Louisiana State University
,
Baton Rouge, LA
.
IPCC
2007
Climate change 2007: Impacts, Adaptation and Vulnerability
. In:
Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change
(
Parry
M. L.
,
Canziani
O. F.
,
Palutikof
J. P.
,
van der Linden
P.
&
Hanson
C.
, eds).
Cambridge University Press
,
Cambridge
,
UK
, pp.
1
27
.
IUSS Working Group WRB
2006
World Reference Base for Soil Resources 2006
.
World Soil Resources Reports No. 103. FAO
,
Rome
.
Jeppesen
E.
,
Kronvang
B.
,
Meerhoff
M.
,
Søndergaard
M.
,
Hansen
K. M.
,
Andersen
H. E.
,
Lauridsen
T. L.
,
Liboriussen
L.
,
Beklioglu
M.
,
Ozen
A.
&
Olesen
J. E.
2009
Climate change effects on runoff, catchment phosphorous loading and lake ecological state, and potential adaptations
.
J. Environ. Qual.
38
,
1930
1941
.
Katsuyama
M.
,
Shibata
H.
,
Yoshioka
T.
,
Yoshida
T.
,
Ogawa
A.
&
Ohte
N.
2009
Applications of a hydro-biogeochemical model and long-term simulations of the effects of logging in forested watersheds
.
Sustain. Sci.
4
,
189
198
.
Leathwick
J. R.
,
Moilanen
A.
,
Francis
M.
,
Elith
J.
,
Taylor
P.
,
Julian
K.
,
Hastie
T.
&
Duffy
C.
2008
Novel methods for design and evaluation of marine protected areas in offshore waters
.
Conserv. Lett.
1
,
91
102
.
Moilanen
A.
,
Leathwick
J.
&
Elith
J.
2008
A method for spatial freshwater conservation prioritization
.
Freshw. Biol.
53
,
577
592
.
Moilanen
A.
,
Arponen
A.
,
Stockland
J. N.
&
Cabeza
M.
2009a
Assessing replacement cost of conservation areas: how does habitat loss influence priorities?
Biol. Conserv.
142
,
575
585
.
Naidoo
R.
,
Balmford
A.
,
Ferraro
P. J.
,
Polasky
S.
,
Ricketts
T. H.
&
Rouget
M.
2006
Integrating economic costs into conservation planning
.
Trends Ecol. Evol.
21
(
12
),
681
687
.
Neitsch
S. L.
,
Arnold
J. G.
&
Krniry
J. R.
2011
Soil and Water Assessment Tool Theoretical Documentation Version 2009
.
Texas Water Resources Institute
,
College Station, Texas
.
Odgaard
M. V.
,
Turner
K. G.
,
Bøcher
P. K.
,
Svenning
J. C.
&
Dalgaard
T.
2017
A multi-criteria, ecosystem-service value method used to assess catchment suitability for potential wetland reconstruction in Denmark
.
Ecol. Indic.
77
,
151
165
.
Pachepsky
Y. A.
,
Timlin
D. J.
&
Rawls
W. J.
2001
Soil water retention as related to topographic variables
.
Soil Sci. Am. J.
65
,
1787
1795
.
Pan
X. L.
,
Xu
L. Y.
,
Yang
Z. F.
&
Yu
B.
2017
Payments for ecosystem services in China: policy, practice, and progress
.
J. Clean. Prod.
158
,
200
208
.
Possingham
H. P.
,
Ball
I. R.
,
Andelman
S.
2000
Mathematical methods for identifying representative reserve networks
. In:
Quantitative Methods for Conservation Biology
(
Ferson
S.
&
Burgman
M.
, eds).
Springer-Verlag
,
New York
, pp.
3
30
.
Qiu
L. J.
,
Wu
Y. P.
,
Wang
L. J.
,
Lei
X. H.
,
Liao
W. H.
,
Hui
Y.
&
Meng
X. Y.
2017
Spatiotemporal response of the water cycle to land use conversions in a typical hilly-gully basin on the Loess Plateau, China
.
Hydrol. Earth Syst. Sci.
21
,
6485
6499
.
Rodríguez
J. P.
,
Bread
T. D.
,
Bennett
E. M.
&
Cumming
G. S.
2006
Trade-offs across space, time, and ecosystem services
.
Ecol. Soc.
11
(
1
),
1
15
.
Schuol
J.
,
Abbaspour
K. C.
,
Yang
H.
,
Srinivasan
R.
&
Zehnder
A. J. B.
2008
Modeling blue and green water availability in Africa
.
Water Resour. Res.
44
(
7
),
W07406
,
1-18
.
Sharpley
A. N.
&
Williams
J. R.
1990
EPIC-Erosion Productivity Impact Calculator, 1. Model Documentation
.
Tech Bull, 1768. Agricultural Research Service
,
U.S. Department of Agriculture
,
Washington, DC
,
USA
.
Shriner
S. A.
,
Wilson
K. R.
&
Flather
C. H.
2006
Reserve networks based on richness hotspots and representation vary with scale
.
Ecol. Appl.
16
,
1660
1673
.
Singh
V. P.
&
Woolhiser
D. A.
2002
Mathematical modeling of watershed hydrology
.
J. Hydrol. Eng.
7
,
270
292
.
Tsukamoto
Y.
&
Ohta
T.
1988
Runoff process on a steep forested slope
.
J. Hydrol.
102
,
165
178
.
Wallace
K. J.
2007
Classification of ecosystem services: problems and solutions
.
Biol. Conserv.
139
(
3–4
),
235
246
.
Warman
L. D.
,
Sinclair
A. R. E.
,
Scudder
G. G. E.
,
Klinkenberg
B.
&
Pressey
R. L.
2004
Sensitivity of systematic reserve selection to decisions about scale, biological data, and targets: case study from Southern British Columbia
.
Conserv. Biol.
18
,
655
666
.
Yang
W. X.
&
Li
L. G.
2017
Analysis of total factor efficiency of water resource and energy in China: a study based on DEA-SBM model
.
Sustainability
9
,
1
21
.