# Land Use Zoning for Conserving Ecosystem Services under the Impact of Climate Change: A Case Study in the Middle Reaches of the Heihe River Basin

Authors: Chenchen Shi; Jinyan Zhan; Yongwei Yuan; Feng Wu; Zhihui Li
Journal: Advances in Meteorology
 (2015)

---

## Abstract

Ecosystem services are the benefit human populations derive directly and indirectly from the natural environment. They suffer from both the human intervention, like land use zoning change, and natural intervention, like the climate change. Under the background of climate change, regulation services of ecosystem could be strengthened under proper land use zoning policy to mitigate the climate change. In this paper, a case study was conducted in the middle reaches of the Heihe River Basin to assess the ecosystem services conservation zoning under the change of land use associated with climate variations. The research results show the spatial impact of land use zoning on ecosystem services in the study area which are significant reference for the spatial optimization of land use zoning in preserving the key ecosystem services to mitigate the climate change. The research contributes to the growing literature in finely characterizing the ecosystem services zones altered by land use change to alleviate the impact of climate change, as there is no such systematic ecosystem zoning method before.


---

## Body

## 1. Introduction

Ecosystem services are the natural environment provision and utility the ecosystem provides upon which human-beings survival lie [1]. It is the material basis to maintain the human survival and support earth life system [2]. According to the Millennium Ecosystem Assessment (MA) [3], the ecosystem services can be mainly divided into four categories, namely, provisioning, supporting, regulating, and cultural services. The degree of human demand influences the sustaining of the ecosystem services [4]. Meanwhile, the services provided by the ecosystem directly affect the human well-being [5, 6]. Both natural and human interventions will exert great influence on the ecosystem services, which will hinder the social and economic development and have an impact on the human well-being.Among those natural interventions of ecosystem services, climate change is likely to affect the abundance, production, distribution, and quality of ecosystems [7]. Therefore, ecosystem services, such as climate stabilization through carbon sequestration, the provision of nonirrigated forage for livestock and wildlife species, the delivery of water which supports fish for commercial and recreational sport fishing, the provision of critical habitat for biodiversity, and many other types of ecosystem services, are likely to be impacted by a changing climate [8]. And in return the ecosystem service conservation and restoration under proper land use zoning policy play a key role in climate mitigation [9].As for human interferences that influence the ecosystem services, numerous studies have showed that land use change associated with climate variations could affect the structure and function of ecosystem and then affect the supply of ecosystem services functions [2, 10, 11]. The spatial land use allocation may affect the provision of a wide range of ecosystem services and is instrumental to their conservation and enhancement [12].Previous studies have shown that, in order to characterize the ecological changes, monitoring landscape pattern changes provides an indirect approach, as such changes would influence a variety of ecological processes and functions [13–15]. However, ecosystem services zoning, which aims at identifying major regional ecological environment problems and therefore establishes ecosystem services zones, seems to be a more direct way of showing the ecosystem services change associated with land use change.Ecological function zoning is one of the hot issues among academia in the field of ecology [16]. There are various researches of ecosystem service function zoning from different research aspects [17–19]. However, these researches either focus on only one type of ecosystem [19] or are confined to land use zoning [20]. The spatial changes of the regional ecosystem service have not been reflected. And case studies of grid scale space dynamic zoning are absent.A scientifically sound guidance is always needed for planners to reveal the mechanism of the ecological and climatic effect of land use change and additionally to promote sustainable development. In order to achieve that, it is necessary to conduct zoning characterization of the ecosystem service under the comprehensive consideration of land use and land cover change and the site condition. While there is lack of a successful model and systematic method to accomplish this goal in the academia, this paper addressed the issue by presenting a case study research aimed at empirically exploring how land use zoning will affect the ecosystem services. Different land use zoning through the time period from 1988 to 2008 was analyzed and the ecosystem services were identified and zoned by using SIZES. The SIZES can automatically generate the boundaries of the ecosystem service function zone in the case study area after the input of the raster data and therefore visually reveals the temporal and spatial variation characteristics of the regional ecosystem services. Finally, through the spatial and temporal variation of the key ecosystem services generated by SIZES, effects of the land use zoning on selected ecosystem services are then qualitatively analyzed.The study area is located in the Heihe River Basin (HRB), China. As one of the largest inland river basins in Northeast China, it is characterized by the multiple natural landscapes and is composed by various natural geographic units. As a typical arid area, it comprises almost all the major natural characters. And the living of the residents or farmers there depends largely upon the provision of the ecosystem services. While in the middle reaches, there are 95% of the cultivated land, 91% of the population and over 80% of the production value of the whole watershed. It is the main agricultural production area of the whole region. However, in recent decades, with the strengthening of human economic activities, a series of ecological environment problems appear, such as grassland degradation, natural forest decrease, species number and glacier area reduction, land desertification and secondary salinization, and water scarcity. These problems are all associated with the land use and land cover change and the quality changes of ecological environment. The degraded ecosystem will weaken the ecosystem services function, thus influencing the social-economic development of the whole region. Therefore we selected the three cities and five counties along the middle reaches as the study area to assess the impact of LUC on ecosystem services. And the method we used will lend some credence to the related research in similar areas.The following content of this paper is divided into four parts, namely, methodology, case study, results, and discussion and conclusions. In the methodology part, we will explicitly introduce SIZES. Thereafter, a case study of the middle reaches of the HRB is presented to elaborate the land use zoning for conserving ecosystem services under the impact of climate change in this ecological fragile region. Major results of this study are stated in the “Results” section. Finally several discussions are made to conclude this research.


## 2. Methodology

We built a system for identifying and zoning ecosystem services (SIZES) to identify the key ecosystem services and their spatial and temporal variation. The system is composed of the following four principal components: indicator system for ecosystem services evaluation based on analytic hierarchy process (ISESE-AHP), ranking matrix of ecosystem services based on factor analysis (RMES-FA), identification of the core ecosystem services based on principal component evaluation (ICES-PCE), and identification and zoning of ecosystem services based on unsupervised fuzzy clustering (IDES-UFC).


### 2.1. Principle of SIZES

SIZES delimitates the functional zones for regional ecosystem services based on core ecosystem services in the minimum consistent units. The indicator system of ecosystem services in SIZES is developed with the conceptual framework provided by MEA (2003) and the influencing factors at different spatiotemporal scales are accordingly selected. The core ecosystem services, the indicators which comprehensively represent the most important ecosystem services in a functional zone, are identified with factor analysis according to their contribution to the variance of indicators of ecosystem services. The core ecosystem services are interpreted according to the factor loading matrix obtained with the principal component method and further calculated on the 1 km × 1 km grid cells. Finally, the functional zones for ecosystem services are delimitated with fuzzy clustering analysis and critical areas for ecosystem services are identified. The procedures of SIZES are illustrated more specifically as the following.(
1) AHP-Based Indicator System for Ecosystem Services Assessment (ISESA-AHP). An assessment hierarchy is constructed to break down ecosystem services into a series of interrelated factors. Since the regional ecosystem services are influenced by various factors [8], measurable primary indicators for ecosystem services have been selected based on the analytic hierarchy process so as to integrate these various influencing factors into one model to determine their interdependencies and the perceived consequences interactively. With a conceptual framework provided by MEA [21], important influencing factors of ecosystem services are elicited according to the driving mechanism of regional ecosystem services. There are four levels in this assessment hierarchy, that is, the overall goal at the top, criteria and subcriteria in the middle, and the measurable primary indicators at the bottom. During the process of selecting primary indicators of regional ecosystem services, three principles are followed: (1) the indicators characterizing ecosystem services are easy to access; (2) the selected indicators should cover all the aspects of regional ecosystem services; and (3) the selected indicators should be coherent with the upper-layer indicators.(
2) FA-Based Extraction of the Core Ecosystem Services (ECES-FA). Factor analysis is used to summarize and aggregate the aforementioned primary indicators of ecosystem services into a few principal components (PCs) that represent certain core ecosystem services. It describes the variance of primary indicators of ecosystem services in terms of potentially fewer comprehensive indicators, that is, PCs. The PCs for a data set are defined as linear combinations of variables that account for the maximum variance within the entire data set and are orthogonal to each other. With the correlation matrix of standardized data of primary indicators of ecosystem services, the eigenvalues, that is, the contribution of core ecosystem services to the variance of primary indicators of ecosystem services, are calculated and their corresponding eigenvectors are obtained. The contribution rate of each PC and its cumulative contribution rates are also calculated. It is assumed that the PCs with higher contribution rate represent better core ecosystem services. Only the PCs with the cumulative contribution rate reaching a certain threshold are retained. The threshold is generally set to be above 60%. Therefore, only the PCs with the cumulative contribution reaching above 60% are retained and used to represent the core ecosystem services in SIZES.(
3) PCM-Based Ranking Matrix of Ecosystem Services (RMES-PCM). The relative importance of primary indicators to the core ecosystem services is analyzed with RMES-PCM, a ranking matrix based on the factor loading matrix obtained with principal component method (PCM). The factor loading reflects the correlation between one primary indicator and one certain PC. The core ecosystem services are identified and interpreted with the primary indicators with the heavy factor loading; a correlation matrix involving these primary indicators was used as an input for the analysis instead of a covariance matrix, resulting in normalized PCM. Since there may be no representative indicators of the PCs, it is necessary to conduct factor rotation so as to obtain satisfactory PCs. The eigenvalues are used as a criterion for selecting PCs for the factor rotation. Only PCs with eigenvalues above 1 are retained and subjected to varimax rotation for maximizing the correlation between primary indicators and PCs [22]. The RMES-PCM is obtained by sorting the primary indicators of ecosystem services according to their factor loadings in descending order. The relative importance of primary indicators of ecosystem services to PCs is ranked to represent certain core ecosystem services.With RMES-PCM, the core ecosystem services can be identified and interpreted according to the most important primary indicators of ecosystem services. The factor scores of the retained PCs are calculated at each 1 km × 1 km grid cell to spatially represent the core ecosystem services. The factor score can be calculated with many methods, such as the regression estimation method, Bartlett estimation method, and Thomson estimation method. SIZES adopts Bartlett estimation method and spatially represents the core ecosystem services on 1 km × 1 km grid cells in the study area with the factor score.(
4) UFC-Based Zonation of Ecosystem Services (ZES-UFC). With the core ecosystem services represented at 1 km × 1 km grid cells, the functional zones for ecosystem services are delimitated with unsupervised fuzzy clustering (UFC) algorithm, which is based on the modified fuzzy c-means algorithm and maximum likelihood estimation (MLE). The minimum consistent units (samples) are coherent in the type of core ecosystem services and they are classified into clusters with unsupervised optimal fuzzy clustering algorithm. The core ecosystem services of samples within one cluster are more similar to each other than that belonging to different clusters, where the similarity is defined by a distance measure. This partition approach defines the cluster to which each sample belongs, using the elements of the membership matrix. Each sample is attached with different degrees of membership to each cluster. The initial centers of each partition are found with the fuzzy c-means with the Euclidean distance and the fuzzy c-means with the exponential (or MLE) distance are used for precise modeling of the Gaussian mixture in the samples. The number of clusters is determined by finding the partition with minimum number of Gaussian clusters. Specifically, the number of clusters is determined by beginning with one group at the center of data and adding another group with each iteration process at a location in which membership in the previous groups is small, while examining the partition validity criteria until reaching the maximal number of groups. Finally the number of clusters providing maximal validity criteria measurements is chosen for optimal partition.


