Hotspot identification and interaction analyses of the provisioning of multiple ecosystem services: Case study of Shaanxi Province, China

Ecological Indicators 107 (2019) 105566

Contents lists available at ScienceDirect

Ecological Indicators journal homepage: www.elsevier.com/locate/ecolind

Hotspot identification and interaction analyses of the provisioning of multiple ecosystem services: Case study of Shaanxi Province, China

T

⁎

Lulu Liua,b, Hanbing Zhanga, Yang Gaoa,b, , Wenjie Zhua, Xing Liuc, Qiaodi Xua a

College of Land Science and Technology, China Agricultural University, Beijing 100193, China Shaanxi Key Laboratory of Land Consolidation, Shaanxi 710075, China c College of Resources and Environment, Southwest University, Chongqing 400716, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Ecosystem services Trade-offs/synergies Multifunctionality Hotspot identification Shaanxi Province

The hotspot identification and interaction analyses of the provisioning of multiple ecosystem services (ESs) are basic prerequisites for sustainable ecosystem management. ESs and their interrelationships are determined by ecological processes over time and space. The high-value ES areas is crucial for regional land use management. However, most of hotspots identification are limited to a single point in time and/or in space. Here we proposed a framework for identifying a long-term stable supply area for multiple ESs, namely, an ES multifunctional area (ESMA) to recognize the region that can stably provide high-value ESs. Using Shaanxi Province as an example, four ESs, i.e., net primary productivity (NPP), soil conservation (SC), water yield (WY), and habitat quality (HQ), were quantified from 2000 to 2015 at a pixel scale of 1 km. Then, we identified ESMA in Shaanxi and explored the trade-offs and synergies among ESs within and outside ESMA. Results of the study provided the following conclusions: (1) With the exception of WY, ESs exhibit apparent spatial differentiation in areas bounded by the Qinling Mountains due to terrain and climate differences between the north and south of these mountains. (2) The ES supply areas of NPP accounted for 41.47% of the total area of Shaanxi Province. This finding indicates the NPP provision was relatively stable. (3) The ESMA of Shaanxi Province is largely situated in southern Shaanxi. The primary land use types of ESMA are forest and grassland, and the areas that provide NPP, HQ, and SC services are the largest among all service combinations. (4) The trade-offs and synergies of ESs are more stable in ESMA than those in non-ESMA. This study proposes to provide a long-term provision of multiple stable highvalue ESs to a priority protected area, thereby offering support and services for the development of ecosystem differential spatial protection policies.

1. Introduction Ecosystem services (ESs) refer to natural conditions and utilities provided by ecosystems and ecological processes that can sustain human life (Daily, 1997). They not only provide human beings with food, medicine, and other raw materials for production and maintenance but also create and maintain Earth’s ecological life-support system and develop environmental conditions necessary for human survival (Fisher et al., 2009). The Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES) reports that ecosystem health is deteriorating at an unprecedented rate, undermining the foundations of our livelihoods and global quality of life (IPBES, 2019). Against this background, the protection and management of ESs is an indispensable link for the sustainable development (Vihervaara et al., 2010; Qiao et al., 2019). Hotspot areas of ESs can

⁎

provide high-level services both on-site and off-site. Therefore, hotspots identification and interaction analyses of the provisioning of multiple ESs are of great significance to enhance the human well-being. Recently, enhancing the benefits from ES synergistic promotion has become a popular topic (Foley et al., 2005; Wu et al., 2013). A series of excellent results has emerged to explore how ESs correlate with one another and to search the hotspots of ES supply for conservation. Raudsepp-Hearne et al. (2010) identified six types of ES bundles by analyzing the spatial patterns of 12 ESs in a mixed-use landscape that consisted of 137 municipalities in Quebec, Canada. They demonstrated that landscape-scale trade-offs between provisioning and nearly all regulating and cultural ESs. Palomo et al. (2013) used the participatory survey method to analyze the supply and benefiting areas of the ecological supply, regulation, and cultural services of national parks in the southwestern coast of Spain. They also developed a conceptual model of

Corresponding author at: College of Land Science and Technology, China Agricultural University, Beijing 100193, China. E-mail address: [email protected] (Y. Gao).

https://doi.org/10.1016/j.ecolind.2019.105566 Received 31 March 2019; Received in revised form 10 June 2019; Accepted 15 July 2019 1470-160X/ © 2019 Elsevier Ltd. All rights reserved.

Ecological Indicators 107 (2019) 105566

L. Liu, et al.

service space flow based on the spatial location of the service supply and benefiting areas and proposed to protect hotspots of ES provision. In terms of research methods, descriptive analysis methods (e.g., spatial cartographic analysis) and statistical analysis methods based on mathematical theory (e.g., regression analysis) are widely used in the study of complex interactions among multiple ESs (Qiu and Turner, 2013; Raudsepp-Hearne et al., 2010; Dobbs et al., 2014). The trade-offs/synergies of ESs have gradually evolved from qualitative to quantitative with in-depth research. Moreover, methods for quantifying the differences between an individual ES and the average ES, such as the rootmean-square error (RMSE), have been developed to describe the magnitude of scattering from the mean. However, ESs and their interrelationships are determined by ecological processes over time and space (Rodriguez et al., 2006). In existing studies, the analysis of the relationship between ESs is frequently limited to a particular point in time (Hou et al., 2017; Tilman, 2000; Qiao et al., 2019), and few empirical papers on ESs are dealing with temporal dynamics (Rau et al., 2018). In addition, most studies of ESs supply focus on mapping the supply of ESs at a single point in time and/ or in space (Lopes and Videira, 2017). However, ESs supply can strongly vary at the time frame of decades (Hein et al., 2016). The identification of high-value ES areas, which is crucial for land use management (Hamann et al., 2015), will lack reliability and accuracy, if it is recognized only based on certain time point (Jaligot et al., 2019). In addition, the interactions of ESs are also interrelated with regional scope. The spatial heterogeneity of ESs for distribution is attributed to the differences in ecological background and resource conditions within a varied geographical range, which will further affect the relationships among ESs (Dai et al., 2016). A long-time series analysis can help us to obtain the stability of the interactions between ESs, then promoting our understanding of the potential driving factors (Raynolds et al., 2008; Hao et al., 2019). Therefore, studying the relationship among ESs based on a certain time span and geographical division is critical to break through the limitations of static studies and achieve better management of regional ecosystems (Han et al., 2017). In the current study, we propose the concept of multifunctional areas of ESs (i.e., ESMA) to explore hotspots of ES supply from the spatiotemporal perspective. These hotspots are defined as areas that can provide a variety of high-value ESs over a long period. From the concepts of multifunction and ESs, ESMA is the inevitable result of a stable spatial cooperation among diverse ESs (Vejre et al., 2007; De Groot et al., 2002), and its identification is based on the spatial heterogeneity of ESs within a certain time span rather than a simple overlay analysis of a single year. It should be noted that, different from landscape multifunction, ESMA emphasizes on the areas where natural ecosystem can provide products and services to meet human needs from the perspective of donor-side. Landscape is a mix and repeated of local ecosystem or land use types over one land (Forman, 1995). Landscape multifunctionality is often conceived and evaluated as the combined supply of multiple ESs at the landscape level (Mastrangelo et al., 2014), which emphasizes human use as well as regional natural locality and attributions. The conservation of ESMA, as an area that can steadily supply multiple ESs, is an effective strategy to meet the unlimited development demands of humans within a limited land and to increase the benefits obtained by humans from ESs. Furthermore, considering interactions among ESs and conducting targeted ecological management are likely to produce considerably better outcomes for societies because human management preferences and decisions can change the type, scale, and relative combination of services provided by ecosystems (Swinton et al., 2007; Carpenter et al., 2009). Consequently, we also investigated the similarities and differences between the trade-offs/synergies of ESs inside and outside the boundaries of ESMA. We aim to determine whether the characteristics of the relationships among ESs in ESMA differ from those in other places. Moreover, we intend to guide local governments in efficiently protecting regional ecosystems and in formulating

