Spatial priority conservation areas for water yield ecosystem service under climate changes in Teshio watershed, northernmost Japan

Water yield provision as one ecosystem service plays an important role in economic and agricultural development, including water supply for irrigation and drinking water, energy production, flood control and so on (Millennium Ecosystem Assessment 2005). According to the Intergovernmental Panel on Climate Change (IPCC), climate change can influence key temporal and spatial patterns and fluxes of water yield in a watershed, particularly in regions where snowmelt is the primary component of the water balance budget (IPCC 2007). It also further affects the potential availability of freshwater for both ecosystems and humans in a river basin (e.g. the sustainability of community water supplies and hydroelectric power potential) (Somura et al. 2009; Fan & Shibata 2014). Effects of changes in the energy balance of the climate system have been observed in Japan. Especially, in Hokkaido, an island in the northern part of Japan, the increase of temperature is predicted to reach nearly 4 °C, where the study area (Teshio River watershed) is located (JMA 2009). The climate of this area is classified as a temperate zone, showing a prominent seasonal cycle (Yasunaka & Hanawa 2006).

It has been recognized that climate change would affect the hydrological cycle and various processes of a watershed system in many previous studies including changes in total water yield, nutrient and sediment productions (Jha et al. 2004; Ding et al. 2016; Jiang et al. 2016). In particular, some studies have suggested that an increase or decrease in winter water yield and spring water yield can be caused by seasonal shifts in the snow accumulation pattern because the water yield of winter precipitation as opposed to the formation of snow pack is affected by climate change (Bouraoui et al. 2004; Brown et al. 2005; Marshall & Randhir 2008; Fan & Shibata 2016).

The potential influence of climate change is therefore a major concern for watershed management and policy. In the face of these changes, it is important to integrate and extend systematic conservation planning for monitoring and reporting water yield ecosystem services to mitigate the influences of climate change. Systematic conservation planning focuses on locating, designing and managing the protected areas that comprehensively represent water yields at a range of spatial and temporal scales under natural and anthropogenic disturbances (Boteva et al. 2004; Chan et al. 2006; Nelson et al. 2008). Systematic conservation planning is a framework for conservation prioritization and large-scale spatial conservation planning. It identifies areas of landscapes that are important for retaining ecosystem services and connectivity simultaneously for multiple ecosystem services, thus which can provide enhanced balance between ecosystem services conservation and economic development in the long term (Moilanen 2007, 2008; Moilanen & Kujala 2008). Furthermore, systematic conservation planning has been applied in large-scale terrestrial data (Moilanen et al. 2005; Nelson et al. 2009; Fan & Shibata 2016), high resolution planning at the scale of an entire biodiversity hotspot (Chan et al. 2006; Kremen et al. 2008; Egoh et al. 2011), identification of marine protected areas (Leathwick et al. 2008), and in forestry (Mikusinski et al. 2007) and urban planning (Gordon et al. 2009). In the systematic conservation planning of ecosystem services, the systematic conservation model can take into account distribution of ecosystem services, connectivity responses, variable emphasis placed on highest local quality of ecosystem services and cost efficiency analysis for ecosystem services planning (Naidoo et al. 2008; Nelson et al. 2008).

The principal challenges in managing water yield ecosystem services across different periods are that they are not independent of each other, and that the relationships between them may be highly non-linear (Tilman et al. 2002; Rodríguez et al. 2006). Because the hydrological flow and cycles largely fluctuate with various spatial and time scales, it is often difficult to depict their characteristics in local and regional scales with heterogeneous environments. Besides, water yield always changes in different water periods mainly influenced by meteorological factors (e.g. precipitation and temperature); the heterogeneity of water yield is also attributed to topography, land use and soil property.

Knowledge and awareness of interactive relationships among water yields in typical periods under climate change are necessary for making sound decisions about how to manage water resources appropriately. However, few studies have defined temporal water yield ecosystem service planning processes under climate change by coupling a hydrology model with a systematic conservation model (e.g. Guo et al. 2000; Naidoo & Adamowicz 2005; Nelson et al. 2009; Rodriguez et al. 2011; Tortajada & Joshi 2013; Fan & Shibata 2016). This study focuses on the water yields of typical periods as key features of ecosystem services to construct their spatial priority conservation areas under climate change from the temporal and spatial scales, which allows stakeholders to understand the conservation criteria. Moreover, spatially optimal priority conservation from the perspective of water yields could increase the potential for securing funding implementation including economic values of hydropower electricity production, irrigation for crop production and resident water use (as payoff fees) (Pagiola et al. 2010; Fan & Shibata 2014). Therefore, this study aims at evaluating the role of climate change on water yields at temporal and spatial scales and discusses mitigation strategies by constructing optimal priority conservation areas for water yields.

