Variations in ecosystem service value in response to land use/land cover changes in Central Asia from 1995–2035

Study area

Our study area covered five Central Asian countries (from 35°08′N, 55°25′N to 46°28′E, 87°29′E), including Tajikistan, Kyrgyzstan, Turkmenistan, Uzbekistan and Kazakhstan (Fig. 1). With a land area of 4 million km², Central Asia extends from Russia in the north to Afghanistan in the south and from the Caspian Sea in the west to western borders of China in the east (Ozturk et al., 2017). The terrain gradually changes from the western Caspian lowlands to the Altai Mountains, Tianshan Mountains and Pamirs (Beurs et al., 2015). Because this area is located far from the ocean in the hinterland of the Eurasian continent, it has a distinctive continental arid and semiarid climate with low precipitation and intensive evaporation (Lioubimtseva & Henebry, 2009). In recent decades, potential evapotranspiration has increased in Central Asia, especially in the Aral Sea region and western Kazakhstan, with an annual increase of 7.42 mm/year (Jiang et al., 2019). Central Asia has a heterogeneous landscape, including diverse land-cover types such as temperate deserts (e.g., Kyzylkum, Karakum Desert), forests, lakes (e.g., the Aral Sea, Balkhash Lake) and vast grasslands. Major transboundary rivers, such as the Syr Darya, the Amu Darya, the Irtysh River, and the Ili River, are critical water sources for regional ecosystems and agriculture (Zhou et al., 2019). With the large-scale development of irrigated agriculture in the Aral Sea Basin in the early 1960s, the irrigated cultivated land area increased by 60%, and the planted cotton area doubled between the 1960s and the 1990s. Now, irrigated land accounts for almost half of the total cultivated land area (Suo et al., 2009).

Data collection and land use classification