differential spatial protection policies. As a typical region with obvious spatial distribution of climate, vegetation and topography, Shaanxi Province, China, has representative of the spatial distribution and hotspots of ESs. Therefore, we selected Shaanxi Province, located in Northwest China, as an example. We identified its ESMA by selecting four key ESs that consider local conditions, namely, net primary productivity (NPP), soil conservation (SC), water yield (WY), and habitat quality (HQ). Then, we explored the temporal variation and spatial scale dependence of the trade-offs and synergies among ESs within and outside ESMA. We focused on the following issues: (i) quantifying the spatial distribution of regional critical ESs, (ii) identifying ESMA based on spatiotemporal stability, and (iii) comparing and analyzing the trade-offs/synergies between ESs within and outside ESMA. This study intends to enhance the efficiency of eco-environment management under the premise of preserving the value of regional ESs. It also aims to provide a decisionmaking basis for land management and landscape planning. 2. Materials and methods 2.1. Study area Shaanxi Province (105°29′–111°15′E, 31°42′–39°35′N), which has an area of 210,000 km2, has a long and narrow span from north to south. Shaanxi is bounded by the Qinling Mountains and exhibits a distinct diversity of natural and geographical features (Fig. 1). On the basis of its geographical environment, the province can be divided into three natural regions, namely, Northern Shaanxi, Guanzhong, and Southern Shaanxi, by the Beishan Mountains and the Qinling Mountains from north to south. The spatial distributions of precipitation, temperature, and solar radiation in the territory exhibit heterogeneity. Mountains, plateaus, basins, deserts, and large and small rivers are interlaced; and the area has temperate, warm temperate, and northern subtropical climates. Altitude is mostly distributed between 500 m and 2000 m. The average annual temperature is 7 °C–16 °C, and the mean annual rainfall is 320–1258 mm. The land use types in the area are primarily farmlands, grasslands, shrubs, sparse woodlands, and woodlands. The hilly area in Northern Shaanxi probably has the highest soil erosion rate in the world (Fu et al., 2000). The Qinling Mountains are rich in mineral resources and high in forest coverage, but are threatened by illegal mining and overcutting. The Guanzhong Plain has a flat terrain, urban agglomeration, and a concentrated distribution of dryland resources. However, its soil is relatively soft and prone to drought. Since 2000, Shaanxi Province has implemented a series of ecological restoration projects, such as the Grain-for-Green Program (GFG) and Natural Forest Protection Program (NFPP). These projects have considerably improved vegetation status in the province, particularly in Northern Shaanxi. However, Shaanxi remains as one of the regions with the most serious ecological and environmental problems, such as soil and water loss (Li et al., 2017), deforestation (Liu et al., 2013), and desertification (Kateb et al., 2013). Moreover, Shaanxi has lost its natural protective and productive functions in China (Hou et al., 2015). To sum up, Shaanxi Province has a unique geographical location and fragile ecological environment with obvious spatial heterogeneity in climate, topography and vegetation types. Taking Shaanxi Province as a research area has typicality and representativeness. Exploring the relative changes of ES hotspots and their trade-off/synergies is essential for the scientific construction of the local land space pattern and the sustainable development of the region. 2.2. Data sources The four key ESs, namely, NPP, SC, WY, and HQ, were quantitatively calculated on the basis of the actual geographical environment and the government’s development policy priorities in the region. The landscape pattern maps used in this study were derived from land use 2

Ecological Indicators 107 (2019) 105566

L. Liu, et al.

Fig. 1. Location and regional division of Shaanxi Province, China.

vegetation greenness obtained through remote sensing satellites and the spatial interpolation of climate data (Field et al., 1995). Combined with the distribution characteristics of terrestrial vegetation in China, we used Zhu’s improved CASA model (Zhu et al., 2007) to calculate NPP in this study.