The objectives of the study are the following: (1) to estimate changes in water yields across typical periods caused by climate change and depict their spatial distribution patterns; and (2) to identify spatial priority conservation areas where conservation efforts should be directed for water yields under climate change at watershed scale and to quantify the interactive relationships among water yields in these typical periods.

The overall analytical framework includes modeling water yields across different periods (water yields in February, April and October, composite and annual average water yields) under climate change in the Teshio River watershed, and simulating optimal priority conservation areas for water yields in typical periods. We simulated the temporal and spatial changes in water yields using the hydrology model (Soil and Water Assessment Tool, SWAT) (Arnold & Allen 1996; Abbaspour et al. 2007). We then simulated the optimal priority conservation areas for water yields using the systematic conservation model (MARXAN model) (Airame et al. 2003; Cook & Auster 2005). The flow-chart of the methodology adopted in this study is shown in Figure 1.

The study site is in the Teshio River catchment of northern Hokkaido. This river is the fourth longest (256 km) in Japan, originating from Mount Teshio and emptying into the Sea of Japan. The agricultural land, population and three hydropower plants are concentrated in the upper and middle watershed, implying that the relationships between water yield and anthropogenic activity are more active in the upper and middle watershed than in the lower watershed. Therefore, this study focused on the upper and middle watershed which is located at 44.33° north, 142.25° east (Figure 2).

The catchment area of the study site is 2,908 km². Average water flow is 138.6 m³ s^–1, with maximum 1,278 m³ s^–1 in April from snowmelt, minimum 19.8 m³ s^–1 in February during mid-snowpack season, and a second smaller peak at 597.4 m³ s^–1 in October (between late summer and early autumn) at Bifuka monitoring station from 2001 to 2009 (Ileva et al. 2009; Katsuyama et al. 2009). The standard deviation and skewness of monthly average water flow during this period are 42.62 and 1.03, respectively. Approximately 78% of the catchment is covered by forest categorized as cool-temperate mixed forest, including deciduous broadleaf and evergreen coniferous species with the dense understory of Sasa dwarf bamboo. Other land uses are mainly farmland and paddy fields, with area percentages 13% and 4%, respectively (Figure 3(a)). The farmland and paddy fields are concentrated on both sides of the main Teshio channel and in the southwestern watershed where the field water capacity in soil is low (Figure 3(b)). The slope ranges from 0 to 83.8% with an average of 14.5% (Figure 3(c)). The soil is dominated by brown forest soil (Cambisol; IUSS Working Group WRB 2006) which is located in steep slope regions; others are gray lowland soil (Gleyic Fluvisols; IUSS Working Group WRB 2006), brown lowland soil (Haplic Fluvisols; IUSS Working Group WRB 2006), gray soil (Gleyic Fluvisols; IUSS Working Group WRB 2006), and peat soil (Histosols; IUSS Working Group WRB 2006).

We used the climate change scenario originated from Fan & Shibata (2015). Observed climate data of precipitation and temperature for 34 years (1976–2009) were used as the baseline climate conditions. We used a general circulation model (GCM) prediction model to simulate future climate 2010–2039 (hereafter called short-term change), 2040–2069 (mid-term change) and 2070–2099 (long-term change) (Figure 4 and Table 1).

The future climate change scenario was derived from the Model for Interdisciplinary Research on Climate (MIROC) model (MIROC3.2-HI) under IPCC Special Report on emissions scenarios B1 (SRB1) scenario (IPCC 2007). The MIROC3.2-HI model, developed at the National Institute for Environmental Studies (NIES) of Japan, had the highest spatial resolution of approximately 1.1° among the selected models in the IPCC Fourth Assessment Report (AR4). The scenario used for the model was created using warming predictions by the Center for Climate System Research/National Institute for Environmental Studies (CCSR/NIES) for modifying regional historical precipitation and temperature observations. The SRB1 includes improved emission baselines and the latest information on economic restrictions throughout the world. SRB1 examines different rates and trends of technological change and expands the range of various economic development pathways, including narrowing of the income gap between developed and developing countries (IPCC 2007).

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 AR4 SRB1 scenario (IPCC 2007) for the study site are given in Table 1 and Figure 4. According to the above procedures, we calculated average data for short-term, mid-term and long-term climate changes for the following model simulation.

