Spatio-temporal variation of ecosystem services value in the Northern Tianshan Mountain Economic zone from 1980 to 2030

Rapid agricultural land expansion and urbanization have accelerated land use and land cover changes (LUCC) in the Northern Tianshan Mountain Economic Zone and have significantly impacted on the ecosystem services (ESs). However, the spatiotemporal variations of ecosystem service value (ESV) to LUCC are not well understood. Based on the land use and land cover (LULC) data from 1980 to 2019, we used a CA-Markov model to predict LUCC in 2020 and 2030, assess the spatial-temporal changes of ESV and LULC during 1980–2030, and explore the elastic response of ESV to LUCC. We found that cropland and built-up land expanded rapidly by 34.38% and 196.66%, respectively between 1980 and 2030, while grassland and unutilized land decreased significantly by 11.45% and 10.26%, respectively. The ESV of water body, cropland, grassland and forestland accounts for more than 90% of the total ESV. Our research shows that the ESV of cropland increased 32 million yuan from 1980 to 2030, mainly due to the expansion of cropland area. However, the loss caused by the reduction of grassland area was 45 million yuan. Water conservation, waste treatment, soil formation and retention, and biodiversity conservation are the primary ecosystem service function, accounting for 71.82% of the total ESV. Despite notable increases in the ESV from 1980 to 2010, grassland degradation still remains a main ecological and environmental issue from 2010 to 2030. The results suggest that effective land use policies should be developed to control the expansion of croplands and protect water body, grassland and forestland to maintain more sustainable ESs.


INTRODUCTION
Ecosystem service (ES) refers to the benefits directly or indirectly obtained by human through ecosystem structure, process and function, which includes four aspects: provisioning services, regulating services, supporting services and cultural services (Costanza et al., 1997). By evaluating the ecosystem service value (ESV), the benefits obtained from the ecosystem can be quantified to help decision makers make optimal decision on the rational allocation of resources, so as to realize the sustainable management of the ecosystem (Deal, Smith & Gates, 2017). With increasing importance of ES, the evaluation methods of ESV become more and more diversified. Talberth (2015) divided the evaluation methods into four categories, including stated preference method, revealed preference method, cost-based method and benefit transfer method (BTM). A series of models have also been developed to accurately evaluate ESV, including Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) (Redhead et al., 2016), Social Values for Ecosystem Services (SolVES) (Sherrouse, Semmens & Clement, 2014;Zhou et al., 2020), and so on (Sun et al., 2019). In regional or global ESV assessment, BTM is simple and easy to operate and has been widely used. This approach potentially assumes that the value of each ecosystem service function (ESF) is derived from a single or a multiple case study utilizing the specific value of a particular land cover (Costanza et al., 1997). Recently, the evaluation has been updated based on more than 300 case studies worldwide (Costanza et al., 2014;De Groot et al., 2012). Costanza et al. (2014) claimed that the models can be applied at various scales. However, the spatial heterogeneity of ecosystem structure eventually leads to the spatial heterogeneity of ESF (Wilson & Howarth, 2002;Kindu et al., 2016). Since the limitations and uncertainties of land use and land cover changes (LUCC) are ignored, and simply using the global value coefficient to estimate the small-scale ESV limits the accuracy of the assessment . Based on the results of Costanza et al. (2014) and Xie et al. (2010), it is an important issue to adjust the equivalent factor to a small scale and reflect the changing relationship between LUCC and ESV based on the social and economic status quo.
Land use and land cover changes can affect the supply and value of ESs by changing the structure, processes and functions of ecosystems (Daily, 1997). Overuse of land resources may lead to degradation of regional ecosystem and the loss of ESV. Recent studies have shown that with rapidly increasing population, land pressure is gradually increasing (Han, Song & Deng, 2016). Excessive land reclamation and urbanization leads to land degradation, loss of biodiversity, water shortage and carbon loss, and further leads to a significant decline in ESV (Collin & Melloul, 2001).
The Northern Tianshan Mountain Economic Zone (NTMEZ) is an important hub of the "Silk Road" from Asia to Europe, and the largest and most developed economic area in Xinjiang. However, the special geographical environment of the NTMEZ has resulted in limited water resources and a fragile ecosystem (Zhou, Zhao & Yang, 2017). Under the influence of long-term agricultural activities and the development and utilization of water and land resources, the ecological and environmental problems caused by land use in this oasis are becoming more and more prominent . Under China's "One Belt One Road" strategy, the NTMEZ has entered an important stage of development. According to the development plan of urban agglomeration for the NTMEZ in Xinjiang, the urban agglomeration will build an urban system with Urumqi as the core area, which will accelerate population and industrial agglomeration. The rapid urbanization will directly lead to decrease of grassland and forestland, and destruction of ecosystem structure, function, and services in the NTMEZ. Therefore, it is urgent to evaluate the impacts of human disturbances on ES in NTMEZ. In addition to the oasis LUCC and urbanization, there has been no quantitative spatial and temporal assessment of ESV in response to LUCC in NTEMZ.
Based on the research gap identified above, this study quantitatively evaluates the LUCC and ESV of NTMEZ from 1980 to 2030. There are three specific objectives: (1) to reveal the spatial and temporal variations of the land use and land cove (LULC) and ESV from 1980 to 2030; (2) to evaluate changes in ESV in response to LUCC; and (3) to explore the elasticity of the response of ESV to LUCC by 50% adjustment of value coefficients. These results will provide information to government for the sustainable development of land use and ecological environment.