### 2.2. Function Modules Embedded in SIZES

The procedures mentioned above are realized with four function modules in SIZES, that is, MIES (module for identifying ecosystem services), MRES (module for ranking ecosystem services), MECES (module for identifying the core ecosystem services), and MZES (module for zoning ecosystem services).


#### 2.2.1. Module for Identifying Ecosystem Services (MIES)

MIES mainly consists of an indicator system for ecosystem services, representing interaction processes between ecosystem services and influencing factors, and enabling interactive access to SIZES. There are four indicators that fit broadly at the top level in this indicator system, namely, supporting, provisioning, regulating, and cultural services. Since the indicators at lower levels may vary among different researches, users of SIZES may recode the program and set these indicators by specifying input parameters and the system will dynamically track the input indicators. Users may also access the online help system, clarifying the underlying assumptions and formal definitions of parameters used in this software. Within MIES of SIZES, each cell (1 km × 1 km) is assigned with multilayer values to represent attributes of each factor (e.g., NPP and land quality). With the 1-kilometer (km) grid pixel data model, MIES can combine, analyze, and interpret the multiple and hierarchical spatial factors efficiently in consideration of maintenance and improvement of ecosystem services.


#### 2.2.2. Module for Extracting the Core Ecosystem Services (MECES)

MECES aims to extract the PCs that represent the core ecosystem services according to established criteria with factor analysis. The core ecosystem services are extracted with the following steps.(1) Assessment matrix for ecosystem services is calculated with(1)
X
=
x
i
j
n
×
p
=
x
11
x
12
⋯
x
1
p
x
21
x
22
⋯
x
2
p
⋮
⋮
⋯
⋮
x
n
1
x
n
2
⋯
x
n
p
,where x
i
j is the value of the jth (
j
=
1
,
2
,
…
,
p
) variable of ecosystem services on the ith (
i
=
1
,
2
,
…
,
p
) grid cell. The vector on the ith (
i
=
1
,
2
,
…
,
n
) grid cell can be presented by X
i
=
(
x
i
1
,
x
i
2
,
…
,
x
i
p
)
′.
(2) Sample deviation matrix is constructed with(2)
E
=
∑
i
=
1
n
X
i
-
X
-
X
i
-
X
-
′
≜
e
k
l
,where X
-
=
1
/
n
∑
i
=
1
n
X
i
=
(
x
-
1
,
x
-
1
,
…
,
x
-
1
)
′.
(3) Sample correlation matrix is calculated as(3)
R
=
e
k
l
e
k
k
e
l
l
,
k
,
l
=
1
,
2
,
…
,
p
.
(4) The eigenvalues and corresponding eigenvectors are obtained with following steps.
The eigenvalueλ
i (i
=
1,2
,
…
,
p) of the equation (λ
I
-
R
=
0) is obtained with the Jacobi method and the values are ranked as λ
1
≥
λ
2
≥
⋯
≥
λ
p
≥
0. Then l
i
(
i
=
1,2
,
…
,
p
), the eigenvector of the corresponding eigenvalue λ
i, is calculated. Concurrently, the condition should agree with l
i
=
1; that is, ∑
j
=
1
p
l
i
j
2
=
1, where l
i
j is the jth of vector (
l
i
j
). Then the ith PC is Z
i
=
l
i
X.
(5) The eigenvalues represent the contributions of these PCs to the variance of primary indicators of ecosystem services. The contribution rate of PCs is calculated withm
i
=
λ
i
/
∑
j
=
1
p
λ
j. The PCs are sorted according to their contribution rate in descending order and then their cumulative contribution rate is calculated with ∑
k
=
1
i
m
i. When the cumulative contribution rate of the PCs satisfies ∑
k
=
1
i
m
i
≥
δ, where δ is set to be 60%, the first m PCs are selected to represent the core ecosystem services.


#### 2.2.3. Module for Ranking Ecosystem Services (MRES)

MRES determines the relative importance of primary indicators to the core ecosystem services with RMES-PCM. The interface of MRES embedded in SIZES is shown in Figure2. The factor loading matrix is obtained with the PCM. The factor loading vector is calculated with (4)
a
i
=
λ
i
l
i
i
=
1,2
,
…
,
p
,where λ
i is the eigenvalue corresponding to the ith retained PC and l
i is its corresponding orthonormal eigenvector. The factor loading matrix is A
=
(
a
1
,
…
,
a
p
). The varimax rotation is conducted for maximizing the correlation between primary indicators and PCs for making the representative variables of PCs more prominent. The RMES-PCM is then obtained by sorting the primary indicators of ecosystem services in the factor loading matrix according to their factor loadings in descending order. With RMES-PCM, the core ecosystem services can be identified and interpreted according to the primary indicators of ecosystem services with the factor loading. Finally, SIZES calculates the factor scores of core ecosystem services with Bartlett method and spatially represents the core ecosystem services on 1 km × 1 km grid cells with their factor scores.


#### 2.2.4. Module for Zoning Ecosystem Services (MZES)

MZES utilizes unsupervised fuzzy clustering (UFC) analysis to delimitate and generate functional zones for ecosystem services. The steps of delimitating functional zones for ecosystem services with unsupervised fuzzy clustering analysis are as follows.(1) Calculate the average value of all data records:(5)
P
=
1
m
∑
i
=
1
m
X
iand set the number of clusters (
c
) to be one; that is, c
=
1.
(2) LetU
k
i
-
o
l
d
=
0 (U
k
i
-
o
l
d
=
p
r
e
v
i
o
u
s
U
k
i).The basic fuzzy c-means algorithm (FCM).(3) Calculate the Euclidean distance (for alli
,
k):(6)
d
2
P
k
,
X
i
=
P
k
-
X
i
T
·
P
k
-
X
i
.
(4) Calculate the degree of membership (U
k
i) of all data records in all clusters (for all i
,
k):(7)
U
k
i
=
1
/
d
2
P
k
,
X
i
∑
j
=
1
c
1
/
d
2
P
j
,
X
i
.
(5) Compute the new set of cluster centers,P
1
→
k (for all k):(8)
P
k
=
∑
i
=
1
m
U
k
i
2
X
i
∑
i
=
1
m
U
k
i
2
.
(6) Ifm
a
x
k
,
i
U
k
i
-
U
k
i
-
o
l
d
>
ε, then U
k
i
-
o
l
d
=
U
k
i: go to Step  3, or else a loop is started and continues Step  6.The modified fuzzy c-means algorithm based on MLE (FMLE).(7) Compute the MLE distance.The sum of memberships within thekth cluster (the a priori probability of the kth cluster) is(9)
A
k
=
1
m
∑
i
=
1
m
U
k
i
.The fuzzy covariance matrix of the kth cluster is(10)
F
k
=
∑
i
=
1
m
U
k
i
P
k
-
X
i
P
k
-
X
i
T
∑
i
=
1
m
U
k
i
.The fuzzy hypervolume of the kth cluster is (11)
H
k
=
d
e
t
F
k
1
/
2
.Then the MLE distance (for all i
,
k) is (12)
d
2
P
k
,
X
i
=
H
k
A
k
e
x
p
P
k
-
X
i
T
·
F
k
-
1
·
P
k
-
X
i
2
.
(8) Repeat Steps  4 and  5 using the MLE distance.
(9) Ifm
a
x
k
,
i
U
k
i
-
U
k
i
-
o
l
d
>
ε, U
k
i
-
o
l
d
=
U
k
i, go to Step  7; if not, continue.
(10) Compute the partition and average density validity criteria for choosing the optional number of clusters. For instance, the partition density criterion definition is(13)
P
D
c
=
∑
j
=
1
c
A
j
∑
j
=
1
c
H
j
.
(11) Create an imaginary center that has the same distance (very far) from all data records, for instance, by choosing the initial distance to the new cluster as(14)
d
2
P
k
+
1
,
X
i
=
10
·
v
a
r
i
a
n
c
e
X
i
=
1
,
…
,
m
.
(12) Letc
=
c
+
1. If c
≤
c
m
a
x (c
m
a
x is the maximum feasible number of clusters in the data), U
k
i
-
o
l
d
=
0 and go to Step  4; or else, the UFC algorithm is completed.The matricesU and P created in Step  9 are the required output of the UFC algorithm.