Hydrology model Soil and Water Assessment Tool is a hydrologic and water quality model developed by the United States Department of Agriculture Agricultural Research Service (Arnold & Allen 1996; Narasimhan et al. 2005; Gosain et al. 2006; Abbaspour et al. 2007; Yang et al. 2007; Schuol et al. 2008). The main objective of SWAT is to predict the impact of land-use management on water, sediment and chemical yields 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. Water yield predicted by the SWAT is defined as the net water amount leaving the sub-basin that contributes to streamflow. Water yields on land include surface, lateral, and groundwater in each HRU.

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. Major input data include climate data, a terrain map, soil properties and a land-use map. The climate data were obtained from the Japan Meteorological Agency (JMA). Another requirement of spatial input data for SWAT is a digital elevation model (DEM) and maps of land-use and soil type. The DEM was acquired from the Japanese Geographic Survey Institute (GSI; 50 m × 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; 100 m × 100 m). SWAT was initially set up directly with streamflow during 2001–2002, 2003–2007 and 2008–2009 as warm-up, calibration and validation periods, respectively.

Systematic conservation planning is a widely used approach for designing protected area systems and ecological networks (Kirkpatrick et al. 1983; Nichols & Margules 1993; Humphries et al. 1996; Ball & Possingham 2001; Harwood & Stokes 2003). The MARXAN model is a systematic conservation model, which uses an optimization algorithm to identify protected area planning units that achieve conservation targets while minimizing conservation costs (Pressey et al. 1996; Airame et al. 2003; Cook & Auster 2005). The MARXAN model allows the user to vary many aspects of the problem: the number and types of conservation ecosystem services such as water yield included in the analysis; a target for each conservation ecosystem service; how important it is to meet targets for these conservation ecosystem services; the status of planning units; and the cost of each site that could be in the reserve network (Possingham et al. 2000; Culver et al. 2004; White & Kiester 2008).

The MARXAN has been 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). The evaluation of multiple scenarios is one of its major strengths, so we took the MARXAN model as an effective tool to construct priority conservation areas for water yields across different periods with a set of conservation targets under a climate change scenario. This study used the MARXAN model (Ball & Possingham 2001) to perform simulated annealing; a Monte Carlo procedure was used to minimize multivariate functions associated with a physical or biological system (100 runs with 1,000,000 iterations).

In the MARXAN analysis of water yield, the MARXAN can consider the distribution of water yield, connectivity responses, variable emphases on the highest local quality of water yield, and cost efficiency for water resources conservation planning. The cost can effectively be used to reflect corridors that connect protected areas with one another, or water yield ecosystem service with protected areas, by lowering the cost value of the planning units within these identified corridors (hereafter costs are the reflection of areas in this study). Therefore, we analyzed integration of water yields across different periods into the priority conservation planning. Water yields used in the MARXAN model were obtained from each HRU in the SWAT model, and all calculation in the MARXAN was based on planning unit hexagons (the area size is 100 ha). The basic planning unit hexagons can produce more efficient results, where an efficient result is defined as the extent to which priority areas minimize costs while also meeting targets. The detailed process is as follows: (i) the spatial distributions of simulated water yields of HRUs (shape file) were established based on simulated results of the SWAT model; (ii) the spatial distributions of simulated water yields of HRUs (shape file) were changed into grid-based data using convention tool in the ArcGIS (grid file); (iii) the grid-based data sets of water yields were then changed into planning unit hexagons data by spatial analysis tool (shape file); (iv) the planning unit hexagons data sets were used as input data of MARXAN; and (v) the MARXAN calculated optimal priority conservation areas and selection times for water yields through the whole watershed.

The whole Teshio River watershed was divided into 61 subwatersheds and 648 HRUs. The distribution of the measured average annual water quantity and quality in the calibrated and validated periods was close to the simulated values. For calibration and validation of the SWAT model, NSE and R² values for streamflow were greater than 0.5, which indicate a satisfactory mode. Hence, the simulated results are credible and favorable, suggesting that the SWAT model has outstanding applicability in Teshio watershed.

The calibrated and validated SWAT model was used to simulate water yields under climate changes. The average monthly change of hydrologic components under baseline climate conditions is shown in Figure 5. Annual average precipitation in the study period was 959 mm, and the highest value occurred in October in the form of rainfall. Annual average snowfall was 252 mm, occurring from October through April, and its maximum value appeared in February. Annual average surface runoff was 234 mm, with maxima in April and May derived from spring snowmelt. The lowest surface runoff rates were in January and February when the snowpack was deepest and temperature was coldest. The water yield calculated from the SWAT model within this study watershed include surface runoff and groundwater flows. Annual average potential evapotranspiration (PET) in the Teshio River watershed was 459 mm year⁻¹, and the highest value occurred in July. Particularly, average ET from May through August was 34.5 mm month⁻¹ and less than PET. The highest value of actual evapotranspiration was simulated for July and August. The greatest difference between PET and ET appeared in June and July. Based on the information on climate changes, annual average precipitation increased from 958.7 mm under the baseline climate conditions to 1,099 mm under the long-term climate change scenario (Figure 6). There was a notable reduction of annual average snowfall impacted by warmer climate, because the temperature increased in winter (Table 1).