Study area
The NTMEZ is located in the north of the Tianshan Mountain and the south of the Junggar basin in Xinjiang, China (Fig. 1), covers an area of 9.54 × 10 4 km 2 , accounting for 5.70% of Xinjiang (Xinjiang Bureau of Statistics, 2017), and is one of the 18 national key development areas. The population of the study area was 8.92 million in 2016 (Xinjiang Bureau of Statistics, 2017), accounting for 37.80% of the population of Xinjiang. Notably, the Gross Domestic Product (GDP) of the study area was 64.62 × 10 11 yuan in 2016, accounting for 69.30% of the GDP of Xinjiang in that year (Xinjiang Bureau of Statistics, 2017). Among them, the yields of primary industry, secondary industry and tertiary industry were 6.22 × 10 11 yuan (9.62%), 25.34 × 10 11 yuan (39.21%) and 33.06 × 10 11 yuan (51.17%), respectively. With the implementation of China's "Western Development Strategy", this region has become the most economically developed region in Xinjiang and an important hub of the "Silk Road Economic Belt". The annual average temperature, average annual precipitation, and potential evaporation in the study area are 6.10-8.90 C , 220 mm and 1,210 mm, respectively (Zhao et al., 2010). The evaporation of the study area is much higher than the precipitation, the climate belongs to the temperate arid continental climate, and the ecological environment is extremely fragile. With increase in human activity, the resource and environmental problems in this region have become increasingly prominent.