## 2.1. Principle of SIZES

SIZES delimitates the functional zones for regional ecosystem services based on core ecosystem services in the minimum consistent units. The indicator system of ecosystem services in SIZES is developed with the conceptual framework provided by MEA (2003) and the influencing factors at different spatiotemporal scales are accordingly selected. The core ecosystem services, the indicators which comprehensively represent the most important ecosystem services in a functional zone, are identified with factor analysis according to their contribution to the variance of indicators of ecosystem services. The core ecosystem services are interpreted according to the factor loading matrix obtained with the principal component method and further calculated on the 1 km × 1 km grid cells. Finally, the functional zones for ecosystem services are delimitated with fuzzy clustering analysis and critical areas for ecosystem services are identified. The procedures of SIZES are illustrated more specifically as the following.(
1) AHP-Based Indicator System for Ecosystem Services Assessment (ISESA-AHP). An assessment hierarchy is constructed to break down ecosystem services into a series of interrelated factors. Since the regional ecosystem services are influenced by various factors [8], measurable primary indicators for ecosystem services have been selected based on the analytic hierarchy process so as to integrate these various influencing factors into one model to determine their interdependencies and the perceived consequences interactively. With a conceptual framework provided by MEA [21], important influencing factors of ecosystem services are elicited according to the driving mechanism of regional ecosystem services. There are four levels in this assessment hierarchy, that is, the overall goal at the top, criteria and subcriteria in the middle, and the measurable primary indicators at the bottom. During the process of selecting primary indicators of regional ecosystem services, three principles are followed: (1) the indicators characterizing ecosystem services are easy to access; (2) the selected indicators should cover all the aspects of regional ecosystem services; and (3) the selected indicators should be coherent with the upper-layer indicators.(
2) FA-Based Extraction of the Core Ecosystem Services (ECES-FA). Factor analysis is used to summarize and aggregate the aforementioned primary indicators of ecosystem services into a few principal components (PCs) that represent certain core ecosystem services. It describes the variance of primary indicators of ecosystem services in terms of potentially fewer comprehensive indicators, that is, PCs. The PCs for a data set are defined as linear combinations of variables that account for the maximum variance within the entire data set and are orthogonal to each other. With the correlation matrix of standardized data of primary indicators of ecosystem services, the eigenvalues, that is, the contribution of core ecosystem services to the variance of primary indicators of ecosystem services, are calculated and their corresponding eigenvectors are obtained. The contribution rate of each PC and its cumulative contribution rates are also calculated. It is assumed that the PCs with higher contribution rate represent better core ecosystem services. Only the PCs with the cumulative contribution rate reaching a certain threshold are retained. The threshold is generally set to be above 60%. Therefore, only the PCs with the cumulative contribution reaching above 60% are retained and used to represent the core ecosystem services in SIZES.(
3) PCM-Based Ranking Matrix of Ecosystem Services (RMES-PCM). The relative importance of primary indicators to the core ecosystem services is analyzed with RMES-PCM, a ranking matrix based on the factor loading matrix obtained with principal component method (PCM). The factor loading reflects the correlation between one primary indicator and one certain PC. The core ecosystem services are identified and interpreted with the primary indicators with the heavy factor loading; a correlation matrix involving these primary indicators was used as an input for the analysis instead of a covariance matrix, resulting in normalized PCM. Since there may be no representative indicators of the PCs, it is necessary to conduct factor rotation so as to obtain satisfactory PCs. The eigenvalues are used as a criterion for selecting PCs for the factor rotation. Only PCs with eigenvalues above 1 are retained and subjected to varimax rotation for maximizing the correlation between primary indicators and PCs [22]. The RMES-PCM is obtained by sorting the primary indicators of ecosystem services according to their factor loadings in descending order. The relative importance of primary indicators of ecosystem services to PCs is ranked to represent certain core ecosystem services.With RMES-PCM, the core ecosystem services can be identified and interpreted according to the most important primary indicators of ecosystem services. The factor scores of the retained PCs are calculated at each 1 km × 1 km grid cell to spatially represent the core ecosystem services. The factor score can be calculated with many methods, such as the regression estimation method, Bartlett estimation method, and Thomson estimation method. SIZES adopts Bartlett estimation method and spatially represents the core ecosystem services on 1 km × 1 km grid cells in the study area with the factor score.(
4) UFC-Based Zonation of Ecosystem Services (ZES-UFC). With the core ecosystem services represented at 1 km × 1 km grid cells, the functional zones for ecosystem services are delimitated with unsupervised fuzzy clustering (UFC) algorithm, which is based on the modified fuzzy c-means algorithm and maximum likelihood estimation (MLE). The minimum consistent units (samples) are coherent in the type of core ecosystem services and they are classified into clusters with unsupervised optimal fuzzy clustering algorithm. The core ecosystem services of samples within one cluster are more similar to each other than that belonging to different clusters, where the similarity is defined by a distance measure. This partition approach defines the cluster to which each sample belongs, using the elements of the membership matrix. Each sample is attached with different degrees of membership to each cluster. The initial centers of each partition are found with the fuzzy c-means with the Euclidean distance and the fuzzy c-means with the exponential (or MLE) distance are used for precise modeling of the Gaussian mixture in the samples. The number of clusters is determined by finding the partition with minimum number of Gaussian clusters. Specifically, the number of clusters is determined by beginning with one group at the center of data and adding another group with each iteration process at a location in which membership in the previous groups is small, while examining the partition validity criteria until reaching the maximal number of groups. Finally the number of clusters providing maximal validity criteria measurements is chosen for optimal partition.


## 2.2. Function Modules Embedded in SIZES

The procedures mentioned above are realized with four function modules in SIZES, that is, MIES (module for identifying ecosystem services), MRES (module for ranking ecosystem services), MECES (module for identifying the core ecosystem services), and MZES (module for zoning ecosystem services).


### 2.2.1. Module for Identifying Ecosystem Services (MIES)

MIES mainly consists of an indicator system for ecosystem services, representing interaction processes between ecosystem services and influencing factors, and enabling interactive access to SIZES. There are four indicators that fit broadly at the top level in this indicator system, namely, supporting, provisioning, regulating, and cultural services. Since the indicators at lower levels may vary among different researches, users of SIZES may recode the program and set these indicators by specifying input parameters and the system will dynamically track the input indicators. Users may also access the online help system, clarifying the underlying assumptions and formal definitions of parameters used in this software. Within MIES of SIZES, each cell (1 km × 1 km) is assigned with multilayer values to represent attributes of each factor (e.g., NPP and land quality). With the 1-kilometer (km) grid pixel data model, MIES can combine, analyze, and interpret the multiple and hierarchical spatial factors efficiently in consideration of maintenance and improvement of ecosystem services.


### 2.2.2. Module for Extracting the Core Ecosystem Services (MECES)

MECES aims to extract the PCs that represent the core ecosystem services according to established criteria with factor analysis. The core ecosystem services are extracted with the following steps.(1) Assessment matrix for ecosystem services is calculated with(1)
X
=
x
i
j
n
×
p
=
x
11
x
12
⋯
x
1
p
x
21
x
22
⋯
x
2
p
⋮
⋮
⋯
⋮
x
n
1
x
n
2
⋯
x
n
p
,where x
i
j is the value of the jth (
j
=
1
,
2
,
…
,
p
) variable of ecosystem services on the ith (
i
=
1
,
2
,
…
,
p
) grid cell. The vector on the ith (
i
=
1
,
2
,
…
,
n
) grid cell can be presented by X
i
=
(
x
i
1
,
x
i
2
,
…
,
x
i
p
)
′.
(2) Sample deviation matrix is constructed with(2)
E
=
∑
i
=
1
n
X
i
-
X
-
X
i
-
X
-
′
≜
e
k
l
,where X
-
=
1
/
n
∑
i
=
1
n
X
i
=
(
x
-
1
,
x
-
1
,
…
,
x
-
1
)
′.
(3) Sample correlation matrix is calculated as(3)
R
=
e
k
l
e
k
k
e
l
l
,
k
,
l
=
1
,
2
,
…
,
p
.
(4) The eigenvalues and corresponding eigenvectors are obtained with following steps.
The eigenvalueλ
i (i
=
1,2
,
…
,
p) of the equation (λ
I
-
R
=
0) is obtained with the Jacobi method and the values are ranked as λ
1
≥
λ
2
≥
⋯
≥
λ
p
≥
0. Then l
i
(
i
=
1,2
,
…
,
p
), the eigenvector of the corresponding eigenvalue λ
i, is calculated. Concurrently, the condition should agree with l
i
=
1; that is, ∑
j
=
1
p
l
i
j
2
=
1, where l
i
j is the jth of vector (
l
i
j
). Then the ith PC is Z
i
=
l
i
X.
(5) The eigenvalues represent the contributions of these PCs to the variance of primary indicators of ecosystem services. The contribution rate of PCs is calculated withm
i
=
λ
i
/
∑
j
=
1
p
λ
j. The PCs are sorted according to their contribution rate in descending order and then their cumulative contribution rate is calculated with ∑
k
=
1
i
m
i. When the cumulative contribution rate of the PCs satisfies ∑
k
=
1
i
m
i
≥
δ, where δ is set to be 60%, the first m PCs are selected to represent the core ecosystem services.