All climate changes caused greater monthly precipitation than that of the baseline climate conditions, especially in September, but May and October had different tendencies (Figure 7(a)). However, all climate changes decreased monthly precipitation in October over that of the baseline climate conditions. The maximum value of decreasing snowfall was from October through April for the long-term climate change, in which average temperature increased by 3.88 °C (Figure 7(b) and Table 1). Snowfall also decreased in November, December and March for all climate changes. Snowfall in October and April decreased to zero in the mid-term and long-term climate changes, and changes of snowfall in January and February were not apparent. This suggested that the snow-cover-free period increased from 5 months (May through September) under baseline climate conditions to 7 months (April through October) under future climate changes (Figures 7(b) and 5).

Monthly average surface runoff increased under three climate changes in January, February, March, August, September, November and December, with the peak value in March (Figure 7(c)). Monthly average surface runoff decreased in April and May under all climate changes, and changes of monthly average surface runoff in other months were not obvious. PET and ET increased under all climate changes, especially during snowmelt season (March and April) and summer (June and July) (Figure 7(d)). This indicates that increased temperature could accelerate the evapotranspiration of vegetation in the watershed. ET was less than PET because ET was restricted by soil water availability, especially from April through August. Under baseline climate conditions, PET was estimated at 1.2 times greater than ET in March and April. Under the long-term climate change scenario (3.88 °C increase of air temperature), PET was double ET in March and April. Despite changes in the seasonal distribution of ET, the change of annual average ET was relatively small. Annual average ET under long-term climate change was 1.14 times higher than that under baseline climate conditions.

Average monthly water yields increased in January, February, March, September and December under all climate changes, whereas changes in other months were not transparent (Figure 7(e)). Annual average surface runoff, lateral flow, groundwater, ET, and water yield varied with multiple climate change conditions (Table 2). All the hydrologic components in Table 2 were increased by climate change relative to the baseline climate (long-term > mid-term > short-term > baseline). As a whole, annual average water yield under long-term climate change was 1.16 times higher than that under the baseline climate conditions.

In this study, the water yields in February, April and October can represent typical snowpack, snowmelt and rainy seasons in the Teshio River watershed (Figure 5). We took the water yields across these different periods under long-term climate change as an example to explore their spatial characteristics and construct spatially optimal priority conservation areas for water resources, because the largest values of water yields occurred under long-term climate change. The water yields across different periods under long-term climate change fluctuated seasonally and spatially (Figure 8).

The water yield in February ranged from 8.13 to 29.33 mm across the Teshio River watershed, and the higher value in the whole of the landscape occurred in the southwestern and riverine regions (Figure 8(a)). The water yield in April ranged from 39.89 to 179.67 mm across study watershed, with the largest value in the whole watershed located in the northern area (Figure 8(b)). The water yield in October ranged from 37.74 to 132.37 mm across the study catchment, and the largest value in the whole watershed appeared in the southeastern area (Figure 8(c)). The annual average water yield ranged from 437.33 to 833.66 mm in the Teshio watershed, with the highest value in the northern regions (Figure 8(d)). In order to discuss the importance of water yield in April, we further calculated the composite water yield which was obtained from summation of water yields in February, April, and September. In this study, we chose the largest values of each water yield in different periods as the potential conservation targets of water resources conservation planning.

We chose 10% to 90% of the largest values of water yields corresponding to different typical periods (water yields in February, April and October, annual average water yield) in Figure 7(e) as the conservation target levels according to the agricultural and forest land priority conservation areas for all water yields in typical periods under climate changes. The agricultural and forest land priority conservation areas for water yield ecosystem services increased almost linearly with the conservation targets (Figures 9 and 10). The varied trends and ranges of forest land priority conservation areas for all water yields in typical periods were more regular than those of agricultural land priority conservation areas because forest is the dominant land use type in Teshio River watershed which determines the main hydrology and ecological processes.

The overall trends in agricultural land priority conservation areas for water yields were characterized by volatile change processes because the agricultural lands in this study site showed complex interlaced distributions. Agricultural land priority conservation areas for water yield in February increased almost linearly with conservation target levels until conservation targets reached 20%, 40% and 70% target levels under baseline climate conditions, short-term, and mid-term climate changes, respectively (Figure 9(a)).