Multiscale satellite observations were used to characterize LULC change (Song, 2018). In a preliminary study, we compared the multiple LULC datasets for the study area, including the GLC2000 (Corresponding & Belward, 2005), the GlobCover 2009 (Arino et al., 2008), MODIS land cover (Friedl et al., 2010), and the newly released annual European Space Agency Climate Change Initiative land cover maps (CCI-LC) (http://maps.elie.ucl.ac.be/CCI/viewer) (Radoux et al., 2013). The comparison showed that the CCI-LC has the highest spatial resolution and has better accuracy in the study area (Hartley et al., 2017; Wei et al., 2015). This product was developed using the GlobCover unsupervised classification chain and combined a variety of earth observation products based on ESA’s GlobCover products. Unlike many remote sensing products based on single-sensor methods, this dataset is generated using multiple sensors, such as the Advanced Very High Resolution Radiometer (AVHRR), System Probatoire d’Observation de la Terre Vegetation (SPOT-VGT), and PROBA-V. To validate the data, Defourny et al. (2016) compared the CCI LC product for the 2010 with the certain and homogeneous points of the GlobCover 2009, with a total accuracy of 73.2% (Georgievski & Hagemann, 2019). More high-accuracy land classifications were found: water bodies (92%–96%), irrigated cropland (89%–92%), rainfed irrigated cropland (89%–83%), bare land (89%–88%) and urban areas (88%–86%) (Defourny et al., 2016). However, natural vegetation and sparse vegetation have low user accuracy values, but all errors relate to confusion between these and other LULC classes, which very much limits the impact of the low values. Details about the CCI-LC dataset, including its accuracy and the confusion matrix, can be found in the study by Defourny et al. (2016). To match Costanza’s biomes, plant functional types were classified into seven major LULC types, including cropland, forestland, grassland, wetland, urban areas, bare land and water bodies (Table S1).

Model projection of future LULC change

Future LULC change from 2015 to 2035 was projected by combining a Geographic Information System (GIS) with a CA-Markov model, which is a robust approach for simulation of the spatial—temporal change patterns of LULC (e.g., Fu, Wang & Yang, 2018; Hartley et al., 2017; Muller & Middleton, 1994). The CA-Markov model merges cellular automata (CA) and Markov chain (Ji et al., 2009). The Markov chain is constructed based on the probability of change matrices for LUCC from t to t + 1: ${\displaystyle s_{\left(t+1\right)}=s_{\left(t\right)}\times A}$

where $s_{\left(t+1\right)}$ is the state probability of any time, and $s_{\left(t\right)}$ is the initial state probability. A is the transition probability matrix, and the formula is as follows: ${\displaystyle A=\left[\begin{array}{cccc}\hfill A_{11}\hfill & \hfill A_{12}\hfill & \hfill \dots \hfill & \hfill A_{1n}\hfill \\ \hfill A_{21}\hfill & \hfill A_{22}\hfill & \hfill \dots \hfill & \hfill A_{2n}\hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill A_{n1}\hfill & \hfill A_{n2}\hfill & \hfill \dots \hfill & \hfill A_{nn}\hfill \end{array}\right]}$

where A_ij is the sum of areas from the ith land cover category to the jth category from the initial period to forecast periods, and n is the number of LULC categories. We accomplished this process by using the MARKOV module in the IDRISI, a software that integrates GIS and image processing functions.

The CA model is used to predict the space–time dynamic change pattern using a transition map of the LULC (Yang, Zheng & Lv, 2012), and the model can be defined as follows (Wang, Zheng & Zang, 2012): ${\displaystyle 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): ${\displaystyle kappa=\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.

By incorporating the advantages of these two methods, the CA–Markov model is able to accomplish a better simulation of LULC changes both in quantity and space (Yang et al., 2014). In this study, the transition probability matrix was performed for the time period between 1995 and 2005 to predict the LULC map of 2015, which would be used to verify model accuracy. Then, using the 2015 classified map as the LULC baseline and the 2005 and 2015 maps for the transition probability matrix, we predicted the 2025 and 2035 LULC maps with the CA-Markov model.

Assessment of ecosystem service values

In this study, ESV were estimated based on benefit transfer method (BTM) proposed by Costanza et al. (1997). Nine ecosystem service functions generated from Xie et al. (2008), derived from the 17 ecosystem services listed by Costanza et al. (1997), were selected in this assessment (Sannigrahi et al., 2018). Cultivated land, temperate forest, grassland/shrubland, swamps/floodplains, urban areas, and desert/tundra/ice/bare-rocks in Central Asia were matched to cropland, forests, grassland, wetland, urban and bare land in Costanza et al.’s (2014) model, respectively (Table S1). The biomes that we used as proxies for the 7 LULC categories are not perfectly matched to those of Costanza et al.’s (1997) ESV model in some cases (Table S1), but they are closely related (Kreuter et al., 2001). The equivalent value coefficient of each ES updated by Costanza et al. (2014) was used to calculate ecosystem services values (Table 1). The formula is as follows: ${\displaystyle ES{V}_k={\Sigma }_f{A}_k\times V{C}_{kf}}$ ${\displaystyle ES{V}_f={\Sigma }_k{A}_k\times V{C}_{kf}}$ ${\displaystyle ESV={\Sigma }_f{\Sigma }_k{A}_k\times V{C}_{kf}}$ where ESV_k refers to the ecosystem service value of LULC type ‘k’, A_k represents the area (ha) for LULC category ‘k’, and VC_fk is the value coefficient (US $/ha/year) of function  f for the LULC type ‘k’ (Kreuter et al., 2001). ESV_f is the ecosystem service value of service function f, and ESV is the total ecosystem service value. We use the following formula to evaluate changes in ESV: ${\displaystyle ES{V}_{cr}=\frac{ES{V}_{t2}-ES{V}_{t1}}{ES{V}_{t1}}\times 100\text{\%}.}$

In this expression, ESV_cr refers to the change rate of ESV from the initial year to the final year, and ESV_t1 and ESV_t2 represent the total ESV at the start and end of the study period, respectively.

Elasticity for the response of ESV to LULC change

The biomes that we used as proxies for the 7 LULC categories are not perfectly matched to Costanza et al.’s (1997) ESV model in some cases (Table S1), which results in uncertainties in the assessment of the ESV. Thus, we used sensitivity analysis to evaluate the changes in ESV in response to 50% adjustments of the ESV coefficients for each LULC type (Kindu et al., 2016). The standard economic concept of elasticity was used to calculate the coefficient of sensitivity (CS) using the following formula (Kreuter et al., 2001): ${\displaystyle CS=\frac{\left(ES{V}_j-ES{V}_i\right)∕ES{V}_i}{\left(ES{V}_{jk}-ES{V}_{ik}\right)∕ES{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 (Kreuter et al., 2001).

Analysis of LULC dynamics

Combining the GIS technology with the CA-Markov model, we used the LULC base map from 2005 and transition probabilities from 1995 to 2005 to simulate the LULC for 2015. The kappa statistic of 0.93 indicates that there is good consistency with the actual value of the LULC types and the predicted results for the base year. Then, the future LULC in 2025 was predicted with the CA-Markov model in IDRISI using the LULC base map from 2015 and the transition probabilities from 2005 to 2015. Following the above process, the future LULC in 2035 was predicted.

Figure 2 presents the patterns of spatial distribution of LULC in Central Asia from 1995 to 2035, and Table 2 shows the magnitude of changes for the same periods. In 1995, grassland occupied approximately 51.39% of the study area, followed by bare land and cropland, which occupied 23.90% and 18.98% of the study area, respectively (Fig. 3). During 1995–2015, cropland and urban area increased substantially. The cropland expanded at a rate of 0.76% per annum, increasing 1,219.56  × 10⁴ ha by the year 2015 (Table 2). Rapid urbanization also caused the proportion of urban build-up to increase from 27.57 × 10⁴ ha in 1995 to 60.21  × 10⁴ ha in 2005 and then to 89.19  × 10⁴ ha in 2015, with an average growth rate of 10.64% per year (Table 2). Cropland and urban areas were expected to continue to increase in the 2025 and 2035 periods. However, the coverage of bare land, mainly concentrated in the Kyzylkum Desert and Karakum Desert, decreased from 23.90% to 22.95% during 1995–2015 and was projected to further decrease to 21.65% by 2035.

Wetlands and water bodies only account for 3% of the total study area. Wetlands increased from 121.14 ×10⁴ ha in 1995 to 125.25 ×10⁴ ha in 2015 and 128. 41 ×10⁴ ha in 2035 (Table 2), mainly near Balkhash Lake (Fig. 2). The area of water bodies decreased 292.59  × 10⁴ ha from 1995 to 2015 and is expected to further decrease 221.90 × 10⁴ ha by 2035. Among water bodies in Central Asia, the Aral Sea has shrunk most dramatically (Fig. 2).

Changes in total ecosystem services values

According to our estimation, the total ESV in Central Asia was approximately 1,505.31 billion US$ in 1995 (Table 3). Grassland had the highest contribution of 56.90%, followed by cropland and water bodies (28.09% and 11.15%, respectively) (Fig. 3). Due to LULC change, the regional ESV increased 5.68 billion US$ during 1995–2005, mainly due to the increased ESV in cropland and urban build-up, which overcompensated the ESV loss in grassland and water bodies. The regional ESV further increased 5.23 billion US$ from 2005–2015. Overall, the ESV in Central Asia increased 10.91 billion US$ during 1995–2015. It is noteworthy that the proportion of water bodies decreased sharply by 21.80% from 1995–2015, causing a loss of 36.61 billion US$. These trends will continue to occur in 2025 and 2035 (Table 3).

We further analyzed the ESV for administrative units in 1995 (Fig. S1). The highest ESV was found in the Karaganda state (174.09 billion US$), followed by East Kazakhstan (131.69 billion US$), Aktobe (121.09 billion US$) and Almaty (116.70 billion US$) states. The ESV for these administrative units were mainly contributed by grassland in Karaganda (86.13%), East Kazakhstan (67.14.8%) and Aktobe (86.39%). Andijon had the lowest ESV of 2.17 billion US$, 93.73% of which was contributed by cropland. In addition, we calculated the ESV change rates during 1995–2015, 2005–2015, 2015–2025 and 2025–2035 for the administrative units (Fig. 4). From 1995–2005, the ESV of Karakalpakistan, Uzbekistan declined 31.08%, with a total loss of 17.36 billion US$ (Fig. 4A). The ESV in Kyzylorda also decreased substantially, which was mainly caused by shrinking of water bodies. In contrast, the ESV in Pavlodar and Karaganda increased substantially, mainly due to the increased cropland ESV. During 2005–2015, the ESV in Karakalpakistan, Kyzylorda and Manghystau further decreased (Fig. 4B), with the highest ESV loss found in Karakalpakistan (−24.77%), followed by Kyzylorda (−8.83%) and Manghystau (−1.05%). In contrast, Karaganda had the highest increase rate of ESV (5.72%), followed by Aktobe (1.41%) and Almaty (1.25%). It is worth noting that the trend of decreasing ESV will continue to occur in Karakalpakistan and Kyzylorda from 2015 to 2035 (Figs. 4C and 4D).

Changes in values of ecosystem service functions

Table 4 shows the changes in individual ecosystem functions (ESV_f). The most important ESV_f in Central Asia were biodiversity, food production and water regulation, which contributed to 40.44%, 28.30% and 11.78% of the total ESV in 1995, respectively, 40.03%, 29.47% and 10.21% of the total ESV in 2015, and 40.51%, 30.14% and 8.93% of the total ESV in 2035. Most of the ESV_f decreased during 1995–2015 except for food production, raw materials, climate regulation, soil formation and waste treatment, which increased by 4.87%, 7.92%, 12.11%, 12.01%, and 2.91%, respectively (Fig. 5). It is noteworthy that the ESV of water regulation declined more rapidly than other ecosystem services (−12.70%), followed by gas regulation (−3.00%), cultural and tourism (−3.14%) and biodiversity (−0.29%). However, most of the ESV_f were projected to increase from 2015–2035 (Fig. 5). Only the ESV of water regulation and cultural service/tourism were expected to decrease in the future.

Ecosystem sensitivity analysis

In the observed (1995–2015) and projected (2025–2035) study periods, CS for grassland was the highest (0.55) due to the high service value coefficient and large grassland area (Table 5). Meanwhile, CS for cropland increased from 0.28 in 1995 to 0.32 in 2015 and to 0.34 in 2035. Compared with grassland and cropland, the CS (0.02) for forestland was relatively stable. The CS for water bodies decreased from 0.11 in 1995 to 0.09 in 2015 and to 0.07 in 2035. In this study, all CS were far less than “1”, indicating that the total estimated ecosystem values are inelastic with respect to the ecosystem value coefficients.

Impact of LULC change on ecosystem services in Central Asia

Globally, patterns of LULC change are themselves changing rapidly due to the acceleration of processes such as population growth, expansion of urban areas, and agricultural intensification (Lambin & Geist, 2006), which could affect regional ecosystem services (Hu, Liu & Min, 2008; Nahuelhual et al., 2014); these effects are particularly prominent in the Central Asia region. The United Nations statistics show that Central Asia’s urban population increased by 37.97% between 1995 and 2015, from 24.01 million to 33.13 million, which led to the overexploitation of water and land resources to meet people’s needs for water, food, and energy (Granit et al., 2012). We found significant expansions of cropland (+22.10%) and urban areas (+322.40%) and shrinking of water bodies (−38.43%) and bare land (−9.42%) during 1995–2035 (Table 2). Correspondingly, cropland ecosystem services value increased by 93.45 billion US$ from 1995 to 2035, which was mainly caused by the expansion of cropland area (Table 3). However, the area of water bodies decreased sharply during 1995–2035, causing a loss of 64.38 billion US$ (Table 3).

Cropland changes may seem economically profitable, but a large increase in agricultural land can result in the loss of natural ecosystem services. More specifically, the impact of complex human agricultural activities on the ecosystem services of water quantity and quality is considerably overlooked by current agricultural management models that maximize food production and minimize other ecosystem services (Coupe, Barlow & Capel, 2012). We found that the expansion of agricultural land has largely led to an increase in the service functions of food production (+32.23 billion US$), raw materials (+32.23 billion US$), soil formation (+10.86 billion US$) and biodiversity (+9.72 billion US$) during 1995–2035 in Central Asia (Table 4). However, these findings are consistent with the results of numerous studies around the world showing that agriculture and urban expansion have a negative impact on the provision of other important ES, such as water regulation (−41.03 billion US$) (Schröter et al., 2005), gas regulation (−0.04 billion US$) (Wang et al., 2006), and recreation, cultural and tourism (−21.91 billion US$) (Nahuelhual et al., 2014). Especially in Amu River Basin, irrigation water consumption of crops such as wheat, cotton and maize increased by at least 60% from 1962 to 2002, leading to a significant reduction in the amount of water flowing into the Aral Sea (McDermid & Winter, 2017). Now the Aral Sea has shrunk to approximately 1/10 of its original area (Ablekim et al., 2017). Our results show that the ESV in Karakalpakistan, Uzbekistan, decreased by more than 50% during 1995–2035 (Fig. 4), which was mainly caused by the shrinkage of the Aral Sea area. Therefore, wider cultivation of crops with lower water consumption and higher added value would reduce the shrinkage of the Aral Sea and thereby enhance ecosystem services.

In addition, LULC changes are interrelated with other global processes such as climate change and land degradation (Kertész, Nagy & Balázs, 2019), which directly or indirectly affect local ecosystem services (Hossu et al., 2019; Qiao et al., 2019). For example, approximately 23% of farmland in the Amu Lowland has been degraded due to inappropriate land use activities (Djanibekov et al., 2018), resulting in the loss of biodiversity and the decline of carbon storage capacity, which will further affect agricultural production, disturbance regulation and climate regulation (Sutton et al., 2016). Namely, driven by common factors (e.g., land use changes and land degradation, etc.), synergies and trade-offs exist between system services between ESs over time (Cord et al., 2017).

Limitations and areas of future research

Remote sensing images are the most important data source for research in ES, and LULC is the most widely used variable for assessment in ESV (Song, 2018). However, the limitations of global land cover data arise from the product generation process, including satellite sensor characteristics (spectral, temporal and spatial resolutions), definition and classification methods of land cover. Because of error propagation in the process of quantifying land cover data (Sexton et al., 2015), LULC change may be considerably underestimated or overestimated (Table 2). Such uncertainties are inevitably introduced in the analysis of spatial change patterns when the biomes used as proxies match LULC types (Kreuter et al., 2001). In this study, shrubs and sparse vegetation are used as grasslands (Table S1), which obviously exaggerates the area of grasslands and leads to an overestimation of ESV (Table 3). Therefore, to address these limitations, higher spatial resolution remote sensing data and more precise LULC classification will be used to further accurately assess ecosystem services in Central Asia.

Nine ecosystem service functions proposed by Xie et al. (2008) were selected to calculate the total ESV in Central Asia (Table 1), which could mask the value of other key ecosystem service functions. For example, sand fixation is an important function provided by desert ecosystems in arid and semi-arid areas, which has a close relationship with ecological benefits (Li & Xu, 2019; Musa, Deming & Cunyang, 2014; Taylor et al., 2017). Hao et al. (2013) found that the ESV of sand fixation in the Ulan Buh desert ecosystem was 0.68 billion US$ in 2009. However, there is no comprehensive assessment of ecosystem service functions for dryland ecosystem services in this study. To address this limitations, Schild et al. (2018) have identified and compiled a comprehensive database to estimate the ESV in arid areas. However, their findings indicate that there are still important limitations in the valuation of ESV in arid region. In further research, we will critically evaluate and further improve monetary valuation techniques to promote sustainable land management in Central Asia.

Our approach for estimating the ecosystem services value by using simple benefits transfer methods has certain limitations (Schmidt, Manceur & Seppelt, 2016). For example, we assumed homogeneity of ESV in the entire LULC category and generalized the unit values of one area as the average unit value of all areas. As Costanza et al. (2014) estimated, farmland ecosystem services provided $ 5567/ha in 2011 (Table 1). Although the average unit value provides a common value coefficient for all nations of the world to compare the relative differences, different types of farmland (e.g., rainfed cropland in northern Kazakhstan and irrigate cropland in Amu River Delta) have different crops planting structures and provide different ecosystem services and functions (Zheng, Zhang & Cao, 2018). In addition, the costs of ecosystem services (e.g., soil salinization, the loss of genetic resources, and water eutrophication) were ignored in assessing the value of ecosystem services, which might exaggerate the natural ecosystem services value (Cao et al., 2018). As Cao et al. (2018) revealed, the value of ecosystem services decreased by 52.66% in 2014 when the costs of water resources, investments in ecological protection, land rent and management were considered to estimate the net farmland ecosystem services value in China. In Central Asia, especially in Uzbekistan, irrigated agriculture remains the basis of the Uzbek economy, consuming more than 90% of the water resources and 228.40 kg/ha of fertilizers (FAO, 2015) and employing 33% of the labor force (UZSTAT, 2010). All of these factors suggest that farmland ecosystems have paid massive costs while providing a large number of services. However, because we did not incorporate these negative factors into the evaluation system, the ESV of farmland was overestimated (Table 3). Therefore, development of a framework that can assess the net value of ecosystem services is urgently needed to make reasonable arrangements for resource allocation.

We determined that the CS for the ESV was quite low (Table 5) when the value coefficients were adjusted by 50% to estimate the effects. However, we only focused on the directional change in ESV but ignored the magnitude of the ecosystem services values at specific points in time. It is worth noting that the interaction of some ES changed dramatically with different levels of economic development. For example, ESV of urban areas was overestimated (Table 3) in this study, masking the loss of essential ES provided by other LULC types. Yi et al. (2017) also believed that the value coefficient of urban built-up was assigned Costanza et al. (2014) seems unreasonable. We suggest that the value coefficients should be adjusted by collecting satellite or radar imagery and products obtained through efficient processing (e.g., soil moisture content, vegetation carbon stocks and chlorophyll content) to accurately reflect local ecosystem services of Central Asia in further studies.