### 2.2.3. Module for Ranking Ecosystem Services (MRES)

MRES determines the relative importance of primary indicators to the core ecosystem services with RMES-PCM. The interface of MRES embedded in SIZES is shown in Figure2. The factor loading matrix is obtained with the PCM. The factor loading vector is calculated with (4)
a
i
=
λ
i
l
i
i
=
1,2
,
…
,
p
,where λ
i is the eigenvalue corresponding to the ith retained PC and l
i is its corresponding orthonormal eigenvector. The factor loading matrix is A
=
(
a
1
,
…
,
a
p
). The varimax rotation is conducted for maximizing the correlation between primary indicators and PCs for making the representative variables of PCs more prominent. The RMES-PCM is then obtained by sorting the primary indicators of ecosystem services in the factor loading matrix according to their factor loadings in descending order. With RMES-PCM, the core ecosystem services can be identified and interpreted according to the primary indicators of ecosystem services with the factor loading. Finally, SIZES calculates the factor scores of core ecosystem services with Bartlett method and spatially represents the core ecosystem services on 1 km × 1 km grid cells with their factor scores.


### 2.2.4. Module for Zoning Ecosystem Services (MZES)

MZES utilizes unsupervised fuzzy clustering (UFC) analysis to delimitate and generate functional zones for ecosystem services. The steps of delimitating functional zones for ecosystem services with unsupervised fuzzy clustering analysis are as follows.(1) Calculate the average value of all data records:(5)
P
=
1
m
∑
i
=
1
m
X
iand set the number of clusters (
c
) to be one; that is, c
=
1.
(2) LetU
k
i
-
o
l
d
=
0 (U
k
i
-
o
l
d
=
p
r
e
v
i
o
u
s
U
k
i).The basic fuzzy c-means algorithm (FCM).(3) Calculate the Euclidean distance (for alli
,
k):(6)
d
2
P
k
,
X
i
=
P
k
-
X
i
T
·
P
k
-
X
i
.
(4) Calculate the degree of membership (U
k
i) of all data records in all clusters (for all i
,
k):(7)
U
k
i
=
1
/
d
2
P
k
,
X
i
∑
j
=
1
c
1
/
d
2
P
j
,
X
i
.
(5) Compute the new set of cluster centers,P
1
→
k (for all k):(8)
P
k
=
∑
i
=
1
m
U
k
i
2
X
i
∑
i
=
1
m
U
k
i
2
.
(6) Ifm
a
x
k
,
i
U
k
i
-
U
k
i
-
o
l
d
>
ε, then U
k
i
-
o
l
d
=
U
k
i: go to Step  3, or else a loop is started and continues Step  6.The modified fuzzy c-means algorithm based on MLE (FMLE).(7) Compute the MLE distance.The sum of memberships within thekth cluster (the a priori probability of the kth cluster) is(9)
A
k
=
1
m
∑
i
=
1
m
U
k
i
.The fuzzy covariance matrix of the kth cluster is(10)
F
k
=
∑
i
=
1
m
U
k
i
P
k
-
X
i
P
k
-
X
i
T
∑
i
=
1
m
U
k
i
.The fuzzy hypervolume of the kth cluster is (11)
H
k
=
d
e
t
F
k
1
/
2
.Then the MLE distance (for all i
,
k) is (12)
d
2
P
k
,
X
i
=
H
k
A
k
e
x
p
P
k
-
X
i
T
·
F
k
-
1
·
P
k
-
X
i
2
.
(8) Repeat Steps  4 and  5 using the MLE distance.
(9) Ifm
a
x
k
,
i
U
k
i
-
U
k
i
-
o
l
d
>
ε, U
k
i
-
o
l
d
=
U
k
i, go to Step  7; if not, continue.
(10) Compute the partition and average density validity criteria for choosing the optional number of clusters. For instance, the partition density criterion definition is(13)
P
D
c
=
∑
j
=
1
c
A
j
∑
j
=
1
c
H
j
.
(11) Create an imaginary center that has the same distance (very far) from all data records, for instance, by choosing the initial distance to the new cluster as(14)
d
2
P
k
+
1
,
X
i
=
10
·
v
a
r
i
a
n
c
e
X
i
=
1
,
…
,
m
.
(12) Letc
=
c
+
1. If c
≤
c
m
a
x (c
m
a
x is the maximum feasible number of clusters in the data), U
k
i
-
o
l
d
=
0 and go to Step  4; or else, the UFC algorithm is completed.The matricesU and P created in Step  9 are the required output of the UFC algorithm.


## 2.2.1. Module for Identifying Ecosystem Services (MIES)

MIES mainly consists of an indicator system for ecosystem services, representing interaction processes between ecosystem services and influencing factors, and enabling interactive access to SIZES. There are four indicators that fit broadly at the top level in this indicator system, namely, supporting, provisioning, regulating, and cultural services. Since the indicators at lower levels may vary among different researches, users of SIZES may recode the program and set these indicators by specifying input parameters and the system will dynamically track the input indicators. Users may also access the online help system, clarifying the underlying assumptions and formal definitions of parameters used in this software. Within MIES of SIZES, each cell (1 km × 1 km) is assigned with multilayer values to represent attributes of each factor (e.g., NPP and land quality). With the 1-kilometer (km) grid pixel data model, MIES can combine, analyze, and interpret the multiple and hierarchical spatial factors efficiently in consideration of maintenance and improvement of ecosystem services.


## 2.2.2. Module for Extracting the Core Ecosystem Services (MECES)

MECES aims to extract the PCs that represent the core ecosystem services according to established criteria with factor analysis. The core ecosystem services are extracted with the following steps.(1) Assessment matrix for ecosystem services is calculated with(1)
X
=
x
i
j
n
×
p
=
x
11
x
12
⋯
x
1
p
x
21
x
22
⋯
x
2
p
⋮
⋮
⋯
⋮
x
n
1
x
n
2
⋯
x
n
p
,where x
i
j is the value of the jth (
j
=
1
,
2
,
…
,
p
) variable of ecosystem services on the ith (
i
=
1
,
2
,
…
,
p
) grid cell. The vector on the ith (
i
=
1
,
2
,
…
,
n
) grid cell can be presented by X
i
=
(
x
i
1
,
x
i
2
,
…
,
x
i
p
)
′.
(2) Sample deviation matrix is constructed with(2)
E
=
∑
i
=
1
n
X
i
-
X
-
X
i
-
X
-
′
≜
e
k
l
,where X
-
=
1
/
n
∑
i
=
1
n
X
i
=
(
x
-
1
,
x
-
1
,
…
,
x
-
1
)
′.
(3) Sample correlation matrix is calculated as(3)
R
=
e
k
l
e
k
k
e
l
l
,
k
,
l
=
1
,
2
,
…
,
p
.
(4) The eigenvalues and corresponding eigenvectors are obtained with following steps.
The eigenvalueλ
i (i
=
1,2
,
…
,
p) of the equation (λ
I
-
R
=
0) is obtained with the Jacobi method and the values are ranked as λ
1
≥
λ
2
≥
⋯
≥
λ
p
≥
0. Then l
i
(
i
=
1,2
,
…
,
p
), the eigenvector of the corresponding eigenvalue λ
i, is calculated. Concurrently, the condition should agree with l
i
=
1; that is, ∑
j
=
1
p
l
i
j
2
=
1, where l
i
j is the jth of vector (
l
i
j
). Then the ith PC is Z
i
=
l
i
X.
(5) The eigenvalues represent the contributions of these PCs to the variance of primary indicators of ecosystem services. The contribution rate of PCs is calculated withm
i
=
λ
i
/
∑
j
=
1
p
λ
j. The PCs are sorted according to their contribution rate in descending order and then their cumulative contribution rate is calculated with ∑
k
=
1
i
m
i. When the cumulative contribution rate of the PCs satisfies ∑
k
=
1
i
m
i
≥
δ, where δ is set to be 60%, the first m PCs are selected to represent the core ecosystem services.


## 2.2.3. Module for Ranking Ecosystem Services (MRES)

MRES determines the relative importance of primary indicators to the core ecosystem services with RMES-PCM. The interface of MRES embedded in SIZES is shown in Figure2. The factor loading matrix is obtained with the PCM. The factor loading vector is calculated with (4)
a
i
=
λ
i
l
i
i
=
1,2
,
…
,
p
,where λ
i is the eigenvalue corresponding to the ith retained PC and l
i is its corresponding orthonormal eigenvector. The factor loading matrix is A
=
(
a
1
,
…
,
a
p
). The varimax rotation is conducted for maximizing the correlation between primary indicators and PCs for making the representative variables of PCs more prominent. The RMES-PCM is then obtained by sorting the primary indicators of ecosystem services in the factor loading matrix according to their factor loadings in descending order. With RMES-PCM, the core ecosystem services can be identified and interpreted according to the primary indicators of ecosystem services with the factor loading. Finally, SIZES calculates the factor scores of core ecosystem services with Bartlett method and spatially represents the core ecosystem services on 1 km × 1 km grid cells with their factor scores.


## 2.2.4. Module for Zoning Ecosystem Services (MZES)