For the water yield in April, there were two obvious groups: agricultural land priority conservation areas under baseline climate conditions and short-term climate change, and agricultural land priority conservation areas under mid-term and long-term climate changes based on the variation in trends and shapes of curves (Figure 9(b)). Agricultural land priority conservation areas for water yield in October overlapped as the conservation targets increased from 10% to 40% target levels under all climate changes (Figure 9(c)). When the conservation targets continued to increase, the curves for agricultural land priority conservation areas were upward convex under baseline climate conditions and short-term climate change, and higher than the other two climate changes. Agricultural land priority conservation areas increased with conservation target levels until conservation targets were 60% and 70% target levels under baseline climate conditions and short-term climate change, respectively. The curves for agricultural land priority conservation areas were slightly downward convex under mid-term and long-term climate changes. For composite water yield (including summation water yields derived from February, April and October), the curves for agricultural land priority conservation areas were slight upward convex under all climate changes (Figure 9(d)).

Agricultural land priority conservation areas for water yield ecosystem services increased almost linearly with conservation target levels until the conservation target was 40% target level under short-term climate change. There was an intersection point of curves for agricultural land priority conservation areas at 60% conservation target under all climate changes. Agricultural land priority conservation areas for annual average water yield increased as the conservation target increased from 10% to 90% under all climate changes (Figure 9(e)). The curves for agricultural land priority conservation areas overlapped each other until the conservation target was 60% target level.

The spatial patterns of forest land priority conservation areas were different from those of agricultural land priority conservation areas, and the distances among these curves for forest land priority conservation areas were smaller (Figure 10). Forest land priority conservation areas for water yield in February increased almost linearly with conservation target levels until conservation targets were 20%, 40% and 70% target levels under baseline climate conditions, short-term, and mid-term climate changes, respectively (Figure 10(a)). The characterized tipping points that occurred in forest land priority conservation areas were similar to those that appeared in agricultural land priority conservation areas. For water yield in April, forest land priority conservation areas increased almost linearly with conservation target levels under baseline climate conditions and short-term climate change (Figure 10(b)). The curves for forest land priority conservation areas were slight upward convex under mid-term and long-term climate changes, and higher than the other two climate changes. Forest land priority conservation areas for water yield in October increased almost linearly with conservation target levels under mid-term and long-term climate changes (Figure 10(c)). The curves for forest land priority conservation areas under baseline and short-term climate changes overlapped each other and remained unchanged when the conservation target increased from 70% to 90% levels. When the conservation targets increased from 80% to 90% levels, the curves for forest land priority conservation areas overlapped each other and remained unchanged under baseline climate conditions, short-term and mid-term climate changes.

For composite water yield, forest land priority conservation areas slowly increased as the conservation target increased from 10% and 40% levels under baseline climate conditions (Figure 10(d)). Then the forest land priority conservation areas increased almost linearly with conservation target levels (from 40% to 60% levels) and finally remained unchanged. The curves for forest land priority conservation areas overlapped each other when the conservation targets increased from 70% to 90% levels. Forest land priority conservation areas for annual average water yield increased almost linearly with conservation target levels under all climate changes (Figure 10(e)). There were subtle differences among these curves for forest land priority conservation areas; the curve for forest land priority conservation area under baseline climate conditions was higher than the other three climate changes.

According to the previous results, there were three obvious tipping points in the agricultural and forest land priority conservation areas for water yields with 20%, 40% and 60% maximum total amount target level under climate changes, respectively. The spatial priority conservation areas for water yields in February, April and October and annual average water yield were closely related to spatial distributions of their corresponding water yields (Figures 8 and 11). The places with higher water yield amounts would be a priori involved to reserve networks. Furthermore, we took the optimal priority conservation areas for water yields with 40% maximum total amount target as an example to explore their spatial distribution characteristics and interactive relationships in detail (Figure 11). Spatial priority conservation areas identified for water yields in February, April and October, composite and annual average water yield with 40% maximum total amount target level under long-term climate change were 28%, 56%, 30%, 56% and 32% of landscape, respectively (Figure 11(a)).

There were similarities between the priority conservation areas spatial pattern and the water yield spatial pattern. The spatial pattern of priority conservation areas for water yield in April was close to that of composite water yield, suggesting the water yield in April played an important role in water supply in this watershed. The Teshio catchment is a winter-dominated watershed and is covered by snowpack for about half the year; in particular, the water yield in April derived from snowmelt contributed much more to the total annual average water yield. Although there were the same spatial priority conservation areas for water yields in typical periods to some extent, the selection times (ranking) of priority conservation areas were different (Figure 11(b)). There were some differences among the selection times of water yields across different typical periods because their spatial distributions were totally different. The selection times in the southwestern (agricultural land and populated area), southeastern (headwater of this watershed), and northern (agricultural land) areas were always higher than in the other areas through comprehensively taking into account water yields in typical periods. This indicates that these areas needed to have priority for conservation because of the high water supply for hydropower electricity production, irrigation for crop production and resident water use in Teshio River watershed.