Data preparation
The Landsat TM/ETM remote sensing image data of years: 1980, 1990, 2000, 2010, 2015 and 2019, were obtained from the U.S. Geological Survey (https://www.usgs.gov/) with a spatial resolution of 30 m. The radiometric correction, geometrical correction, and atmospheric correction of the images were done using the ERDAS 9.2 software. In this article, we referred to the ESV classification system presented by Xie et al. (2010). Meanwhile, considering the practical situation of the study area, supervised classification and visual interpretations were combined in order to classify the LULC types into six categories: cropland, forestland, grassland, water body, built-up land and unutilized land. The classification basis of LULC is shown in Table S1. There were more than 40 training sites for each LULC type for classification process, which were obtained using a GPS device; Google Earth was applied for validation of the LULC type (Elkhrachy, 2015). The confusion matrix method was used to evaluate the interpretation accuracy, and the Kappa coefficient was 87%. We also collected the digital elevation modal, meteorological data, the distribution data of the population and GDP, basic geographic information data, and grain yield data. Data descriptions and resources are shown in Table 1.

CA-Markov model
CA model is a lattice dynamic model with discrete time and space states, which focuses on the interaction between cells with different spatial and temporal characteristics . It has strong spatial computational simulation ability and is especially suitable for dynamic simulation and spatial display. The Markov model focuses on the quantitative prediction of LUCC, but it cannot carry out spatial expression and show the spatial distribution of LUCC (Wickramasuriya et al., 2009). However, CA can express the spatial-temporal dynamic evolution process of complex space system, which can make up the deficiency of the Markov model. By combining the CA model with Markov model, the CA-Markov model has the advantage of expressing the temporal and spatial patterns of LUCC, so as to better simulate the spatial and temporal patterns of future regional LUCC, and further provide support for the future land use and ecological environment sustainable development (Mitsova, Shuster & Wang, 2011). The Markov chain was constructed based on the probability of change matrices for LUCC from t to t+1: is the state probability of any time, and S (t) is the initial state probability. A is the transition probability matrix, and the formula is as follows: where A ij is the sum of areas from the ith to the jth land cover category from the initial to the forecast period, and n is the number of LULC categories. We accomplished this process by using the Markov model with the IDRISI software that integrates geographic information system and image processing functions.
The CA model was used to predict the spatial and temporal dynamic change pattern using a transition map of the LUCC (Yang, Zheng & Lv, 2012), and the model can be defined as follows (Wang, Zheng & Zang, 2012): where S is a set of cellular states, N is the cellular field, t and t+1 represent different time periods, and f is the local transition rule of the cell.
To ensure the reliability of the simulation results, we used the kappa index to test the consistency level of the simulated and observed land cover maps (Mitsova, Shuster & Wang, 2011): where Kappa is the index of simulation accuracy, P c is the expected simulation accuracy in a random state, and P 0 is the actual simulation accuracy.
In this study, the transition probability matrix from 2005 to 2010 was used to predict the LUCC from 2015 to 2019. The transition probability matrix is showed (Table S2). Then, the forecast data in 2015 and 2019 were compared with the respective actual data to verify the reliability of the CA-Markov model (Fig. 2). The Kappa coefficient was 0.79 and 0.72 in 2015 and 2019, respectively, which were above the required standard of 70%. The simulation results were accurate and reliable, can be used for prediction and simulation of future LUCC. Finally, the LUCC spatial data in 2020 and 2030 were predicted.

Assessment of ESVs
In this article, the estimation coefficients of ESV were referred to Costanza et al. (1997) and Xie et al. (2010) were referred to. Given the spatial and temporal effects of ES, the ESV for LULC type was revised according to the actual conditions of the study area by consulting 30 ecologists. According to the specific conditions of the research area, the economic value of the annual natural grain yield per unit area of cropland was modified as follows: the average annual grain yield of the research area from 1990 to 2010 was used to replace the grain yield per unit area of each year; the average purchase price of grain (1.98 yuan/kg) in China in 2010 was used as the grain price for the research period (China Food Industry Association, 2010); the equivalent factor of the ESV of the research area was calculated as 1,304.86 yuan/hm 2 . Meanwhile, considering the negative impacts of built-up land on the ecosystem in terms of water supply and waste treatment, the replacement cost method was used to estimate the ESV of built-up land (Mayila et al., 2018). On this basis, the ESV coefficients of LULC in the study area were calculated ( Table 2; Table S3). Notably, the equivalent value of ESs per standard unit is defined as 1/7 of the economic value of the average annual grain yield of 1 hm 2 of cropland. The formula is given as follows: where E a is the economic value of 1 hm 2 farm's annual grain crops (yuan/hm 2 ); P i is the national average price of i food crops (yuan/t); q i represents the yield of grain crop i (t/hm 2 ); M i is the area of i food crops (hm 2 ); M is the total area of food crops (hm 2 ); i is the crop type, including rice, wheat and corn. The calculation formula of the ESV is as follows: where A i is the distribution area of land use type i in the research area (hm 2 ); VC i is the ESV coefficient (yuan·hm −2 ·a −1 ); ESV f represents the functional value of individual f ecosystem services; VC if represents the functional value of individual f ecosystem services of land use type i (yuan·hm −2 ·a −1 ).

Elasticity for the response of ESV to LUCC
We used ESV coefficients of LULC to replace the six LULC categories that did not exactly match the report by Costanza et al. (1997) (Table S3), which resulted in uncertainties in the assessment of the ESV. Therefore, we used sensitivity analysis to evaluate the changes in ESV in response to 50% adjustments of the ESV coefficients for each LULC type  (Mamat et al., 2019). The standard economic concept of elasticity was used to calculate the coefficient of sensitivity (CS) using the following formula: where ESV is the estimated total value of ecosystem services, VC is the value coefficient, and "i", "j" and "k" represent the initial, adjusted values, and LULC categories, respectively. If CS > 1, then the estimated ESV is elastic with respect to that coefficient; if CS ≤ 1, the estimated ESV is inelastic. Thus, when CS < 1, even if the accuracy of VC values used as proxy biomes is low, the results of estimation of ESV are credible.

Changes of LULC
According to the simulation results of the CA-Markov model, the spatial distribution of LULC is showed in Fig. 3 and the dynamic change of land use area is showed in Table 3.
In 1980, LULC was dominated by grassland and unutilized land, which accounted for 42.73% and 38.83% of the total area, respectively, followed by cropland ( Due to water shortage in the arid areas, the water body only accounted for 3% of the total area. Water body increased from 3,093.59 km 2 in 1980 to 3,392.38 km 2 in 2010. However, with the increase of cropland area, water consumption will further increase, and water body will further decrease 1.03% during 2010-2030.

Changes of ESVs
According to the ESV of different types of LULC per unit area (Table 2), we calculated the ESV of different LULC in NTMEZ from 1980 to 2030 (Table 4). The total ESV in NTMEZ was approximately 7.83 × 10 8 yuan in 1980. Grassland had the highest contribution of 49.43%, followed by cropland, water body and forestland with a contribution of 23.75%, 13.54% and 11.75%, respectively. Due to LUCC, the ESV increased by 2.35% from 1980 to 2010. Cropland and water body increased and played an important role in improving ESV, which over compensated the ESV loss in grassland and built-up land (Fig. 4). Notably, with the acceleration of urbanization, grassland and unutilized land gradually decreased, and the ESV will decrease by 0.3 × 10 8 yuan from 2010 to 2030.
The spatial distribution of the ESV in the study area was unbalanced. The high ESV area was mainly distributed in Aibi Lake Basin, Selimu Lake, Baiyang River, Toutun River and Manas River Basin (Fig. 5). In addition to the higher ESV distribution along river water, the northern part of the Tianshan Mountain also had a higher distribution of ESV. The main distribution of the lowest ESV areas was located in the edge of the Gurbantunggut Desert. Notably, compared with the period from 1980 to 2010, the decrease of ESV between 2010 and 2030 shifted from Shihezi, Hutubi, Changji and Jimusaer to Karamay, Fukang, Urumqi and Manas (Table 5). With the urban population growth and industrial agglomeration, the ecological environment will gradually deteriorate during 2010-2030.

Changes in values of ESF
The contribution and change of individual ESF are summarized in Table 6. Water conservation, waste treatment, soil formation and protection, and biodiversity conservation contributed the most to the ESV, with a contribution of 20.14%, 20.67%, 17.22% and 13.80% in 2030. Most of individual ESFs increased except food production during 1980-2010 with cropland expansion. Water conservation, food production, raw material and recreational culture increased by 5.73%, 4.92%, 5.70% and 7.66%, respectively (Fig. 6). However, most of individual ESFs were projected to decrease from 2010 to 2030. Due to the increase of cropland area, only the ESV of food production will increase during 2010-2030.

Ecosystem sensitivity analysis
In this article, the value index of LULC type was adjusted by 50% to measure the change of total ESV (Table 7). The results showed that the sensitivity index of ESV coefficient of all LULC types was less than 1, and the ESV of LULC was inelastic to VC, indicating that the ESV coefficient had little impact on the change of ESV. The sensitivity index of unutilized land was the smallest and the sensitivity index of grassland was the largest, indicating that grassland contributed the most to ESV. In addition, from the change of the sensitivity index in different periods, the CS value increased significantly due to the large increase of the cropland area, while other sensitivity index did not change significantly. The sensitivity index of unutilized land and built-up land was close to zero, and the VC change of the two land types had little impact on the change of ESV in NTMEZ. Overall, this study constructed the ESV coefficient suitable for the local actual situation.

Impact of LUCC change on ESV
Due to the acceleration of urbanization, population growth and agricultural intensification globally, the LULC is changing rapidly, which may lead to changes in ES (Nahuelhual et al., 2014). These effects are particularly significant in arid regions. Zhang et al. (2020) shows that the population showed a sharp increase in NTMEZ from 1980 to 2015, especially during the "Western Development" period, which led to a sharp increase in the population in 2001. Population growth led to the overexploitation of water body, cropland and built-up land to meet demands of water, food, energy and land (Granit et al., 2012). This study found rapid expansion of cropland (34.38%) and built-up land (196.66%), and shrinking of grassland (11.45%) and unutilized land (10.26%) during 1980-2030. Correspondingly, the ESV of cropland increased 0.32 × 10 8 yuan from 1980 to 2030, which was mainly caused by the expansion of cropland area. However, the area of grassland decreased sharply during 1980-2030, causing a loss of 0.45 × 10 8 yuan. Increase in cropland may yield short-term economic benefits, but large increases in cropland may result in loss of ES (Li et al., 2019). Although short-term economic benefits can be gained from expanding agricultural production to meet the needs of population growth, it also reduces the ESV of grassland and unutilized land (Coupe, Barlow & Capel, 2012). In this study, we found the increase of cropland has led to increase in ESV of raw material by 9.12% from 0.307 × 10 8 yuan to 0.335 × 10 8 yuan and food production by 3.00 % from 0.167 × 10 8 yuan to 0.172 × 10 8 yuan in NTMEZ during 1980-2030. However, the ESV of soil formation and retention, biodiversity conservation and gas regulation were projected to decrease. These findings are consistent with many previous studies around the world (Sutton et al., 2016;Qiao et al., 2019). Agriculture and urban sprawl have a negative impact on the provision of ESFs, especially in Shihezi, Hutubi and Fukang. Our results showed that the ESV in Shihezi, Hutubi and Fukang, decreased by more than 20% from 1980 to 2030, which was mainly caused by the increase of cropland area and shrinkage of forestland and grassland area. Therefore, a series of protection measures should be carried out urgently. We suggest that a "Red Line for Ecological Protection" should be drawn and ecological land should be integrated into the land use classification system in order to restrict the encroachment of cropland on ecological land.

Limitations and future work
Remote sensing images are the most important data source in our ESV evaluation. However, during the data processing, the spatial and temporal resolution, spectral resolution and classification method of LULC data will affect the interpretation accuracy. ESV is inevitably influenced by errors in the process of quantifying LULC data (Sexton et al., 2015). In this study, a land with more than 5% covered by plants is classified as grassland, which includes shrub pastures and mixed pastures with shrub canopy coverage less than 10% (Table S1). This significantly overstates the grassland area, leading to an overestimation of ESV. To address these limitations, higher spatial resolution remote sensing data and more accurate LULC classifications will be used to further accurately assess ESV. Some ecosystem service costs are ignored when evaluating the ESV, such as soil salinization (Schäfer et al., 2012), loss of genetic resources (Leroy et al., 2018), and eutrophication of water (Wen & Théau, 2020), so our evaluation results may overestimate   Cao et al. (2018) found that the ESV in China decreased by 52.66% in 2014 by considering such factors as water resource cost, ecological protection investment and land management. In the future research, we need to consider the multi-dimensional ES category to improve the evaluation accuracy by combining with the current situation of ecological environment in the research area.
In this study, the CA-Markov model was used to predict the spatial pattern of the ESV, and the reliability of the model was verified. The advantages of using the CA-Markov model in predicting large-scale ESV were demonstrated. However, the study is based on the assumption that the land use policy is basically unchanged , and the natural and human driving factors of LUCC are selected based on the existing research results and experience (López-Marrero et al., 2011). The lack of quantitative analysis of the assumed conditions and the influence of various factors may cause certain e rrors in the prediction results (Kindu et al., 2016;Cao et al., 2018).
The simple benefit transfer approach used in this article to estimate the ESV has some limitations (Schulp et al., 2014). This equivalent factor method is only a static evaluation method, which lacks consideration of spatial and temporal differences of ecosystem types and quality, and the estimation results are insufficient to reflect the dynamic changes of ESFs in spatial and temporal scale, which limits the accuracy of the evaluation (Schmidt, Manceur & Seppelt, 2016;Cao et al., 2018). For example, in the study, we assumed that the cropland ESV had homogeneity in the whole region, and extended the unit value of one region to all regions. In fact, cropland with different crop growing structures provides different ecosystem services functions (Zheng, Zhang & Cao, 2018). In a recent study, a series of ES evaluation models were used to reveal the relationship between land use structure and ecosystem, including InVEST (Ouyang et al., 2016), SolVES (Sherrouse, Semmens & Clement, 2014;Zhou et al., 2020), and Coupling Coordination Model (Sun et al., 2019). For example, Ouyang et al. (2016) quantitatively evaluated the types of ESV in China based on the InVEST model and LUCC, revealing the imbalance between supply and demand. Zhang et al. (2017) used the coupling coordination model to evaluate the coordination degree of LUCC and ESV in Beijing from 2006 to 2015. These models can accurately simulate the dynamic changes of ES such as carbon storage, water and soil conservation, and cultural services in terrestrial ecosystems according to the spatial and temporal changes of LULC (Leh et al., 2013;Zhou et al., 2020;Wen & Théau, 2020). Due to the large amount of data and diverse types of data needed for ESF assessment models, the difficulty of data acquisition limits the application of these models in arid areas (Zhang, Yushanjiang & Jing, 2018). Therefore, there are uncertainties in our assessment results.
In the future, we will find more advanced approach to integrate different data sources to improve the accuracy in estimating ESV.

CONCLUSIONS
The study results showed that the cropland ESV increased about 14 million yuan during 1980-2010. This increase was mainly attribute to the increase of cropland area (15.31%). However, between 1980 and 2030, the grassland area decreased sharply by 3.61%, resulting in a loss of 14 million yuan. It must be pointed out that changes in cropland may yield short-term economic benefits but unsustainable agricultural production and rapid urbanization may lead to loss of natural ecosystem services. This trend will continue to intensify in the future, and ESV will continue to gradually decrease 30 million yuan with the urban population growth and industrial agglomeration from 2010 to 2030. Ecological protection red lines should be made to control farmland expansion and protect water bodies, grassland and forestland for more sustainable ecosystem services.