MZES utilizes unsupervised fuzzy clustering (UFC) analysis to delimitate and generate functional zones for ecosystem services. The steps of delimitating functional zones for ecosystem services with unsupervised fuzzy clustering analysis are as follows.(1) Calculate the average value of all data records:(5)
P
=
1
m
∑
i
=
1
m
X
iand set the number of clusters (
c
) to be one; that is, c
=
1.
(2) LetU
k
i
-
o
l
d
=
0 (U
k
i
-
o
l
d
=
p
r
e
v
i
o
u
s
U
k
i).The basic fuzzy c-means algorithm (FCM).(3) Calculate the Euclidean distance (for alli
,
k):(6)
d
2
P
k
,
X
i
=
P
k
-
X
i
T
·
P
k
-
X
i
.
(4) Calculate the degree of membership (U
k
i) of all data records in all clusters (for all i
,
k):(7)
U
k
i
=
1
/
d
2
P
k
,
X
i
∑
j
=
1
c
1
/
d
2
P
j
,
X
i
.
(5) Compute the new set of cluster centers,P
1
→
k (for all k):(8)
P
k
=
∑
i
=
1
m
U
k
i
2
X
i
∑
i
=
1
m
U
k
i
2
.
(6) Ifm
a
x
k
,
i
U
k
i
-
U
k
i
-
o
l
d
>
ε, then U
k
i
-
o
l
d
=
U
k
i: go to Step  3, or else a loop is started and continues Step  6.The modified fuzzy c-means algorithm based on MLE (FMLE).(7) Compute the MLE distance.The sum of memberships within thekth cluster (the a priori probability of the kth cluster) is(9)
A
k
=
1
m
∑
i
=
1
m
U
k
i
.The fuzzy covariance matrix of the kth cluster is(10)
F
k
=
∑
i
=
1
m
U
k
i
P
k
-
X
i
P
k
-
X
i
T
∑
i
=
1
m
U
k
i
.The fuzzy hypervolume of the kth cluster is (11)
H
k
=
d
e
t
F
k
1
/
2
.Then the MLE distance (for all i
,
k) is (12)
d
2
P
k
,
X
i
=
H
k
A
k
e
x
p
P
k
-
X
i
T
·
F
k
-
1
·
P
k
-
X
i
2
.
(8) Repeat Steps  4 and  5 using the MLE distance.
(9) Ifm
a
x
k
,
i
U
k
i
-
U
k
i
-
o
l
d
>
ε, U
k
i
-
o
l
d
=
U
k
i, go to Step  7; if not, continue.
(10) Compute the partition and average density validity criteria for choosing the optional number of clusters. For instance, the partition density criterion definition is(13)
P
D
c
=
∑
j
=
1
c
A
j
∑
j
=
1
c
H
j
.
(11) Create an imaginary center that has the same distance (very far) from all data records, for instance, by choosing the initial distance to the new cluster as(14)
d
2
P
k
+
1
,
X
i
=
10
·
v
a
r
i
a
n
c
e
X
i
=
1
,
…
,
m
.
(12) Letc
=
c
+
1. If c
≤
c
m
a
x (c
m
a
x is the maximum feasible number of clusters in the data), U
k
i
-
o
l
d
=
0 and go to Step  4; or else, the UFC algorithm is completed.The matricesU and P created in Step  9 are the required output of the UFC algorithm.


## 3. Case Study in the Heihe River Basin

### 3.1. Study Area

The Heihe River Basin, with an altitude between 2000 m and 5500 m, is one of the largest inland river basins in Western China, which stretches over 821 km and covers an area of 142.9 thousand km2. It is characterized by the multiple natural landscapes of ice, river, lake, oasis, and desert and is composed by the natural geographic units of high mountain snow-ice zone, plain oasis zone, and Gobi desert zone. As a typical arid area, it comprises almost all the major natural characters and therefore possesses various ecosystem services. And the living of the residents or farmers depends largely upon the provision of the ecosystem services. In the middle reaches, the runoff has been fully used, while 95% of the cultivated land, 91% of the population, and over 80% of the production value of the whole watershed gathered here. It is the main agricultural production area of the whole region. Both human and natural interventions will exert great influence on the ecosystem services, which will hinder the social and economic development and thus impact the human well-being.In this situation, climate change plays a key role in altering the provision service of the middle reaches and the region depends largely on the regulating services of ecosystem to mitigate this change. Therefore investigating how the land use zoning could conserve the ecosystem service and thus counter the impact of climate change is of great significance in this area. So we selected the three cities and five counties along the middle reaches as the study area (Figure1).Figure 1
Location and geographic boundaries of the middle reaches of the HRB.Figure 2
(a) Land use patterns of middle reaches of the HRB in 1988 and (b) land use patterns of middle reaches of the HRB in 2008.
(a)
(b)Identifying the general characteristics of the regional ecological environment is of great significance to the rational regional ecosystem service function identification. In the middle reaches, the key ecological environment problems are land desertification and salinization and water pollution. Besides, due to the lack of monitoring and control measures, the condition is getting worse. The reduction of grassland area caused by the massive land reclamation as well as overgrazing leads to the grassland degradation. The unreasonable industrial structure not only limited the regional economy development, but also worsened the supply and demand contradictions of the water resources, leading to further problems such as land desertification and salinization.


### 3.2. Data Selection

The selection of the ecosystem services is based on the characteristics of the study area mentioned above and is under the guidance of MA, which has been systematically elaborated in Section2.1, “(
1
) AHP-based indicator system for ecosystem services evaluation” part. In this study, we divided the ecosystem services in the middle reaches of HRB into three major categories, namely, provisioning services, regulating services, and supporting services. And in each category, there are subcategories and further items, each being relevant to the land use change process that characterizes this study; thus an index system of the key ecosystem services of the HRB was built (Tables 1 and 2). Given the fact that the cultural services function is difficult to quantify, this research did not take the cultural function into consideration.Table 1
Key ecosystem services indicators in the HRB.
Categories
Subcategories
Indicators and description
Provisioning  services
Food
Cultivated land percentage %
Fresh water
Water area percentage %
Fuelwood
Forestry area percentage %
Regulating services
Climate regulating
Air temperature: mean annual temperature (0.1°C)
Temperature departure
Precipitation departure
Cumulated temperature (>10°C) Annually cumulated temperature of daily meanAir temperature above 10°C (0.1°C)
Sunshine hours: sunshine hours rectifying the spatial variability of solar  radiation
Sunshine hours departure
Species diversity
Landscape diversity index (Shannon diversity index)
Supporting services
Soil formation
Loam proportion %
Slope: terrain slope derived from DEM (0.01 degrees)
Nutrient cycling
Soil phosphorus content %
Primary production
Net primary production gC m−2yr−1
Grassland area
Grassland percentage %Table 2
Descriptive statistics of key ecosystem services indicators in the Heihe River Basin.
Variable
Unit
Mean
Std. Dev.
Min
Max
Cultivated land
%
5.52
1873.59
0
100
Water area
%
4.49
1616.67
0
100
Grassland
%
16.79
3091.04
0
100
Forestry area
%
1.12
660.01
0
100
Air temperature
°C
0.51
2.54
0
0.96
Temperature departure
°C
0
2.52
−5.59
4.27
Precipitation departure
mm
0
75.04
−137.93
213.55
Cumulated temperature
°C
2479.10
733.29
461.29
3942.56
Sunshine hours
h
8.33
0.35
7.31
9.30
Sunshine hours departure
h
-
2.85
e
-
15
0.31
−1.00
0.76
Landscape diversity index
—
−37.81
116.46
−460.52
460.52
Loam proportion
%
22.50
6.37
4
34
Slope
0.01°
4.72
5.53
0
60.57
Soil phosphorus content
%
0.07
0.01
0.02
0.08
Net primary production
gC m−2yr−1
97.54
112.16
0
897
Note:  total observation: 161584.In this study, the information about land use change was derived from the remote sensing images of the HRB in the years of 1988, 1955, 2000, and 2008 through the human-computer interactive interpretation (Figure2). And the information can be used as the basic data for analyzing the different land use zoning.


### 3.3. Data Processing

Our study uses a land cover classification data set developed by the Chinese Academy of Sciences (CAS) [23, 24] as the background data for spatially identifying the ecosystem services. The data set is derived from Landsat TM/ETM images of bands 3, 4, and 5 with a spatial resolution of 30 m × 30 m. And the classification contains 6 types of land uses/covers, including cultivated land, forestry, grassland, water area, built-up area, and unused land [25].


#### 3.3.1. Raster Data Processing

The identification of the key ecosystem services function is based on the raster data model. The indicators of the ecosystem services function are characterized by multi-source and multi-scale. Therefore, the 1 km area percentage data [26, 27] (Figure 3) is an expressive way of data integration and fusion in monitoring, forecasting, and analyzing the regional ecosystem services changes. The procedure to generate the 1 km area percentage data set includes four steps. The first step is to generate a vector map of land cover and land use changes during the study period at a scale of 1 : 100,000 based on the remote sensing Landsat TM/ETM data. The second step is to generate a 1 km FISHNET vector map georeferenced to a boundary map of the study area at a scale of 1 : 100,000. Each cell of the generated 1 km FISHNET vector map owns a unique ID. The third step is to overlay the vector map of land cover and land use changes with the 1 km FISHNET vector map. This is done by aggregating converted areas in each 1 km grid identified by cell IDs of the 1 km FISHNET vector map in the TABLE module of ArcGIS software. Finally, the area percentage vector data are transformed into grid raster data to identify the conversion direction and intensity. The design of working flow insists on no loss of area information. Without special notification, the statistical area of cultivated land according to the GRID data is survey area by satellite remote sensing data, which can be called “gross area.”Figure 3
Flowchart of generating 1 km area percentage data set of land use changes.