Priority conservation areas identified for water yields in February, April and October, composite and annual average water yields with 20% maximum total amount target level under long-term climate change were 12%, 24%, 14%, 24% and 15% of landscape, respectively. Moreover, priority conservation areas identified for water yields in February, April and October, composite and annual average water yields with 60% maximum total amount target level under long-term climate change were 46%, 92%, 47%, 92% and 59% of landscape, respectively.

The overlay analysis of priority conservation areas for water yields in different periods with 10% to 90% conservation targets under long-term climate change also showed that core protected areas for water yields occurred in the southwestern, southeastern and northern areas (Figure 12). This useful information complemented selection times of conservation areas, and the value indicates how many times the planning units are involved in priority conservation planning under different targets (e.g. the value of 9 showed that the planning units were included nine times in conservation areas from 10% to 90% target levels). As the conservation targets increased, the total boundary length of priority conservation areas increased (Figure 13), which revealed that the number of conservation patches increased and became more fragmented. The pair-wise Jaccard's indexes and correlations of selection times of optimal priority areas for water yields across different periods with 40% conservation targets under long-term climate change are shown in Table 3. The Jaccard's indexes of water yields in February and April and water yield in October were low (0.11–0.15), the Jaccard's indexes of water yields in February, April and October and annual average water yield were medium (0.27–0.42), and the Jaccard's indexes of water yields in April and composite water yield were high (0.91). The correlations between water yields in February and April and water yield in October were negative (−0.0030 and −0.38, respectively). The correlation between water yield in April and composite water yield was positive and highest (0.99).

In this study watershed, the water yield in February increased from 1.68 to 13.89 mm under climate changes, owing to increasing precipitation in the form of rainfall caused by an increase in temperature (more snowfall shifted into rainfall). Bouraoui et al. (2004) also found a similar increase in water yield in winter under climate change. The highest water yield was located in the southwestern and riverine regions (Figure 8(a)).

Spatial variability of climate factors including precipitation and temperature, soil and vegetation causes spatial variability of the water balance (Brown et al. 2005). The water yield in February was mainly provided from the melting water as baseflow because of the lowest temperature and thickest snowpack, which is a typical phenomenon in snow-dominated cool regions (e.g. Marshall & Randhir 2008). In this area, the land is covered by farmland and paddy fields, but there are no crops growing in winter (Figure 3(a)). The dominant soils were brown lowland and gray lowland soils which can easily saturate because of low field capacity water contents and flat topography (Figure 3(b) and 3(c)). It was suggested these characteristics of land cover (no crop), soil (low field capacity) and topography (flat) create a spatial increase of water yield in the southwestern region.

The peak water yield occurred in April under the baseline climate conditions, short-term and mid-term climate changes. The water yield in April increased from baseline climate condition to short-term climate change, but decreased in mid-term and long-term climate changes. Under long-term climate change, the maximum value of water yield was in March, earlier than in the other scenarios. This suggests that climate changes shifted the snowmelt peak from April to March and further decreased water yield in April. The other reason is that there was a marked reduction of annual snowfall under a warmer climate owing to the temperature increase during winter and snow-cover-free period increase, which also decreased water yield from the snowmelt process. Early snowmelt could change the annual distribution of the water balance for the watershed, which could cause a further decline in water yield during late spring and summer months (Marshall & Randhir 2008; Fan & Shibata 2016).

The maximum value of water yield in April occurred in northern areas (Figure 8(b)). The water yield in April was mainly originated from snowmelt, contributing the peak flow during the snowmelt period. The water yield in western parts was higher than in other areas, which corresponded to the land covers of pasture and farmland, and soil types of brown lowland and peat soils. It was suggested that higher precipitation in winter (Wassamu, Shibetsu and Bifuka in Figure 2) and deeper snowpack in the western and northern areas contributed higher water yield by snowmelt in these areas in April. The larger value of water yield in October appeared in southeastern areas (Figure 8(c)).

The precipitation in the north and southeast was higher than in the other areas, which provided the largest water yield in the southeast region where the original water source of Teshio watershed is located. The slopes were steeper in the southern area than in the other areas, suggesting that the steeper topography can enhance the rapid response of water yield to the autumn storm rainfall (Figure 3(c)). The steeper slope might also increase the soil water availability in the lower slopes, causing rapid water saturation in soil during a storm event. Thus, the spatial redistribution of surface runoff due to the topographic factor would result in higher soil water availability on lower slope positions, which contributes to the spatial water yield gradient in slope (Tsukamoto & Ohta 1988; Tani 1997). The higher value of annual average water yield was located in the northern regions (Figure 8(d)). It was suggested that lower temperature and evapotranspiration and higher precipitation (including the rainfall and snowfall) contributed to the higher annual average water yield in the northern area.

