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

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⁴ km², 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¹¹ 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¹¹ yuan (9.62%), 25.34 × 10¹¹ yuan (39.21%) and 33.06 × 10¹¹ 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 (Chen et al., 2020). 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: $$S_{\left(t+1\right)}=S_{\left(t\right)}\times A$$where S_(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: $$\mathrm{A}=\left[\begin{array}{c}\\ A_{11}{A}_{12}\dots A_{1n}\\ A_{21}{A}_{22}\dots A_{2n}\\ \dots \dots \dots \dots \\ A_{n1}{A}_{n2}\dots A_{nn}\end{array}\right]$$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): $$S_{\left(t,t+1\right)}=f\left(S_{\left(t\right)},N\right)$$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): $$\mathrm{K}\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{a}={\displaystyle \frac{P_0-P_c}{1-P_c}}$$where Kappa is the index of simulation accuracy, P_c is the expected simulation accuracy in a random state, and P₀ 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². 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² of cropland. The formula is given as follows: $$E_a={\displaystyle \frac{1}{7}\sum _{i=1}^n{\displaystyle \frac{m_i\times p_i\times q_i}{M}(\mathrm{i}=1,2,3\dots \mathrm{n})}}$$where E_a is the economic value of 1 hm² farm’s annual grain crops (yuan/hm²); P_i is the national average price of i food crops (yuan/t); q_i represents the yield of grain crop i (t/hm²); M_i is the area of i food crops (hm²); M is the total area of food crops (hm²); i is the crop type, including rice, wheat and corn.

The calculation formula of the ESV is as follows: $$\mathrm{E}\mathrm{S}\mathrm{V}=\sum A_i\times \mathrm{V}{\mathrm{C}}_{\mathrm{i}}$$ $$\mathrm{E}\mathrm{S}{\mathrm{V}}_f=\sum A_i\times \mathrm{V}{\mathrm{C}}_{\mathrm{i}\mathrm{f}}$$where A_i is the distribution area of land use type i in the research area (hm²); VC_i is the ESV coefficient (yuan·hm⁻²·a⁻¹); 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⁻²·a⁻¹).

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: $$\mathrm{C}\mathrm{S}={\displaystyle \frac{\left(\mathrm{E}\mathrm{S}{\mathrm{V}}_j-\mathrm{E}\mathrm{S}{\mathrm{V}}_i\right)/\mathrm{E}\mathrm{S}{\mathrm{V}}_i}{\left(\mathrm{E}\mathrm{S}{\mathrm{V}}_{jk}-\mathrm{E}\mathrm{S}{\mathrm{V}}_{ik}\right)/\mathrm{E}\mathrm{S}{\mathrm{V}}_{ik}}}$$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 (10.65%) and forestland (3.90%). During 1980–2010, cropland and built-up land increased substantially. The cropland expanded 1,562.31 km² with an increase rate of 15.31% by 2010. With the increase of urbanization, built-up land also increase from 643.15 km² in 1980 to 1,095.13 km² in 2010, with an increase rate of 70.28%. Cropland and built-up land were expected to continue to increase during 2020–2030. However, the spatial distribution of unutilized land mainly concentrated in the Gurbantunggut Desert, which decreased from 37,206.00 km² to 35,713.36 km² during 1980–2010, and it will be further decrease to 33,390.00 km² in 2030. The coverage of grassland mainly concentrated in the north slope of Tianshan Mountain, which decreased from 42.73% in 1980 to 32.83% in 2030 with the expansion of 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² in 1980 to 3,392.38 km² 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⁸ 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⁸ 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⁸ 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⁸ 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⁸ yuan to 0.335 × 10⁸ yuan and food production by 3.00 % from 0.167 × 10⁸ yuan to 0.172 × 10⁸ 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 the ESV. For example, 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 (Zhang et al., 2020), 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.