#### 3.3.2. Indicator Data Processing

Among the indicators that are needed in the identification of the key ecosystem services (Figure4), the raster data of the net primary production are provided by the CAS. The data sources and processing methods of other indicators are listed below.(1) Food, water, and fuelwood supply: in order to precisely characterize the supply of grain, water, and fuelwood in each cell, in this research, we use the percentage of the cultivated land, water area, and forestry to indicate these provisions. The 1 km area percentage data processing method is shown above in Figure3.
(2) Climate variables: all the climatic variables are generated from the site-based observations from the China Meteorological Administration. The spline interpolation algorithm is employed to make the surface data of climatic variables acquired at observation stations [28, 29]. The values for the climatic variables during simulation period are estimated using the space-time stochastic model [30].
(3) Landscape diversity index: we selected the Shannon diversity index to measure the landscape diversity index [31].
The Shannon diversity index(
H
) is expressed as [32] (15)
H
=
∑
i
p
i
ln
⁡
p
i
,where p
i is the share of land use type i in total area of the landscape.
(4) Soil data [33]: soil data were from multiple sources. Firstly, the specific soil survey points were calibrated on local administrative division base map, based on carefully comparing the information on National Secondary Soil Survey data and 1 : 250000 topographic maps. Secondly, an Albert projection value for specific soil survey point was assigned onto the map and generated a soil survey base map. After that, soil attribute data were linked onto soil survey base map by using the Kring interpolation algorithm and generated spatial soil data encompassing soil type and physicochemical characters. Finally, these soil data were saved as ArcGIS Grid format, with soil type, soil depth, soil organic content, and soil texture information.
(5) Elevation data: the elevation data were derived from China 1 : 250000 digital elevation models, with a spatial resolution of 50 m. Average slope information was extracted on each 1 km × 1 km plot by GIS spatial analyzing function; then it was used as basic information for land suitability evaluation.Figure 4
Spatial dispersion and characterization of the ecosystem services indicators, (a) cultivated land percentage, (b) grassland area percentage, (c) air temperature, (d) precipitation, (e) soil phosphorus content, and (f) net primary production.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.1. Study Area

The Heihe River Basin, with an altitude between 2000 m and 5500 m, is one of the largest inland river basins in Western China, which stretches over 821 km and covers an area of 142.9 thousand km2. It is characterized by the multiple natural landscapes of ice, river, lake, oasis, and desert and is composed by the natural geographic units of high mountain snow-ice zone, plain oasis zone, and Gobi desert zone. As a typical arid area, it comprises almost all the major natural characters and therefore possesses various ecosystem services. And the living of the residents or farmers depends largely upon the provision of the ecosystem services. In the middle reaches, the runoff has been fully used, while 95% of the cultivated land, 91% of the population, and over 80% of the production value of the whole watershed gathered here. It is the main agricultural production area of the whole region. Both human and natural interventions will exert great influence on the ecosystem services, which will hinder the social and economic development and thus impact the human well-being.In this situation, climate change plays a key role in altering the provision service of the middle reaches and the region depends largely on the regulating services of ecosystem to mitigate this change. Therefore investigating how the land use zoning could conserve the ecosystem service and thus counter the impact of climate change is of great significance in this area. So we selected the three cities and five counties along the middle reaches as the study area (Figure1).Figure 1
Location and geographic boundaries of the middle reaches of the HRB.Figure 2
(a) Land use patterns of middle reaches of the HRB in 1988 and (b) land use patterns of middle reaches of the HRB in 2008.
(a)
(b)Identifying the general characteristics of the regional ecological environment is of great significance to the rational regional ecosystem service function identification. In the middle reaches, the key ecological environment problems are land desertification and salinization and water pollution. Besides, due to the lack of monitoring and control measures, the condition is getting worse. The reduction of grassland area caused by the massive land reclamation as well as overgrazing leads to the grassland degradation. The unreasonable industrial structure not only limited the regional economy development, but also worsened the supply and demand contradictions of the water resources, leading to further problems such as land desertification and salinization.


## 3.2. Data Selection

The selection of the ecosystem services is based on the characteristics of the study area mentioned above and is under the guidance of MA, which has been systematically elaborated in Section2.1, “(
1
) AHP-based indicator system for ecosystem services evaluation” part. In this study, we divided the ecosystem services in the middle reaches of HRB into three major categories, namely, provisioning services, regulating services, and supporting services. And in each category, there are subcategories and further items, each being relevant to the land use change process that characterizes this study; thus an index system of the key ecosystem services of the HRB was built (Tables 1 and 2). Given the fact that the cultural services function is difficult to quantify, this research did not take the cultural function into consideration.Table 1
Key ecosystem services indicators in the HRB.
Categories
Subcategories
Indicators and description
Provisioning  services
Food
Cultivated land percentage %
Fresh water
Water area percentage %
Fuelwood
Forestry area percentage %
Regulating services
Climate regulating
Air temperature: mean annual temperature (0.1°C)
Temperature departure
Precipitation departure
Cumulated temperature (>10°C) Annually cumulated temperature of daily meanAir temperature above 10°C (0.1°C)
Sunshine hours: sunshine hours rectifying the spatial variability of solar  radiation
Sunshine hours departure
Species diversity
Landscape diversity index (Shannon diversity index)
Supporting services
Soil formation
Loam proportion %
Slope: terrain slope derived from DEM (0.01 degrees)
Nutrient cycling
Soil phosphorus content %
Primary production
Net primary production gC m−2yr−1
Grassland area
Grassland percentage %Table 2
Descriptive statistics of key ecosystem services indicators in the Heihe River Basin.
Variable
Unit
Mean
Std. Dev.
Min
Max
Cultivated land
%
5.52
1873.59
0
100
Water area
%
4.49
1616.67
0
100
Grassland
%
16.79
3091.04
0
100
Forestry area
%
1.12
660.01
0
100
Air temperature
°C
0.51
2.54
0
0.96
Temperature departure
°C
0
2.52
−5.59
4.27
Precipitation departure
mm
0
75.04
−137.93
213.55
Cumulated temperature
°C
2479.10
733.29
461.29
3942.56
Sunshine hours
h
8.33
0.35
7.31
9.30
Sunshine hours departure
h
-
2.85
e
-
15
0.31
−1.00
0.76
Landscape diversity index
—
−37.81
116.46
−460.52
460.52
Loam proportion
%
22.50
6.37
4
34
Slope
0.01°
4.72
5.53
0
60.57
Soil phosphorus content
%
0.07
0.01
0.02
0.08
Net primary production
gC m−2yr−1
97.54
112.16
0
897
Note:  total observation: 161584.In this study, the information about land use change was derived from the remote sensing images of the HRB in the years of 1988, 1955, 2000, and 2008 through the human-computer interactive interpretation (Figure2). And the information can be used as the basic data for analyzing the different land use zoning.


## 3.3. Data Processing

Our study uses a land cover classification data set developed by the Chinese Academy of Sciences (CAS) [23, 24] as the background data for spatially identifying the ecosystem services. The data set is derived from Landsat TM/ETM images of bands 3, 4, and 5 with a spatial resolution of 30 m × 30 m. And the classification contains 6 types of land uses/covers, including cultivated land, forestry, grassland, water area, built-up area, and unused land [25].


### 3.3.1. Raster Data Processing

The identification of the key ecosystem services function is based on the raster data model. The indicators of the ecosystem services function are characterized by multi-source and multi-scale. Therefore, the 1 km area percentage data [26, 27] (Figure 3) is an expressive way of data integration and fusion in monitoring, forecasting, and analyzing the regional ecosystem services changes. The procedure to generate the 1 km area percentage data set includes four steps. The first step is to generate a vector map of land cover and land use changes during the study period at a scale of 1 : 100,000 based on the remote sensing Landsat TM/ETM data. The second step is to generate a 1 km FISHNET vector map georeferenced to a boundary map of the study area at a scale of 1 : 100,000. Each cell of the generated 1 km FISHNET vector map owns a unique ID. The third step is to overlay the vector map of land cover and land use changes with the 1 km FISHNET vector map. This is done by aggregating converted areas in each 1 km grid identified by cell IDs of the 1 km FISHNET vector map in the TABLE module of ArcGIS software. Finally, the area percentage vector data are transformed into grid raster data to identify the conversion direction and intensity. The design of working flow insists on no loss of area information. Without special notification, the statistical area of cultivated land according to the GRID data is survey area by satellite remote sensing data, which can be called “gross area.”Figure 3
Flowchart of generating 1 km area percentage data set of land use changes.


### 3.3.2. Indicator Data Processing