The determined 20%, 40% and 60% targets could reflect the spatial characteristics of priority conservation, especially for watershed ecosystem conservation involving more efficient agricultural and forest lands (Egoh et al. 2009; Fan & Shibata 2015). The selected typical targets could be used as the reference values on conservation thresholds for ensuring high water supply in the Teshio River watershed. Results of this study suggest a great deal of flexibility in the selection of priority conservation areas for equal representation of water yields in typical periods in the Teshio watershed. Not only do conservation target levels of representation of water yields across different periods, expressed in terms of agricultural and forest land priority conservation areas, correlate strongly with the percentage of the study area needed to represent them in a minimal amount of space, but these areas remained fairly constant under a variety of conservation targets (Figures 9 and 10). These patterns were consistent over a range of representation levels as well.

Many planning units that were selected for the most optimal or efficient solutions for water yields across typical periods were located in the southwest, southeast and north of the study site, which also demonstrated a great deal of overlap across a variety of target levels (Figure 11). Areas that recurred in the outcomes for the range of target levels for water yields across different periods represented a good starting point for water resources management in Teshio watershed design discussion. There is less agreement about optimal conservation areas and selection times between the smaller and larger conservation target levels. The most optimal outcomes for the smaller water yield conservation targets exhibited more spatial clustering than for the larger target levels. This may be attributed to a combination of conservation cost surface (planning unit area in this study) and the inherent challenge of the optimization method to find the lowest-cost outcome while meeting targets. As conservation targets increased, it became more difficult for MARXAN to find both the lowest cost and most compact outcome using the planning unit area cost surface as a basis.

In this study, each of the 100 runs of a simulation represented an optimal or near-optimal outcome for a set of representation targets given constraints including the area cost and connectivity defined by the optimal function in the MARXAN model. The summed outcomes (selection times) provided an indication of the amount of inflexibility and priority ranking in number and location of these outcomes.

Measuring the ecological, economic and climate effectiveness of a protected area network is a major challenge in systematic conservation planning (Margules & Pressey 2000; Naidoo et al. 2006), but the priority conservation areas and selection times are widely used as surrogates to measure the efficiency and ecological viability of reserve network designs (Kark et al. 2009; Grantham et al. 2010; Weeks et al. 2010).

In this empirical study, we also used the pair-wise Jaccard's index and correlations (corresponding to conservation areas and selection times, respectively) of optimal priority conservation areas of water yield ecosystem service under long-term climate change to explore relationships among water yields across different periods (Khan et al. 1997; Chape et al. 2005). The results suggested that there were tradeoffs between water yields in February and April and water yield in October owing to lower Jaccard's indexes and negative correlations. Understanding tradeoffs of priority conservation areas of individual ecosystem services through Jaccard's index and correlation of selection times could promote construction of highly compact priority areas of a watershed conservation system to satisfy conservation targets. These two indicators also could quantitatively assess the spatial concordance among water yields across different periods on a regional scale. Finding out how water yields across different periods interactions change as climate and conservation targets change may help avoid unnecessary losses by focusing on finding the most efficient solutions to mitigate the tradeoffs or to enhance synergies (Moilanen & Wintle 2007; Smith et al. 2010). Identifying the tradeoffs and synergies among water yields in typical periods is likely to improve water resources management practices and strengthen the decision-making processes to achieve specific conservation objectives (Carreno et al. 2012).

The tradeoffs are evident in the difference of spatial pattern and priority conservation areas for water yields across different periods under climate change, where priority conservation areas for water yields in February and April, and composite water yield were traded off against those for water yield in October under minimum protected areas constraints, and highlighted a final cautionary note in water resources conservation planning (Figures 11 and 12; Table 3). Planning of composite water yield involving water yields in February, April and October can maximize the priority conservation areas for water yields in February and April, but the implementation of such plans may be challenging for conservation of water yield in October or vice versa. The spatial distribution characteristic of priority conservation areas for composite water yield was very close to that of water yield in April, indicating that the water yield in April significantly affected composite water yield conservation on a typical seasonal or monthly basis. Meanwhile, the priority conservation areas for annual average water yield could involve the important places with higher water yields across different periods to the greatest extent. Both tradeoffs and synergies of water yields across different periods should be managed to either reduce their associated costs to society or enhance water resources multifunctionality and net human welfare. Systematic conservation that places planning units into management zones of differing costs, actions and contributions to targets will be valuable in recognizing the different management requirements of water resources across the watershed under climate change.

The effectiveness of a fixed water resource system of reserves may be compromised under climate change, a problem that may be compounded when water yields fluctuate and they are readily impacted by geographical, soil, vegetation and meteorological driving factors from temporal and spatial scales (Guo et al. 2000; Somura et al. 2009; Fan & Shibata 2015). Therefore, the spatial distribution of priority conservation areas for water yields across different periods that coincidentally overlap under current climate may diverge in distribution with a set of target levels under climate changes.