data from the Resource and Environment Data Cloud Platform (http:// www.resdc.cn/) in 2000, 2005, 2010, and 2015. The area was reclassified into six categories based on research requirements: cultivated land, forest, grassland, built-up land, water body, and other lands. The spatial resolution was 1000 m × 1000 m. The MOD16 evapotranspiration product required by the water equilibrium model and the Moderate Resolution Imaging Spectroradiometer (MODIS) International Geosphere–Biosphere Program (IGBP) land cover data sets required by the Carnegie–Ames–Stanford Approach (CASA) model were produced by NASA (http://modis.gsfc.nasa.gov/). The topographical information used in this study was derived from a digital elevation model (DEM) with a spatial resolution of 30 m × 30 m, which was downloaded from the Geospatial Data Cloud (http://www.gscloud.cn/). Monthly normalized difference vegetation index (NDVI) data with a spatial resolution of 500 m × 500 m were also downloaded from the Geospatial Data Cloud. Meteorological site data were derived from the National Meteorological Information Center and spatialized using the professional meteorological interpolation software Anusplina 4.36. China’s soil maps were obtained from the Harmonized World Soil Database (v1.1), which was provided by the Cold and Arid Regions Sciences Data Center in Lanzhou (http://westdc.westgis.ac.cn). All data were uniformly processed into 1000 m × 1000 m raster data form in accordance with the requirements. Quantitation and spatial analyses were performed using ArcGIS 10.5.

NPP (x , t ) =

∑ [APAR (x , t ) × ε (x , t )]

(1)

APAR (x , t ) = SOL (x , t ) × FPAR (x , t ) × 0.5

(2)

ε (x , t ) = Tε1 (x , t ) × Tε 2 (x , t ) × Wε (x , t ) × εmax

(3)

where NPP( x , t ), APAR( x , t ), and ε (x , t ) represent the NPP (g C/m2), effective absorbed photosynthetic radiation (MJ/m2), and actual light energy utilization rate of pixel x in month t, respectively. SOL( x , t ) represents the total amount of solar radiation (MJ/m2) of pixel x in month t, FPAR( x , t ) represents the proportion at which plants photosynthetically absorb the active radiation of pixel x in month t, and the constant 0.5 indicates the proportion of effective solar radiation accounted for by the total solar radiation. Tε1( x , t ) and Tε2 ( x , t ) refer to the parameters that describe the stress coefficients of the highest and lowest temperatures, respectively. Wε ( x , t ) refers to the parameters that describe the water stress coefficients, which reflect the effects of water conditions. εmax refers to the possible efficiency of different vegetation types under the ideal condition. Potter et al. (1993) and Field et al. (1995) considered that the maximum light energy utilization rate of global vegetation was 0.389 g C/MJ, without distinguishing vegetation types. However, as a physiological property of vegetation itself, the maximum light energy utilization rate of different vegetation types varies. In addition, the value of the maximum light energy utilization rate exerts a considerable influence on the estimation results of NPP. Therefore, this study used the simulation results of the maximum light energy utilization rate of typical vegetation types in China determined by Zhu et al. (2006).

2.3. Methods 2.3.1. NPP NPP refers to the organic matter that is fixed by plants primarily through photosynthesis, and thus, it can reflect the growing status of vegetation and measure the amount of trophic energy flows in food webs and chains (Odum and Barrett, 1971). NPP is an important part of the surface carbon cycle that not only directly reflects the production capacity of vegetation communities under natural environmental conditions and characterizes the quality of terrestrial ecosystems but is also a key factor for determining ecosystem carbon sources/sinks and regulating ecological processes (Field et al., 1998). The CASA model is a terrestrial ecosystem evaluation model that is composed of inversion

2.3.2. SC Soil loss is a major obstacle that restricts regional development in Shaanxi. In this study, the revised universal soil loss equation was used to estimate the potential and actual soil erosion, and the difference was 3

Ecological Indicators 107 (2019) 105566

L. Liu, et al.

used to measure SC service. The equation is expressed as

Dxj = SC = R × K × LS × (1 − C × P )

(4)

Ym

M

∑m =1 ∑

y=1

⎛ ωm ⎞ ⎜ ∑M ω ⎟ m y imxy βx Sjm m m = 1 ⎝ ⎠

(6)

where SC represents SC capacity in t/(ha year); R is the rainfall erosion factor in (MJ mm)/(ha h year); K is the soil erosion index in (t ha h)/(ha MJ mm); LS is the topographic factor (dimensionless); C is the vegetation cover factor (dimensionless, ranges from 0 to 1); and P is the conservation practice factor (dimensionless, ranges from 0 to 1), which refers to the ratio of soil loss under specific conservation measures to soil loss under sloping tillage without conservation measures. Rainfall erosivity describes the likelihood of soil erosion caused by heavy rains; it is determined by rainfall kinetic energy and intensity (Wischmeier and Smith, 1978). In this study, we used the improved daily rainfall erosivity model to calculate the regional rainfall erosion factor by Zhang et al. (2002). Then, we adopted the Kriging method for interpolation to obtain the grating layer of the R factor. Soil erodibility describes the sensitivity of soil particles to water separation and movement; it also reflects soil sensitivity to erosion (Zhang et al., 2008). We used the revised erosion/productivity impact calculator formula for soil properties in China to calculate K (Zhang et al., 2008). The topographic factor reflects the effects of slope length and gradient on soil erosion (Wischmeier and Smith, 1978). The calculation of the topographic and C factors follows those in previous studies (Rao et al., 2014). Given the lack of relevant data, the effect of the P factor is not considered and can be set as 1 (Rao et al., 2014).

d xy ⎞ imxy = 1 − ⎛ ⎝ d mmax ⎠

2.3.3. WY WY is selected as an indicator of the hydrological regulation of ESs on the basis of the balance between precipitation and evapotranspiration. Its calculation formula is as follows:

2.3.5. Identification of ESMA The spatial cartographic method is used to spatialize the evaluation results of ESs and compare the spatial distribution patterns of ESs based on the overlay analysis and map algebra of the geographic information system (GIS) platform and to better express the spatial differentiation characteristics of trade-offs/synergies among ESs (Häyhä et al., 2015; Egoh et al., 2009). Although spatial overlay analysis does not quantify the complex interactions among services, it is an effective method for identifying multifunctional ES hotspots. In the current study, we proposed an identification framework for ES supply areas (ESSAs) from a temporal perspective and for ESMA from a spatiotemporal perspective (Fig. 2). For each ES, if an area can stably provide a certain level of this ES for a certain period, then this area can be regarded as an ESSA of this ES. Meanwhile, a region that can stably provide three or more ESs can be considered ESMA, which is the most effective and appropriate area for multiple ES protection under the constraints of financial and material resources. In accordance with previous studies (Egoh et al., 2008; Wu et al., 2013), we reclassified the intensity values of each ES in 2000, 2005, 2010, and 2015 and defined the top 50% of the grid values in the region as service supply areas for each ES. Then, we superimposed these areas from 2000 to 2015 to identify ESSA for each ES. That is, this area can stably supply this ES over time. Then, we overlaid the raster maps of different ESSAs. When a grid belongs to three or more ESSAs, it is defined as ESMA (Fig. 2).

WY = PPT − ET ± ΔS = PPT − ET

⎜

⎟

(7)

where m is the threat factor, M is the total threat factor, y represents a single raster in the threat factor m raster map, Ym represents the set of grids cells on m’s raster map, andωm represents the normalized threat weight. m y is used to determine whether grid y is the source of the threat factor m, and imxy is the distance function between the habitat and the threat source and the impact of the threat across space. βx represents the reachability level of the threat source to the grid x in the social, legal, and other protection states. Sjm represents the sensitivity of the land cover type j to the threat factor m, and d xy is the linear distance between grid cells x and y. d mmax is the maximum effective distance of threat m’s reach across space. Then, the HQ (HQxj ) of land cover type j at x is z

Dxj ⎞ ⎤ ⎡ HQxj = Hj ⎢1 − ⎜⎛ z ⎟ D + KZ ⎠⎥ xj ⎝ ⎦ ⎣

(8)

where Hj is the habitat suitability of land type j; and k, which is typically 0.5, represents scaling parameters. In addition, the value of HQxj is between 0 and 1. The parameters in the HQ calculation were listed in Supplementary A.

(5)

where WY is the water yield (mm), PPT is the precipitation (mm), and ET is the evapotranspiration (mm). The formula is based on the assumption that water storage change (ΔS) is negligible at a regional scale and long timescale. ET (mm) was derived from the MOD16 data product. The downloaded MOD16 original Hierarchical Data Format (HDF) files were converted into GeoTIFF format files using MODIS reprojection tool software. Operations, such as mosaicking, resampling, and conversion projection, were performed to obtain annual evapotranspiration data. Given that this water equilibrium model for calculating WY is disregarded at a regional scale and for long-term water storage, it differs slightly from the actual figures. However, previous studies have confirmed the usefulness of this model at a regional scale (Jia et al., 2014). 2.3.4. HQ A habitat comprises resources or environmental conditions occupied by species or population and the space for its development. HQ is the ability to provide development conditions, and quality is determined by the richness of regional natural resources (Hall et al., 1997). The HQ module in the InVEST model links land use/cover maps to threat sources and then assesses habitat distribution and degradation under different landscape patterns based on the response of various habitats to threat sources. The calculated HQ can reflect the biodiversity of the region, and the impact of land use pattern change on HQ can be observed. The measurement of HQ consists of four factors: the relative impact of each threat, the relative sensitivity of each habitat type to each threat on the landscape, the distance between habitat and threat sources, and the legal protection level from the disturbance in each cell. The higher the sensitivity of the habitat to the threat factor, the higher the impact of the factor on habitat degradation. Assuming that the land or habitat type at grid x is j, the degree of habitat degradation (Dxj ) at that point can be expressed as

2.3.6. Quantification of trade-offs and synergies Correlation analysis primarily examines whether a mutual dependence occurs between ESs by analyzing the correlation coefficient between two ESs (Willemen et al., 2010). To identify trade-offs/synergies between ESs, we adopted the Spearman rank correlation coefficient using R software to calculate the correlations among NPP, SC, WY, and HQ. If the coefficient is negative, then a trade-off relationship occurs between services, thereby indicating that the strength of one service increases while the other service becomes weak. Conversely, if the coefficient is positive, then potential synergy occurs between two services. The statistical parameter RMSE represents the sample standard deviation of the difference between the predicted and observed values; it is used as an indicator for quantifying the trade-off between two or 4

Ecological Indicators 107 (2019) 105566

L. Liu, et al.

Fig. 2. Identification framework for ESSA and ESMA.

value; and ESmin and ESmax are the minimum and maximum values for this type of ES across the entire region, respectively.

more ESs (Bradford and D'Amato, 2012; Feng et al, 2017; Lu et al., 2014; Yang et al., 2018). That is, RMSE extends the meaning of tradeoffs, not only to characterize negative correlations but also to include imbalances in the same direction of change between different ESs. This approach is a simple but effective means to represent the trade-off degree between two or more ESs. Furthermore, we quantify not only trade-offs but also synergies in this framework. To eliminate differences among ES units, we standardized the value of ESs within the range of 0–1, thereby effectively avoiding the inconsistency of comparative analysis due to different ES units. The standardized ES is calculated using the following equation:

ESi = (ESobs − ESmin )/(ESmax − ESmin )

RMSE =

1 × n−1

n

−

∑ (ESi − ES )2 i=1

(10)

−

where ES represents the expected value of the i number of ESs. RMSE represents the average difference between the standard value of a single ES and the average ES standard value in a given region, which describes the magnitude of the scattering of the mean. RMSE is the distance between the actual point and the diagonal (synergy: x + y−1 = 0 (Fig. 3a); trade-off: x−y = 0 (Fig. 3b)). The average of individual benefits is used to estimate the overall trade-offs and synergies of multiple ESs (Qin et al., 2015). In the synergy between two ESs, the

(9)

where ESi is the standard value for each ES; ESobs is the observation 5

Ecological Indicators 107 (2019) 105566

L. Liu, et al.

Fig. 3. Trade-offs (a) and synergies (b) between two ESs.

service changes considerably with time. By superimposing ESSAs of various ESs, the distribution map of ESMA in Shaanxi Province is shown in Fig. 5. Shaanxi Province’s ESMA is largely situated in Southern Shaanxi, particularly around the Micang Mountains and the Daba Mountains, and the south face of the Qinling Mountains. A small number of distributions also occur in the Huanglong Mountains in Northern Shaanxi. Among these regions, the area that can provide HQ, WY, and SC services comprises 10.84% of the total area of ESMA and is mostly distributed in the Micang Mountains. The areas that can provide NPP, WY, and SC services occupy 12.26% of the total area of ESMA and is largely distributed in the north of the Daba Mountains and the south of the Micang Mountains. The area that provides NPP, HQ, and SC services accounts for 71.01% of the total area of ESMA, which is the largest of all the service combinations, and is largely distributed in the south of the Qinling Mountains and around the Huanglong Mountains. The area that provides NPP, HQ, and WY services comprises 0.65% of the total area of ESMA, which is the smallest of all the service combinations, and its distribution is scattered. Furthermore, 5.24% of the ESMA area can provide all four ESs, and this area is mostly distributed in the southernmost part of Shaanxi.

RMSE value appears as an overall synergistic effect and increases with distance from the inverse 1:1 line. Simultaneously, the trade-offs between two ESs take the RMSE of each ES and increase in distances of 1:1. 3. Results 3.1. Spatiotemporal distribution characteristics of four ESs in Shaanxi Province The spatial patterns of the four ESs varied considerably across the entire region, and their spatial heterogeneity is conspicuous (Fig. 4). With the exception of WY, ESs also exhibit apparent spatial differentiation in areas bounded by the Qinling Mountains due to terrain and climate differences between the north and south of these mountains. High-value NPP is distributed in Hanzhong City, Ankang City, and Shangluo City in Southern Shaanxi. Areas with high SC capacity are mostly concentrated in the three aforementioned cities in Southern Shaanxi, but are also widely distributed in Yan’an City in Northern Shaanxi. Moreover, the heterogeneity of soil retention in spatial distribution is also the largest among the four services. The total amount of SC in Southern Shaanxi is considerably larger than those in Northern Shaanxi and Guanzhong. Although WY is affected by precipitation and interannual variations are evident, most WY remains concentrated and stably distributed in the northern part of Hanzhong City in Southern Shaanxi. The high-value areas of HQ are located in Southern Shaanxi and Yan’an areas in Northern Shaanxi. In terms of the temporal pattern, in addition to WY, ESs in 2000, 2005, 2010, and 2015 exhibit less interannual variations and tend to be stable. Given the development of the governance process of the Loess Plateau and the contribution of ecological restoration projects, the production status of vegetation has improved since 2000, and the annual NPP has increased every year. NPP in Guanzhong and Southern Shaanxi increased with a small fluctuation from 2000 to 2015. SC in each region presented an increasing trend during the 15 years. HQ exhibited a slight decreasing trend across the entire region. WY was considerably affected by the amount of precipitation in different years and demonstrated a fluctuation during the 15 years.

3.3. Correlation between ESs within and outside ESMA

3.2. Identification of ESMA

3.3.1. Differences of ESs within and outside ESMA The standard values of ESs within and outside ESMA are provided in Fig. 6. The values of NPP, HQ, and SC in ESMA are apparently higher than those in the entire region and non-ESMA. This finding is consistent with the result that the area can provide NPP, HQ, and SC services in ESMA, which is the largest area in all the service combinations. From 2000 to 2015, the NPP and HQ values in ESMA always maintain high consistency, and the aggregation area of high-value NPP is consistent with the HQ high-value aggregation area. This finding disagrees with the fact that the NPP and HQ values in the entire region exhibit floating and spatial heterogeneity. The SC value is considerably higher in ESMA than in the entire area and non-ESMA, and it has remained relatively stable from 2000 to 2015. In contrast with the three other services, the difference of WY in ESMA, the entire area, and non-ESMA is inevident. The WY value in ESMA is even slightly lower than that in the entire area in 2000 and 2010. Overall, the interannual variations of ESs in ESMA are smaller and more stable than those in the entire region.

ESSA for each ES is shown in Fig. 4. ESSA of NPP accounted for 41.47% of the total area of Shaanxi Province. This result indicates a small-scale spatial variation of the NPP supply area in Shaanxi Province from 2000 to 2015, and thus, this service is relatively stable. ESSAs of HQ and SC accounted for 36.36% and 32.76% of the total area, respectively. ESSA of WY is the smallest, thereby indicating that this

3.3.2. Quantification of trade-offs and synergies within and outside ESMA All the Spearman rank correlation coefficients between ES pairs passed the significance tests, thereby proving that a trade-off or synergy occurs between the six ES pairs. Simultaneously, we calculated the RMSE index between standardized ES pairs and obtained the intensity of trade-offs/synergies between ESs within and outside ESMA in 2000, 6

Ecological Indicators 107 (2019) 105566

L. Liu, et al.

Fig. 4. Spatiotemporal distribution of ESs in Shaanxi Province from 2000 to 2015.

7

Ecological Indicators 107 (2019) 105566

L. Liu, et al.

2005, 2010, and 2015 (Fig. 7). Within the boundary of ESMA, the relationships between ES pairs remain stable, except for the NPP and SC pair. The trade-off coefficients between SC and HQ in 2000, 2005, 2010, and 2015 are –0.58, –0.53, –0.48, and –0.49, respectively, with weak fluctuations and high tradeoff levels. The trade-off coefficients between WY and HQ from 2005 to 2015 exhibit a certain downward fluctuation with high trade-off levels. The relationship between NPP and WY also present trade-off from 2005 to 2015, but with slightly lower coefficients than the two other pairs. For the SC and WY and the NPP and HQ pairs, the RMSE values demonstrate slight fluctuations with negative values, thereby indicating a steady synergistic relationship between the two ES pairs. The variability of deviations of RMSE during the 15 years suggests that the trade-offs/ synergies among different ES pairs have been relatively stable within the boundary of ESMA. Outside the boundary of ESMA, the relationships between the other ES pairs have changed from 2000 to 2015, except for that of the NPP and SC pair, which has maintained a stable synergy. In general, the relationships between different ES pairs are highly unstable in non-ESMA; some years exhibit trade-offs, whereas other years are synergistic. In ESMA, most of the relationships between ES pairs are highly stable. SC and HQ, NPP and WY, and WY and HQ demonstrate high trade-off levels from 2000 to 2015 and remain stable over time. Meanwhile, the trade-off levels of NPP and WY and SC and WY remain low, but these pairs have highly steady synergy levels from 2000 to 2015. By contrast, the relationship between paired ESs in nonESMA has always changed over time, thereby showing instability. 3.3.3. Distribution of land use types within and outside ESMA To further determine whether multiple high-value ES aggregations exhibit certain relationships with land use/land cover (LULC) types, we calculated the area proportions of each LULC type within and outside ESMA and across the entire region (Table 1). The statistical results show that the primary land use types in ESMA are forest and grassland, which account for 44.87% and 40.83% of the total area, followed by cultivated land, water body, and other lands, which account for 14.16%, 0.09%, and 0.05%, respectively. In contrast with multifunctional areas, cultivated land and grassland account for 41.04% and 36.43% of the total land area in non-ESMA, followed by forest, other lands, and water body. The proportions of forest and grassland in ESMA are considerably higher than those in non-ESMA, and the proportion of cultivated land is considerably lower than that in non-ESMA. Existing studies have supported the view that forest land and grassland contribute the highest ES values. Costanza et al. (1997) estimated that the forest ES value accounted for 38% of the total terrestrial ES values. In terms of the ES value provided per unit area, forests are second only to wetlands and rivers, providing $969/ha y−1. Grasslands provide an ES value of $232/ ha y−1. In China, Xie et al. (2015) found that forest ES provides the highest proportion of the total ES, accounting for 46%. The proportion of grassland ES value accounts for 19.68% in China. Our results also show that forest and grassland are the primary land use types that can provide a variety of high-value ESs at the same pixel scale, which is consistent with the results of previous studies.

Fig. 5. Hotspot distribution of ES provision in Shaanxi Province.

4. Discussion 4.1. Variations in trade-offs/synergies within and outside EMSA The overall estimation of trade-offs/synergies between ESs is an extremely difficult and complex process (Bennett et al., 2009), and the interrelationships of various ESs in different regions may be varied changes or even reversals. In contrast with existing studies that focused on the trade-offs/synergies of ESs at spatial scales (Hou et al., 2017; Felipe-lucia et al., 2014), we explore the relationships of ESs within and outside ESMA and distinguish the similarities and differences in the trade-offs/systemic relationships of ESs in various areas with varying supply strengths.

Fig. 6. Rose map of ES standard values within and outside ESMA.

8

Ecological Indicators 107 (2019) 105566

L. Liu, et al.

Fig. 7. RMSE values of paired ESs within (a) and outside (b) ESMA Note: Positive numbers indicate synergy; negative numbers indicate trade-offs.

that afforestation can increase carbon sequestration and improve HQ, NPP and HQ exhibit a synergistic relationship, which is consistent with Yang’ research results. The relationships between WY and NPP, WY and HQ, and SC and HQ maintain relatively stable trade-off levels in ESMA from 2000 to 2015, whereas the relationships remain unstable in nonESMA. A possible reason for such finding is the fact that afforestation will increase carbon sequestration but also increase regional evapotranspiration, thereby decreasing WY. As an important part of the western development strategy, Shaanxi Province has begun implementing the Grain-for-Green Program in 1999 (Jia et al., 2014), and the forest coverage rate in the entire province increased from 24.1% in 2000 to 43.06% in 2015. Overall, trade-offs/synergies exhibit higher stability in ESMA than in non-ESMA. The ES supply level within ESMA tends to be stable over times. The trade-offs/synergies between ES pairs within ESMA, where ES provision is stable over time, can better illustrate the relationship between each ES pair, with higher instructive value for regional management.

Table 1 Distribution of land use types within and outside ESMA. Land use type

Cultivated land Forest Grassland Water body Other lands

Areas in ESMA (km2)

7747 24,545 22,331 48 27

Proportion (%) Within ESMA

Outside ESMA

Across the entire region

14.16 44.87 40.83 0.09 0.05

41.04 15.29 36.43 1.20 3.11

33.89 23.15 37.60 0.91 2.30

The largest difference within and outside ESMA is the relationship between NPP and SC. Some studies have shown that a strong synergistic relationship is aggregated between NPP and SC in the Loess Plateau region of Northern Shaanxi (Jia et al., 2014) and the Guanzhong–Tianshui region of Shaanxi (Qin et al., 2019), whereas strong trade-offs spatially occur in Southern Shaanxi (Li et al., 2017), thereby indicating that NPP and SC can exhibit opposing relationships in different regions. Previous studies have indicated that the relationship between NPP and SC may be affected by multiple factors (Hao et al., 2017). ES interactions may demonstrate certain relationships with LULC types. The relationship between SC and NPP tends to be trade-offs in the forest area (Li et al., 2017). As mentioned earlier, ESMA is mostly located in mountain areas with plenty of forests; thus, it is prone to exhibiting trade-offs. In addition, areas with high NPP typically have high vegetation coverage; hence, they can prevent soil erosion caused by precipitation, which produces a synergistic relationship between SC and NPP. However, precipitation is generally abundant where NPP value is high, which increases soil erosion, thereby producing a trade-off relationship between NPP and SC (Hao et al., 2017). Non-ESMA is mostly concentrated in the northern area of Shaanxi, and precipitation is considerably less than that in the south with a stable synergy state. By contrast, ESMA is concentrated in the southern mountains, and the relationship between NPP and SC is unstable. Similarly, WY and SC exhibit stable but relatively low synergy in ESMA. The dominant factor that influences WY and SC is precipitation, which leads to the synergistic relationship between these ESs. Wang et al. (2017) reported that the relationship between WY and SC in the western area upstream of the Hanjiang River is mostly synergistic, where the annual precipitation is approximately 700 mm, which nearly equals the precipitation in Southern Shaanxi. NPP and HQ maintain a stable and low level of synergy in ESMA, but the relationship between them is inconsistent in non-ESMA. Given

4.2. Comparison of multifunctional landscape hotspots, land use multifunctional hotspots, and ESMA Landscape multifunctionality and land use multifunctionality are based on the versatility of the landscape and land use, thereby representing multiple function provision in one landscape or piece of land (Table 2). Landscape multifunctionality is commonly understood as the coexistence of different landscape spheres, and it refers to landscape not only playing its primary ecological functions, but also other functions, such as social, economic, cultural, historical, and aesthetic functions, along with the interaction of different functions (Soini, 2001). Multifunctionality of land use refers to the state and performance of the land use function and the economic, environmental, and societal issues of a region. Multifunctionality is an important concept with regard to the relationship between land use and its functions (Pérez-Soba et al., 2008; Zhen et al., 2010); it typically uses land use functions to represent products and services provided by different land use patterns. At present, the identification of multifunctional landscape is mostly based on the overlapping of multiple landscape functions (services) to identify cold spots and hotspots (Willemen et al., 2010; Chan et al., 2006; Peng et al., 2016), mostly through multifunction evaluation, which is achieved by superimposing multiple functions to obtain a comprehensive index. In identifying land use multifunctional hotspots, most methods are directly based on the spatial differentiation of each land function value and on GIS spatial overlay analysis. Multifunctional land 9

10

The multifunctionality of land use refers to the state and performance of the land use function and the economic, environmental, and societal issues in a region (Pérez-Soba et al., 2008; Zhen et al., 2010). The areas of land use multifunctionality refer to hotspot areas with multiple high land function values (Dijst et al., 2005)

ESMA illustrates the versatility of services provided by ecosystems. It refers to a region that can stably provide multiple ESs

ESMA

Landscape multifunctionality is commonly understood as the coexistence of different landscape spheres; multifunctional landscape areas refer to hotspot areas that provide various landscape functions (Soini, 2001)

Multifunctional landscape hotspots

Land use multifunctional hotspots

Concept

Concept

Spatial

Spatial and temporal

ESs

Spatial

Landscape functions

Land functions

Scale

Research subject

The essence of multifunctional ESs is that the services provided by ecosystems are diverse, and ESs are based on the ecological processes of a certain time and space

The essence of land use multifunctionality is the shift from the traditional analysis of land function traits to the benefit of human society (Vereijken, 2003). The fundamental purpose is to achieve sustainable land use (Kates et al., 2001)

Landscape multifunctionality is an intrinsic trait that characterizes landscape (Cassatella and Seardo, 2014). The coordination and optimization of various landscape functions are basic goals in landscape management (Fürst et al., 2010; Carey, 2003)

Essential characteristics

Table 2 Comparison of multifunctional landscape hotspots, land use multifunctional hotspots, and ESMA.

The International Conference on Multifunctional Landscape, held in Roskiller, Denmark in 2000, first proposed the topic of Multifunctional Landscape Research (Brandt, 2007). Recently, the research on multifunctional landscape mainly focuses on conceptual definition, multi-functional evaluation, formation mechanism and spatial identification, planning, design and management and other fields (Peng et al., 2019) The concept of land use versatility is rooted in the concepts of agriculture, ecosystem products and services, and landscape functions. Its research mainly focuses on the definition of the concept and connotation of multi-functionality, the identification and classification of multi-functionality, and the evaluation of multi-functionality of land use (Liu et al., 2016) Research on multiple ESs focuses on hot spot identification, cluster of ESs, and the interaction of multiple ESs, as well as the bundles of multiple ESs (RaudseppHearne et al., 2010; Bennett et al., 2009; Gamfeldt et al., 2013) Landscape multifunctionality is more macroscopic than land versatility (Schößer et al., 2010). The material basis of landscape versatility is the diversity of ESs.

ESs are the basic elements and foundations for realizing landscape, land use, and ES multifunctionalities (Kienast et al., 2009). A multifunctional landscape provides a comprehensive representation of the interactions among ESs at the landscape service level (Turner et al., 2014). Land use multifunctionality is the comprehensive result of the integration of ESs and human activities in supporting land (Wiggering et al., 2006)

The concept of land use multifunctionality is based on land and landscape functions combined with the concept of versatility (Pérez-Soba et al., 2008; Paracchini et al., 2011)

Research progresses

Relationship

L. Liu, et al.

Ecological Indicators 107 (2019) 105566

Ecological Indicators 107 (2019) 105566

L. Liu, et al.

hotspots with multiple high land function values can be identified. The multifunctionality of landscape and land use exhibits evident spatial and temporal heterogeneities. To date, multifunctional landscapes or land use multifunctional hotspots are frequently identified by combining the individual ES maps of an area (Verhagen et al., 2016), thereby disregarding ES changes across time. From the temporal perspective, a long-term stable supply area of multiple ESs in a region, namely ESMA, should be the essential priority area to guarantee ES provision capacity across time. The identification and protection of ESMA are important paths toward prioritizing protected areas and optimizing protection efficiency via spatially explicit land management. 4.3. Limitations The trade-offs/synergies between ESs exhibit distinct and complex scale effects (Feng et al., 2017). Existing studies have suggested that the scale causes the trade-off relationship among ESs to change at different time and space scales (Qiao et al., 2019). Gordon and Enfors’ (2008) research showed that food production and soil conservation are synergistic at the watershed scale. However, this pair of services presents a trade-off relationship across the continent (Maes et al., 2012). In the current study, we consider the trade-offs of ESs at the provincial scale. The comparison of multiscale trade-offs may be further investigated in our subsequent work. In addition, management interested in maintaining locations of high ES supply should consider landscape configuration to effectively identify priority areas (Verhagen et al., 2016). In this study, we observed and compared the differences among various land uses, but did not account for landscape configuration. In fact, clearly non-uniform responses of ESs to the configuration effect are found, particularly at the cell scale. A detailed research on landscape configuration and multiple ES provision is proposed. Lastly, humans exert wide and profound impacts on ecosystems. Therefore, the protection and management of ESs should not only be confined to the delineation of ESMA but should further consider the interference of human activities on ES generation, particularly in population agglomeration areas. In a future research, human activities can be superimposed on the ESMA scope to explore the relationships between various high-value ES supplies and human activities, thereby assessing the impact of human activities and forecasting the risks that threaten ESMA. In order to verify the delimitation accuracy of the ESMA, the location of 22 National Nature Reserves in Shaanxi Province was collected and overlaid with the ESMA (Fig. 8). It can be found that almost all the National Nature Reserves in Shaanxi Province are located within the ESMA region, which indicates that ESMA is represented as critical zones for nature, providing shelters for wildlife habitats with a high level of multiple key ESs. ESMA can support to the delimitation of priority protected areas for optimizing the spatial layout planning of regional ecological protection.

Fig. 8. Overlay of ESMA and National Nature Reserves in Shaanxi in 2015.

changes in 2000, 2005, 2010 and 2015. Shaanxi Province’s ESMA is mostly situated in Southern, particularly around the mountain areas. The proportion of areas that provide high values of NPP, HQ, and SC is the highest among all service combinations. In addition, the major land use types distributed in ESMA are forest land and grassland, and the interannual changes in the trade-offs/synergies between ESs are more stable in ESMA than in non-ESMA. Our results can illustrate the relationships between ES pairs, which offer suggestions for optimized protection of ESs via spatially explicit land management, as well as provide scientifically basis on mapping approaches that have direct policy and decision-making relevance for land use planning. Acknowledgements This research was financially supported by Fund Project of Shaanxi Key Laboratory of Land Consolidation, China (2018-JC14), National Natural Science Foundation of China, China (41501087), Fund of Key Laboratory of Surveying and Mapping Science and Geospatial Information Technology of Ministry of Natural Resources, China (201904), and Chinese Universities Scientific Fund, China (2019TC117 and 2652017100).

5. Conclusions Exploring the temporal and spatial distribution of ESs is the basis and prerequisite for sustainable ecosystem management. From 2000 to 2015, the distribution of ESs exhibit heterogeneous in space in Shaanxi Province. Except for WY, each ES is bounded by the Qinling Mountains, and the other ESs tend to be stable and demonstrate less interannual

11

Ecological Indicators 107 (2019) 105566

L. Liu, et al.

Appendix Appendix A. Parameters for habitat quality in Shaanxi, China Threats

The properties of threats

Sensitivity of different land cover types

Weight Decay Maximum distance of influence Cultivated land Forest Grassland Water body Built-up area Other lands

Habitat

L_crp

L_urb

L_railway

L_highway

– – – 0.6 1 0.8 1 0 0

0.7 Linear 8 0.3 0.8 0.4 0.7 0 0

1 Exponential 10 0.5 1 0.6 0.9 0 0

0.8 Linear 8 0.4 0.6 0.7 0.6 0 0

1 Linear 10 0.4 0.7 0.5 0.5 0 0

of multiple ecosystem services are found in forests with more tree species. Nat. Commun. 4, 1340. Gordon, L.J., Enfors, E.I., 2008. Land degradation, ecosystem services, and resilience of smallholder farmers in Makanya catchment, Tanzania. Conserving Land Protecting Water 6, 33. Hall, L.S., Krausman, P.R., Morrison, M.L., 1997. The habitat concept and a plea for standard terminology. Wildlife Soc. B. 25 (1), 173–182. Hamann, M., Biggs, R., Reyers, B., 2015. Mapping social-ecological systems: identifying ‘green-loop’ and ‘red-loop’ dynamics based on characteristic bundles of ecosystem service use. Glob. Environ. Change 34, 218–226. Han, Z., Song, W., Deng, X.Z., Xu, X.Z., 2017. Trade-offs and synergies in ecosystem service within the three-rivers headwater region, China. Water. 9 (8), 588. Hao, R.F., Yu, D.Y., Sun, Y., Shi, M.C., 2019. The features and influential factors of interactions among ecosystem services. Ecol. Ind. 101, 770–779. Hao, R.F., Yu, D.Y., Wu, J.G., 2017. Relationship between paired ecosystem services in the grassland and agro-pastoral transitional zone of China using the constraint line method. Agr. Ecosyst. Environ. 240, 171–181. Häyhä, T., Franzese, P.P., Paletto, A., Fath, B.D., 2015. Assessing, valuing, and mapping ecosystem services in Alpine forests. Ecosyst. Serv. 14, 12–23. Hein, L., van Koppen, C.K., van Ierland, E.C., Leidekker, J., 2016. Temporal scales, ecosystem dynamics, stakeholders and the valuation of ecosystems services. Ecosyst. Serv. 21, 109–119. Hou, K., Lu, X., Zhang, J., 2015. GIS analysis of changes in ecological vulnerability using a SPCA model in the Loess Plateau of Northern Shaanxi. China. Int. J. Environ. Res. Public. Health 12 (4), 4292–4305. Hou, Y., Lu, Y.H., Chen, W.P., Fu, B.J., 2017. Temporal variation and spatial scale dependency of ecosystem service interactions: a case study on the central Loess Plateau of China. Landscape Ecol. 32 (6), 1201–1217. IPBES, 2019, Summary for policymakers of the global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. https://www.ipbes.net/system/tdf/spm_global_ unedited_advance.pdf?file=1&type=node&id=35245. Jaligot, R., Chenal, J., Bosch, M., Hasler, S., 2019. Historical dynamics of ecosystem services and land management policies in Switzerland. Ecol. Ind. 101, 81–90. Jia, X.Q., Fu, B.J., Feng, X.M., Hou, G.H., Liu, Y., Wang, X.F., 2014. The tradeoff and synergy between ecosystem services in the Grain-for-Green areas in Northern Shaanxi. China. Ecol. Indic. 43, 103–113. Kateb, H.E., Zhang, H., Zhang, P., Mosandl, R., 2013. Soil erosion and surface runoff on different vegetation covers and slope gradients: a field experiment in Southern Shaanxi Province, China. Catena 105, 1–10. Kates, R.W., Clark, W.C., Corell, R., Hall, J.M., Jaeger, C.C., Lowe, I., McCarthy, J.J., Schellnhuber, H.J., Bolin, B., Dickson, N.M., Faucheux, S., Gallopin, G.C., Grübler, A., Huntley, B., Jäger, J., Jodha, N.S., Kasperson, R.E., Mabogunje, A., Matson, P., Mooney, H., Moore III, B., O'Riordan, T., Svedin, U., 2001. Environment and Development: sustainability science. Science 292, 641–642. Kienast, F., Bolliger, J., Potschin, M., De Groot, R.S., Verburg, P.H., Heller, I., Wascher, D., Haines-Young, R., 2009. Assessing landscape functions with broad-scale environmental data: insights gained from a prototype development for Europe. Environ. Manage. 44 (6), 1099–1120. Li, Y.J., Zhang, L.W., Qiu, J.X., Yan, J.P., Wan, L.W., Wang, P.T., Hu, N.K., Cheng, W., Fu, B.J., 2017. Spatially explicit quantification of the interactions among ecosystem services. Landscape Ecol. 32 (6), 1181–1199. Liu, C., Xu, Y.Q., Sun, P.L., Liu, J., 2016. Progress and prospects of multi-functionality of land use research. Prog. Geog. 35 (9), 1087–1099 (in Chinese). Liu, L.Y., Tang, H., Caccetta, P., Lehmann, E.A., Hu, Y., Wu, X.L., 2013. Mapping afforestation and deforestation from 1974 to 2012 using Landsat time-series stacks in Yulin District, a key region of the Three-North Shelter region. China. Environ. Monit. Assess. 185 (12), 9949–9965. Lopes, R., Videira, N., 2017. Modelling feedback processes underpinning management of ecosystem services: the role of participatory systems mapping. Ecosyst. Serv. 28, 28–42. Lu, N., Fu, B.J., Jin, T.T., Chang, R.Y., 2014. Trade-off analyses of multiple ecosystem services by plantations along a precipitation gradient across Loess Plateau landscapes. Landscape Ecol. 29 (10), 1697–1708.

References Bennett, E.M., Peterson, G.D., Gordon, L.J., 2009. Understanding relationships among multiple ecosystem services. Ecol. Lett. 12 (12), 1394–1404. Bradford, J.B., D'Amato, A.W., 2012. Recognizing trade-offs in multi-objective land management. Front. Ecol. Environ. 10 (4), 210–216. Carey, A.B., 2003. Restoration of landscape function: reserves or active management? Forestry 76 (2), 221–230. Carpenter, S.R., Mooney, H.A., Agard, J., Capistrano, D., DeFries, R.S., Díaz, S., Dietz, T., Duraiappah, A.K., Oteng-Yeboah, A., Pereira, H.M., Perrings, C., Reid, W.V., Sarukhan, J., Scholes, R.J., Whyte, A., 2009. Science for managing ecosystem services: beyond The Millennium Ecosystem Assessment. PNAS 106 (5), 1305–1312. Cassatella, C., Seardo, B.M., 2014. In Search for Multifunctionality: The Contribution of Scenic Landscape Assessment. Landscape Planning and Rural Development. Springer, Cham, pp. 41–60. Chan, K.M., Shaw, M.R., Cameron, D.R., Underwood, E.C., Daily, G.C., 2006. Conservation planning for ecosystem services. PLoS Biol. 4 (11), e379. 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 (6630), 253. Dai, E.F., Wang, X.L., Zhu, J.J., Zhao, D.S., 2016. Methods, tools and research framework of ecosystem service trade-offs. Geogr. Res. 35 (6), 1005–1016 (in Chinese). Daily, G.C., 1997. Nature's Service. Island Press, Washington, DC. De Groot, R.S., Wilson, M.A., Boumans, R.M.J., 2002. A typology for the classification, description and valuation of ecosystem functions, goods and services. Ecol. Econ. 41 (3), 393–408. Dijst, M., Elbersen, B., Willis, K., 2005. The challenge of multi-functional land use in rural areas. J. Environ. Plann. Man. 48 (1), 3–6. Dobbs, C., Nitschke, C.R., Kendal, D., 2014. Global drivers and tradeoffs of three urban vegetation ecosystem services. PLoS One 9 (11), e113000. Egoh, B., Reyers, B., Rouget, M., Bode, M., Richardson, D.M., 2009. Spatial congruence between biodiversity and ecosystem services in South Africa. Biol. Conserv. 142 (3), 553–562. Egoh, B., Reyers, B., Rouget, M., Richardson, D.M., Le Maitre, D.C., van Jaarsveld, A.S., 2008. Mapping ecosystem services for planning and management. Agr. Ecosyst. Environ. 127 (1–2), 135–140. Felipe-Lucia, M.R., Comín, F.A., Bennett, E.M., 2014. Interactions among ecosystem services across land uses in a floodplain agroecosystem. Ecol. Soc. 19 (1). Feng, Q., Zhao, W.W., Fu, B.J., Ding, J.Y., Wang, S., 2017. Ecosystem service trade-offs and their influencing factors: a case study in the Loess Plateau of China. Sci. Total Environ. 607, 1250–1263. Field, C.B., Behrenfeld, M.J., Randerson, J.T., Falkowski, P., 1998. Primary production of the biosphere: integrating terrestrial and oceanic components. Science 281 (5374), 237–240. Field, C.B., Randerson, J.T., Malmström, C.M., 1995. Global net primary production: combining ecology and remote sensing. Remote Sens. Environ. 51 (1), 74–88. Fisher, B., Turner, R.K., Morling, P., 2009. Defining and classifying ecosystem services for decision making. Ecol. Econ. 68 (3), 643–653. Foley, J.A., DeFries, R., Asner, G.P., Barford, C., Bonan, G., Carpenter, S.R., Chapin, F.S., Coe, M.T., Daily, G.C., Gibbs, H.K., 2005. Global consequences of land use. Science 309 (5734), 570–574. Forman, R.T., 1995. Some general principles of landscape and regional ecology. Landscape Ecol. 10 (3), 133–142. Fu, B.J., Chen, L.D., Ma, K.M., Zhou, H.F., Wang, J., 2000. The relationships between land use and soil conditions in the hilly area of the loess plateau in northern Shaanxi, China. Catena 39 (1), 69–78. Fürst, C., Volk, M., Pietzsch, K., Makeschin, F., 2010. Pimp your landscape: a tool for qualitative evaluation of the effects of regional planning measures on ecosystem services. Environ. Manage. 46 (6), 953–968. Gamfeldt, L., Snäll, T., Bagchi, R., Jonsson, M., Gustafsson, L., Kjellander, P., Ruiz-Jaen, M.C., Fröberg, M., Stendahl, J., Philipson, C.D., Mikusiński, G., Andersson, E., Westerlund, B., Andrén, H., Moberg, F., Moen, G., Bengtsson, J., 2013. Higher levels

12

Ecological Indicators 107 (2019) 105566

L. Liu, et al.

mapping and map-making. Landscape Urban Plan. 57 (3–4), 225–239. Swinton, S.M., Lupi, F., Robertson, G.P., Hamilton, S.K., 2007. Ecosystem services and agriculture: cultivating agricultural ecosystems for diverse benefits. Ecol. Econ. 64 (2), 245–252. Tilman, D., 2000. Causes, consequences and ethics of biodiversity. Nature 405 (6783), 208–211. Turner, K.G., Odgaard, M.V., Bøcher, P.K., Dalgaard, T., Svenning, J.C., 2014. Bundling ecosystem services in Denmark: trade-offs and synergies in a cultural landscape. Landscape Urban Plan. 125, 89–104. Vejre, H., Abildtrup, J., Andersen, E., Andersen, P.S., Brandt, J., Busck, A., Dalgaard, T., Hasler, B., Huusom, H., Kristensen, L.S., Kristensen, S.P., Præstholm, S., 2007. Multifunctional agriculture and multifunctional landscapes-land use as an interface. Multifunctional land use. Springer, Berlin, Heidelberg, pp. 93–104. Vereijken, P.H., 2003. Transition to multifunctional land use and agriculture. NJASWagen. J. life Sc. 50 (2), 171–179. Verhagen, W., Van Teeffelen, A.J.A., Compagnucci, A.B., Poggio, L., Gimona, A., Verburg, P.H., 2016. Effects of landscape configuration on mapping ecosystem service capacity: a review of evidence and a case study in Scotland. Landscape Ecol. 31 (7), 1457–1479. Vihervaara, P., Kumpula, T., Tanskanen, A., Burkhard, B., 2010. Ecosystem services–A tool for sustainable management of human–environment systems. Case study Finnish Forest Lapland. Ecol. Complex. 7 (3), 410–420. Wang, P.T., Zhang, L.W., Li, Y.J., Jiao, L., Wang, H., Yan, J.P., Lv, Y.H., Fu, B.J., 2017. Spatio-temporal characteristics of the trade-off and synergy relationships among multiple ecosystem services in the upper reaches of Hanjiang River Basin. Acta Geogr. Sin. 72 (11), 131–142 (in Chinese). Wiggering, H., Dalchow, C., Glemnitz, M., Helming, K., Müller, K., Schultz, A., Stachow, U., Zander, P., 2006. Indicators for multifunctional land use-linking socio-economic requirements with landscape potentials. Ecol. Ind. 6 (1), 238–249. Willemen, L., Hein, L., van Mensvoort, M.E.F., Verburg, P.H., 2010. Space for people, plants, and livestock? Quantifying interactions among multiple landscape functions in a Dutch rural region. Ecol. Ind. 10 (1), 62–73. Wischmeier, W.H., Smith, D.D., 1978. Predicting rainfall erosion losses-a guide to conservation planning. Predicting rainfall erosion losses-a guide to conservation planning. Wu, J.S., Feng, Z., Gao, Y., Peng, J., 2013. Hotspot and relationship identification in multiple landscape services: a case study on an area with intensive human activities. Ecol. Ind. 29, 529–537. Xie, G.D., Zhang, C.X., Zhang, C.S., Xiao, Y., Lu, C.X., 2015. The value of ecosystem services in China. Resour. Sci. 37 (9), 1740–1746 (in Chinese). Yang, S.Q., Zhao, W.W., Liu, Y.X., Wang, S., Wang, J., Zhai, R.J., 2018. Influence of land use change on the ecosystem service trade-offs in the ecological restoration area: dynamics and scenarios in the Yanhe watershed. China. Sci. Total Environ. 644, 556–566. Zhang, K.L., Shu, A.P., Xu, X.L., Yang, Q.K., Yu, B., 2008. Soil erodibility and its estimation for agricultural soils in China. J. Arid Environ. 72 (6), 1002–1011. Zhang, W.B., Xie, Y., Liu, B.Y., 2002. Rainfall erosivity estimation using daily rainfall amounts. Sci. Geogr. Sin. 6, 705–711 (in Chinese). Zhen, L., Wei, Y.J., Xie, G.D., Katharina, H., Cao, S.Y., Yang, L., Pan, Y., Hannes, K., 2010. Regional analysis of dynamic land use functions in China. Acta Ecol. Sin. 30 (24), 6749–6761 (in Chinese). Zhu, W.Q., Pan, Y.Z., Yang, X.M., Song, G.B., 2007. Comprehensive analysis of the impact of climatic changes on Chinese terrestrial net primary productivity. Chinese Sci. Bull. 52 (23), 3253–3260.

Maes, J., Paracchini, M.L., Zulian, G., Dunbar, M.B., Alkemade, R., 2012. Synergies and trade-offs between ecosystem service supply, biodiversity, and habitat conservation status in Europe. Biol. Conserv. 155, 1–12. Mastrangelo, M.E., Weyland, F., Villarino, S.H., Barral, M.P., Nahuelhual, L., Laterra, P., 2014. Concepts and methods for landscape multifunctionality and a unifying framework based on ecosystem services. Landscape Ecol. 29 (2), 345–358. Odum, E.P., Barrett, G.W., 1971. Fundamentals of Ecology. Saunders, Philadelphia. Palomo, I., Martín-López, B., Potschin, M., Haines-Young, R., Montes, C., 2013. National Parks, buffer zones and surrounding lands: mapping ecosystem service flows. Ecosyst. Serv. 4, 104–116. Paracchini, M.L., Pacini, C., Jones, M.L.M., Pérez-Soba, M., 2011. An aggregation framework to link indicators associated with multifunctional land use to the stakeholder evaluation of policy options. Ecol. Ind. 11 (1), 71–80. Peng, J., Chen, X., Liu, Y.X., Lü, H.L., Hu, X.X., 2016. Spatial identification of multifunctional landscapes and associated influencing factors in the Beijing-Tianjin-Hebei region. China. Appl. Geogr. 74, 170–181. Peng, J., Hu, X.X., Qiu, S.J., Hu, Y.N., Meersmans, J., Liu, Y.X., 2019. Multifunctional landscapes identification and associated development zoning in mountainous area. Sci. Total Environ. 660, 765–775. Pérez-Soba, M., Petit, S., Jones, L., Bertrand, N., Briquel, V., Omodei-Zorini, L., Contini, C., Helming, K., Farrington, J.H., Mossello, M.T., Wascher, D., Kienast, F., De Groot, R., 2008. Land use Functions: A Multifunctionality Approach to Assess the Impact of Land use Changes on Land use Sustainability. Sustainability impact assessment of land use changes. Springer, Berlin, Heidelberg, pp. 375–404. Potter, C.S., Randerson, J.T., Field, C.B., Matson, P.A., Vitousek, P.M., Mooney, H.A., Klooster, S.A., 1993. Terrestrial ecosystem production: a process model based on global satellite and surface data. Global Biogeochem. Cyst. 7 (4), 811–841. Qiao, X.N., Gu, Y.Y., Zou, C.Y., Xu, D.L., Wang, L., Ye, X., Yang, Y., Huang, X.F., 2019. Temporal variation and spatial scale dependency of the trade-offs and synergies among multiple ecosystem services in the Taihu Lake Basin of China. Sci. Total Environ. 651, 218–229. Qin, K.Y., Li, J., Liu, J.Y., Yan, L.W., Huang, H., J, 2019. Setting conservation priorities based on ecosystem services-A case study of the Guanzhong-Tianshui Economic Region. Sci. Total Environ. 650, 3062–3074. Qin, K.Y., Li, J., Yang, X.N., 2015. Trade-off and synergy among ecosystem services in the Guanzhong-Tianshui Economic Region of China. Int. J. Environ. Res. Pub. He. 12 (11), 14094–14113. Qiu, J.X., Turner, M.G., 2013. Spatial interactions among ecosystem services in an urbanizing agricultural watershed. P. Natl. Acad. Sci. 110 (29), 12149–12154. Rau, A.L., Von Wehrden, H., Abson, D.J., 2018. Temporal dynamics of ecosystem services. Ecol. Econ. 151, 122–130. Rao, E.M., Ouyang, Z.Y., Yu, X.X., Xiao, Y., 2014. Spatial patterns and impacts of soil conservation service in China. Geomorphology 207, 64–70. Raudsepp-Hearne, C., Peterson, G.D., Bennett, E.M., 2010. Ecosystem service bundles for analyzing tradeoffs in diverse landscapes. PNAS 107 (11), 5242–5247. Raynolds, M.K., Comiso, J.C., Walker, D.A., Verbyla, D., 2008. Relationship between satellite-derived land surface temperatures, arctic vegetation types, and NDVI. Remote Sens. Environ. 112 (4), 1884–1894. Rodriguez, J.P., Beard Jr, T.D., Bennett, E.M., Cumming, G.S., Cork, S.J., Agard, J., Dobson, A.P., Peterson, G.D., 2006. Trade-offs across space, time, and ecosystem services. Ecol. Soc. 11 (1), 28. Schößer, B., Helming, K., Wiggering, H., 2010. Assessing land use change impacts–a comparison of the sensor land use function approach with other frameworks. J. Land Use Sci. 5 (2), 159–178. Soini, K., 2001. Exploring human dimensions of multifunctional landscapes through

13