Among the indicators that are needed in the identification of the key ecosystem services (Figure4), the raster data of the net primary production are provided by the CAS. The data sources and processing methods of other indicators are listed below.(1) Food, water, and fuelwood supply: in order to precisely characterize the supply of grain, water, and fuelwood in each cell, in this research, we use the percentage of the cultivated land, water area, and forestry to indicate these provisions. The 1 km area percentage data processing method is shown above in Figure3.
(2) Climate variables: all the climatic variables are generated from the site-based observations from the China Meteorological Administration. The spline interpolation algorithm is employed to make the surface data of climatic variables acquired at observation stations [28, 29]. The values for the climatic variables during simulation period are estimated using the space-time stochastic model [30].
(3) Landscape diversity index: we selected the Shannon diversity index to measure the landscape diversity index [31].
The Shannon diversity index(
H
) is expressed as [32] (15)
H
=
∑
i
p
i
ln
⁡
p
i
,where p
i is the share of land use type i in total area of the landscape.
(4) Soil data [33]: soil data were from multiple sources. Firstly, the specific soil survey points were calibrated on local administrative division base map, based on carefully comparing the information on National Secondary Soil Survey data and 1 : 250000 topographic maps. Secondly, an Albert projection value for specific soil survey point was assigned onto the map and generated a soil survey base map. After that, soil attribute data were linked onto soil survey base map by using the Kring interpolation algorithm and generated spatial soil data encompassing soil type and physicochemical characters. Finally, these soil data were saved as ArcGIS Grid format, with soil type, soil depth, soil organic content, and soil texture information.
(5) Elevation data: the elevation data were derived from China 1 : 250000 digital elevation models, with a spatial resolution of 50 m. Average slope information was extracted on each 1 km × 1 km plot by GIS spatial analyzing function; then it was used as basic information for land suitability evaluation.Figure 4
Spatial dispersion and characterization of the ecosystem services indicators, (a) cultivated land percentage, (b) grassland area percentage, (c) air temperature, (d) precipitation, (e) soil phosphorus content, and (f) net primary production.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.3.1. Raster Data Processing

The identification of the key ecosystem services function is based on the raster data model. The indicators of the ecosystem services function are characterized by multi-source and multi-scale. Therefore, the 1 km area percentage data [26, 27] (Figure 3) is an expressive way of data integration and fusion in monitoring, forecasting, and analyzing the regional ecosystem services changes. The procedure to generate the 1 km area percentage data set includes four steps. The first step is to generate a vector map of land cover and land use changes during the study period at a scale of 1 : 100,000 based on the remote sensing Landsat TM/ETM data. The second step is to generate a 1 km FISHNET vector map georeferenced to a boundary map of the study area at a scale of 1 : 100,000. Each cell of the generated 1 km FISHNET vector map owns a unique ID. The third step is to overlay the vector map of land cover and land use changes with the 1 km FISHNET vector map. This is done by aggregating converted areas in each 1 km grid identified by cell IDs of the 1 km FISHNET vector map in the TABLE module of ArcGIS software. Finally, the area percentage vector data are transformed into grid raster data to identify the conversion direction and intensity. The design of working flow insists on no loss of area information. Without special notification, the statistical area of cultivated land according to the GRID data is survey area by satellite remote sensing data, which can be called “gross area.”Figure 3
Flowchart of generating 1 km area percentage data set of land use changes.


## 3.3.2. Indicator Data Processing

Among the indicators that are needed in the identification of the key ecosystem services (Figure4), the raster data of the net primary production are provided by the CAS. The data sources and processing methods of other indicators are listed below.(1) Food, water, and fuelwood supply: in order to precisely characterize the supply of grain, water, and fuelwood in each cell, in this research, we use the percentage of the cultivated land, water area, and forestry to indicate these provisions. The 1 km area percentage data processing method is shown above in Figure3.
(2) Climate variables: all the climatic variables are generated from the site-based observations from the China Meteorological Administration. The spline interpolation algorithm is employed to make the surface data of climatic variables acquired at observation stations [28, 29]. The values for the climatic variables during simulation period are estimated using the space-time stochastic model [30].
(3) Landscape diversity index: we selected the Shannon diversity index to measure the landscape diversity index [31].
The Shannon diversity index(
H
) is expressed as [32] (15)
H
=
∑
i
p
i
ln
⁡
p
i
,where p
i is the share of land use type i in total area of the landscape.
(4) Soil data [33]: soil data were from multiple sources. Firstly, the specific soil survey points were calibrated on local administrative division base map, based on carefully comparing the information on National Secondary Soil Survey data and 1 : 250000 topographic maps. Secondly, an Albert projection value for specific soil survey point was assigned onto the map and generated a soil survey base map. After that, soil attribute data were linked onto soil survey base map by using the Kring interpolation algorithm and generated spatial soil data encompassing soil type and physicochemical characters. Finally, these soil data were saved as ArcGIS Grid format, with soil type, soil depth, soil organic content, and soil texture information.
(5) Elevation data: the elevation data were derived from China 1 : 250000 digital elevation models, with a spatial resolution of 50 m. Average slope information was extracted on each 1 km × 1 km plot by GIS spatial analyzing function; then it was used as basic information for land suitability evaluation.Figure 4
Spatial dispersion and characterization of the ecosystem services indicators, (a) cultivated land percentage, (b) grassland area percentage, (c) air temperature, (d) precipitation, (e) soil phosphorus content, and (f) net primary production.
(a)
(b)
(c)
(d)
(e)
(f)


## 4. Results

### 4.1. Land Use Change

Land use changes, which are mainly caused by human activities, are the major driving force of the ecosystem services changes. The land use changes in the study area from 1988 to 2008 (Table3) were characterized by a general increase of cultivated land, with a decrease in forestry, grassland, and water area. Built-up area recorded the most significant expansion, being the result of rapid urbanization. Among all six types of land use, the unused land accounts for the largest proportion in all of the four study time periods, which is composed mainly by alp and desert and concentrated mainly in the north part of the region, where the major ecosystem service there is climate regulating services.Table 3
Land use changes during 1988–2008 (unit: km2).
Land use type
1988
1995
2000
2008
Change rate during 1988–2008 (%)
Cultivated land
4650
4290
4814
5194
11.70
Forestry
3509
3575
3538
3505
−0.11
Grassland
13382
13642
13186
13269
−0.84
Water area
739
1075
750
701
−5.14
Built-up area
390
365
419
466
19.49
Unused land
17718
17446
17686
17239
−2.70
Note: change rate during (R) 1988–2008 is calculated based on the following formula: R
=
(
B
-
A
)
/
A
∗
100
%; A represents the area of each type of land in 1988 and B represents the area of each type of land in 2008.As for the transformation of land use types (Table4) characterizing the ecosystem structure changes, we can see that the increased cultivated land mainly comes from grassland and used land. During this process, though the ecosystem service of the cultivated land is strengthened, the grassland’s ecosystem services, which are also crucial in providing provision and supporting services, are impaired. Therefore, in land use zoning, it is important to consider not only the amount of the land use conversion but also the direction of the conversion among various land use types.Table 4
Land use transformation matrix during 1988–2008 (unit: km2).
1
2
3
4
5
6
1
1.84
42.61
2.16
23.41
4.06
2
24.99
14.63
0.06
0.34
2.95
3
345.92
6.46
3.71
7.15
10.28
4
45.23
0
5.44
0.98
5.35
5
6
201.16
2.77
41.41
12.99
31.83
Note: 1: cultivated land, 2: forestry, 3: grassland, 4: water area, 5: built-up area, and 6: unused land.The spatial land use change (Figure5) mainly occurred in the midnorth, midcentral, and midsouth area, which could explain the boundary shifting of the ecosystem services zone in these areas. Through studying the historical land use zoning policy and its impact on ecosystem services, we aim at optimizing the spatial planning of land use to conserve ecosystem services under the impact of climate change.Figure 5
Spatial land use changes during 1988–2008.


### 4.2. Ecosystem Service Zoning Change

In identifying and zoning of ecosystem services, firstly, we applied the analytic hierarchy process to build an indicator system (Table1) with the conceptual framework provided by MEA for ecosystem services assessment with 15 indicators that proved to be the major driving forces of ecosystem services in the middle reaches of the HRB. The data of these indices were then spatially represented with the 1-kilometer (km) grid pixel data model. Then using a correlation matrix, factors with eigenvalues >1 were retained, and factors were subjected to varimax rotation to maximize the correlation between factors and measured attributes. The magnitude of the eigenvalues was used as a criterion for interpreting the relationship between ecosystem services and factors. Afterwards, to summarize and aggregate the selected ecosystem services attributes, a principal component evaluation was performed on the aforementioned variables. We assumed that principal components receiving a high contribution best represented the core ecosystem services. Therefore, we retained only the principal components with a cumulative contribution greater than 60%. An unsupervised fuzzy clustering algorithm was then performed using values for the retained principal components to develop the functional zones for ecosystem services, with subsequent valuing of each 1 km grid cell after the iterations reaching their threshold. Finally, we grouped data in clusters, where within each cluster the data exhibit similarity. Besides, in this paper, the center of gravity of each service zone in each time is generated in the ArcGIS to analyze the shifts of center of gravity (Figure 7).In 1988 and 1995, there are five major clusters, namely, climate regulating function, food provision function, nutrient cycling supporting function, soil supporting function, and supporting function. The name of each cluster is determined according to the above mentioned principal component analysis within each cluster. Similarly, four clusters have been identified in 2000 and 2008, which are climate regulating function, food provision function, soil supporting function, and supporting function.Each functional zone involves one kind of ecosystem services, which is displayed with different zoning boundaries. The functional zones for ecosystem service obtained with SIZES comprehensively reflect the spatial clusters of ecosystem services in the middle reaches of the HRB. There were some uncertainties, but the result obtained with the SIZES still reflected the conditions of local ecosystem services in the HRB; the result indicated the boundaries of functional zones for ecosystem services were conspicuously different from the administrative boundaries (Figure6). The functional zones generally expand across several counties; there is generally one major kind of ecosystem services as well as several kinds of secondary ecosystem services in each county.Figure 6
Ecosystem service zoning of middle reaches of the HRB in the years of (a) 1988, (b) 1995, (c) 2000, and (d) 2008.
(a)
(b)
(c)
(d)Figure 7
Spatial shifts of the center of gravity of ecosystem service zones. Note: 1: supporting function, 2: nutrient cycling supporting function, 3: soil supporting function, 4: fuelwood provision function, 5: climate regulating function, and 6: food provision function.The ecosystem service zoning results (Figure6) showed that, during 1988–2008, the middle reaches maintained the climate regulating function, which is the major function of the area mainly distributed in the north region. It plays a key role in mitigating the climate change. The boundaries of the climate regulation service zone changed slightly and the center of gravity (Figure 7) changed the most among all the other services in this region, which is a result of the changing in land use structure and climate and its own fragile ecological condition in the northern part.The food provision function zone is mainly concentrated in the cultivated land area. From 1988 to 1995, the food provision function is decreasing in response to the cultivated land degradation during the same time period. And by 2008, this provision increased with the expansion of the cultivated land area. The center of gravity also shifted between each two of the study time period, though being the slightest shift among others. In 2000, a large area of fuelwood supply function zone occurred in the southern part, which is considered as the result of forestry land increase.The supporting function showed an overall declining trend characterized by boundary shrinking, which is the result of the net primary production loss and the decreased soil function due to the forestry and grassland reduction. In the southern part, the shift from supporting function to fuelwood provision function after 2000 is the result of the forestry increase in this region; fuelwood provision took place as the major ecosystem service.