Our results showed that identifying optimal priority conservation areas based on knowledge of highly localized water yield amounts, total boundary length and minimum conservation cost could capture core areas of water yields' distribution. These results are one source of qualitative and quantitative insight on regional priority conservation of water resources in the face of climate change through integrating representation of water yields across different periods. In this case, the systematic conservation analysis is used for incorporating changes in hydrological and ecological interactions into conservation prioritization of water yield ecosystem service at temporal scales.

Spatial explicit formulations of these interactions have been promoted by coupling a watershed hydrology model with the MARXAN model (Fan & Shibata 2016). The outcomes produced by the MARXAN model might build on simulated water yields from a hydrology model because we could plausibly evaluate outcomes based on contraction or expansion of reserve areas for water yields in typical periods. Integrating future climate changes into the hydrology model approach used here also could predict current and future water yields based on correlations between complex hydrology mechanisms and climatic variables (Bouraoui et al. 2004; Somura et al. 2009). The proposed research framework makes it possible to value areas not only by the amount of water yield ecosystem service within the area, but also by its connectivity, which could induce connectivity of good quality priority conservation areas in the reserve network of water yields across different periods to account for spatial distribution shifts due to climate change. This framework is a key component in improving the efficiency of localized water yield service conservation through multiscale conservation strategies focusing on both sites and large reserves (Moilanen et al. 2009). Therefore, conservation planners should incorporate temporal and spatially explicit interactions of hydrology cycling into water resources conservation management decisions across the watershed under climate changes.

We made the simplifying assumption that the ‘cost’ of planning units was directly proportional to their size. So we were easily able to directly compare tradeoffs in total priority conservation areas and spatial structure of outcomes equivalent in water yields representation. In future analyses, the cost function could be expanded to include additional social or economic terms, such as those related to commercial and recreational water yields, with appropriate values assigned to individual planning units. Weighting factors could be applied to these terms in order to evaluate the effects on outcomes and the tradeoffs associated with one or more of these factors. Use of some known water resources reserve areas could aid in the perfect and compact outcomes of water yield ecosystem service. For instance, known conservation patches (planning units) that are important for sustaining key water yields (e.g. sustainable water supply) can be identified a priori to be included in the outcome. Such planning units can be selected for water yields of ecological or socio-economic importance as well as those that are sensitive to disturbance or are of particular interest to diverse stakeholders.

Additionally, water resources management problems across the watershed are often characterized by variability over temporal and spatial scales, and sometimes provoke conflict among diverse sets of stakeholders. Our approach here produced an array of planning opinions, but without proper stakeholder engagement, so the outcomes derived from MARXAN may be complementary to a formal stakeholder-driven water resources planning process. In future research, we should work through planning processes with diverse stakeholders by allowing them to become more involved in using MARXAN in a collaborative environment, stimulate discussion, and allow for their own formulation of new socio-economic cost surfaces and multiple objective conservation targets for water yield across different periods under climate changes.

This study provided a framework that integrated a catchment-scale hydrology model to systematic conservation planning of water yield ecosystem service under climate changes. We analyzed the impact of climate changes on hydrology and spatially optimal conservation areas for water yield in the Teshio River watershed of northern Japan. In particular, the water yields in February (snowpack), April (snowmelt), October (rainy seasons) and annual average water yield may gradually increase under climate changes, but the water yield in April only increased under short-term climate change and gradually decreased under mid-term and long-term climate changes. The larger value of water yield in February was concentrated in the southwestern and riverine areas. The higher value of water yield in April was concentrated in the north of the study watershed because the snow melting produces much greater water yield. The higher water yield in October was concentrated in the southeastern areas.

Climate change provides different temporal and spatial patterns in the water yield ecosystem service, thus affecting the optimal priority conservation areas for water resources at spatial and temporal scales. The spatial priority conservation areas for water yields across different periods were concentrated in the northern, southwestern and southeastern areas of Teshio River watershed under climate changes. The total forest and agricultural lands priority areas for water yields increased with conservation target levels, and the spatial structures of priority conservation areas for water yields in typical periods with smaller conservation targets were more compact than those with higher target levels under climate change. There were tradeoffs and synergies among spatial priority conservation areas for water yields across different periods. The temporal and spatial modeling and assessment of water yield integrated with the systematic conservation model may facilitate a better understanding of the climate change challenge. The conservation planning options illustrate how competing needs might be balanced in planning for water yields across different periods under climate change in the Teshio River watershed and are meaningful to guide stakeholders and decision-makers for future water resources planning activities.

This study was supported by National Natural Science Foundation of China (No. 41601088). The work was also supported by Southwest University of Science and Technology (No. 17LZX690; No. 18LZX619). The work was conducted under 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.