## 4.1. Land Use Change

Land use changes, which are mainly caused by human activities, are the major driving force of the ecosystem services changes. The land use changes in the study area from 1988 to 2008 (Table3) were characterized by a general increase of cultivated land, with a decrease in forestry, grassland, and water area. Built-up area recorded the most significant expansion, being the result of rapid urbanization. Among all six types of land use, the unused land accounts for the largest proportion in all of the four study time periods, which is composed mainly by alp and desert and concentrated mainly in the north part of the region, where the major ecosystem service there is climate regulating services.Table 3
Land use changes during 1988–2008 (unit: km2).
Land use type
1988
1995
2000
2008
Change rate during 1988–2008 (%)
Cultivated land
4650
4290
4814
5194
11.70
Forestry
3509
3575
3538
3505
−0.11
Grassland
13382
13642
13186
13269
−0.84
Water area
739
1075
750
701
−5.14
Built-up area
390
365
419
466
19.49
Unused land
17718
17446
17686
17239
−2.70
Note: change rate during (R) 1988–2008 is calculated based on the following formula: R
=
(
B
-
A
)
/
A
∗
100
%; A represents the area of each type of land in 1988 and B represents the area of each type of land in 2008.As for the transformation of land use types (Table4) characterizing the ecosystem structure changes, we can see that the increased cultivated land mainly comes from grassland and used land. During this process, though the ecosystem service of the cultivated land is strengthened, the grassland’s ecosystem services, which are also crucial in providing provision and supporting services, are impaired. Therefore, in land use zoning, it is important to consider not only the amount of the land use conversion but also the direction of the conversion among various land use types.Table 4
Land use transformation matrix during 1988–2008 (unit: km2).
1
2
3
4
5
6
1
1.84
42.61
2.16
23.41
4.06
2
24.99
14.63
0.06
0.34
2.95
3
345.92
6.46
3.71
7.15
10.28
4
45.23
0
5.44
0.98
5.35
5
6
201.16
2.77
41.41
12.99
31.83
Note: 1: cultivated land, 2: forestry, 3: grassland, 4: water area, 5: built-up area, and 6: unused land.The spatial land use change (Figure5) mainly occurred in the midnorth, midcentral, and midsouth area, which could explain the boundary shifting of the ecosystem services zone in these areas. Through studying the historical land use zoning policy and its impact on ecosystem services, we aim at optimizing the spatial planning of land use to conserve ecosystem services under the impact of climate change.Figure 5
Spatial land use changes during 1988–2008.


## 4.2. Ecosystem Service Zoning Change

In identifying and zoning of ecosystem services, firstly, we applied the analytic hierarchy process to build an indicator system (Table1) with the conceptual framework provided by MEA for ecosystem services assessment with 15 indicators that proved to be the major driving forces of ecosystem services in the middle reaches of the HRB. The data of these indices were then spatially represented with the 1-kilometer (km) grid pixel data model. Then using a correlation matrix, factors with eigenvalues >1 were retained, and factors were subjected to varimax rotation to maximize the correlation between factors and measured attributes. The magnitude of the eigenvalues was used as a criterion for interpreting the relationship between ecosystem services and factors. Afterwards, to summarize and aggregate the selected ecosystem services attributes, a principal component evaluation was performed on the aforementioned variables. We assumed that principal components receiving a high contribution best represented the core ecosystem services. Therefore, we retained only the principal components with a cumulative contribution greater than 60%. An unsupervised fuzzy clustering algorithm was then performed using values for the retained principal components to develop the functional zones for ecosystem services, with subsequent valuing of each 1 km grid cell after the iterations reaching their threshold. Finally, we grouped data in clusters, where within each cluster the data exhibit similarity. Besides, in this paper, the center of gravity of each service zone in each time is generated in the ArcGIS to analyze the shifts of center of gravity (Figure 7).In 1988 and 1995, there are five major clusters, namely, climate regulating function, food provision function, nutrient cycling supporting function, soil supporting function, and supporting function. The name of each cluster is determined according to the above mentioned principal component analysis within each cluster. Similarly, four clusters have been identified in 2000 and 2008, which are climate regulating function, food provision function, soil supporting function, and supporting function.Each functional zone involves one kind of ecosystem services, which is displayed with different zoning boundaries. The functional zones for ecosystem service obtained with SIZES comprehensively reflect the spatial clusters of ecosystem services in the middle reaches of the HRB. There were some uncertainties, but the result obtained with the SIZES still reflected the conditions of local ecosystem services in the HRB; the result indicated the boundaries of functional zones for ecosystem services were conspicuously different from the administrative boundaries (Figure6). The functional zones generally expand across several counties; there is generally one major kind of ecosystem services as well as several kinds of secondary ecosystem services in each county.Figure 6
Ecosystem service zoning of middle reaches of the HRB in the years of (a) 1988, (b) 1995, (c) 2000, and (d) 2008.
(a)
(b)
(c)
(d)Figure 7
Spatial shifts of the center of gravity of ecosystem service zones. Note: 1: supporting function, 2: nutrient cycling supporting function, 3: soil supporting function, 4: fuelwood provision function, 5: climate regulating function, and 6: food provision function.The ecosystem service zoning results (Figure6) showed that, during 1988–2008, the middle reaches maintained the climate regulating function, which is the major function of the area mainly distributed in the north region. It plays a key role in mitigating the climate change. The boundaries of the climate regulation service zone changed slightly and the center of gravity (Figure 7) changed the most among all the other services in this region, which is a result of the changing in land use structure and climate and its own fragile ecological condition in the northern part.The food provision function zone is mainly concentrated in the cultivated land area. From 1988 to 1995, the food provision function is decreasing in response to the cultivated land degradation during the same time period. And by 2008, this provision increased with the expansion of the cultivated land area. The center of gravity also shifted between each two of the study time period, though being the slightest shift among others. In 2000, a large area of fuelwood supply function zone occurred in the southern part, which is considered as the result of forestry land increase.The supporting function showed an overall declining trend characterized by boundary shrinking, which is the result of the net primary production loss and the decreased soil function due to the forestry and grassland reduction. In the southern part, the shift from supporting function to fuelwood provision function after 2000 is the result of the forestry increase in this region; fuelwood provision took place as the major ecosystem service.


## 5. Discussion and Conclusions

This paper used a model to identify and zone the ecosystem services to explore the land use zoning on ecosystem services in climate mitigation. The research contributes to the growing literature in finely characterizing the ecosystem services zones via land use zoning to mitigate the impact of climate change, as there is no such systematic ecosystem zoning method before.The application of SIZES proves capable of automatically identifying the boundary of ecosystem service function zone. The method of ecosystem services spatial clustering at watershed scale can be applied in other regions. Though in our study in the middle reaches of HRB we focused on the climate regulation ecosystem services, our method can be further applied in other regions to explore their specific typical ecosystem services. And the results will provide scientific basis for the formulation of regional ecosystem services development planning and ecosystem adaptive management decision. However, the spatial mapping capability of the software needs to be further strengthened.As one of the largest inland river basins, the provision of the ecosystem services in the HRB is of great significance to the production and life of the local residents. And the regulating services are of crucial importance for this ecological fragile region to counter the impact of climate change. In this research, the impact of land use change on ecosystem services was explored in the middle reaches of the HRB. Through the identification of the ecosystem service, we recognize that, in the middle reaches of the HRB, it is important to protect the provision and regulating services, especially given the fragile ecology in this region. We identify that it is a direct and effective way of land use zoning to protect the key ecosystem service and raise the provision and regulating function.Identifying major regional ecological environment problems and thereafter establishing the prioritized ecosystem conservation zone are of practical significance. In the middle reaches of the HRB, the major provision function of ecosystem is food provision function, which is mainly concentrated in the cultivated land area and is strengthened with the expansion of cultivated land during the studied time period. In 2000, a large area of fuelwood provision function zone occurred in the southern part, which is considered as the result of forestry land increase. The supporting services mainly occurred in the grassland and water area and show a general trend of degradation, as a result of the net primary production loss and the decreased soil function due to the forestry and grassland reduction. The major regulation service in this study area is climate regulation services, mainly concentrated in the northern desert, Gobi, and grassland area. The conservation of major ecosystem services should be prioritized. And land use zoning serves as an effective way to protect the key ecosystem services.



---

*Source: 496942-2015-08-09.xml*