Spatiotemporal changes, trade-offs, and synergistic relationships in ecosystem services provided by the Aral Sea Basin

Study area

The Aral Sea Basin (53°37′–78°21′, 33°48′–47°12′) covers 1.7 million square kilometers of Central Asia and is shared by Afghanistan, Iran, Turkmenistan, Kazakhstan, Tajikistan, Uzbekistan, and Kyrgyzstan (Deliry et al., 2020; Harriman, 2014; Fig. 1). The surface water resources in our study area entirely depended on the run-off from the Syr Darya and Amu Darya Rivers, which begin in the Pamir and Tian Shan Mountains (Lioubimtseva, 2015). Amu Darya, the most important river within the Aral Sea Basin, flows nearly 2,400 km from the mountains with an average annual flow of around 79 km³ from the drainage basin (Nezlin, Kostianoy & Lebedev, 2004). The total length of the Syr Darya River is 2,500 km and the average annual flow is 37 km³ (Nezlin, Kostianoy & Lebedev, 2004). The climate of the Aral Sea Basin is semiarid to arid with cold winters and hot summers (Nezlin, Kostianoy & Li, 2005). The annual precipitation in the region reaches 140 mm/year, and most of the region experiences a distinctive spring maximum (Deliry et al., 2020). Due to the obvious continental climate characteristics in this area, there is a large temperature difference between winter and summer. In 1960, the Aral Sea covered an area of 68,000 km², but it has since shrunk significantly because of the construction of canals from the Amu Darya and Syr Darya basins used to irrigate cultivated land. The potential evapotranspiration has increased in this region over recent decades, exacerbating the water shortage problem (Jiang et al., 2019).

Data source

Developments in earth observation technology have led to the use of a variety of land use datasets to study land use changes, but these current generated global land cover datasets have almost no consistency in their observation time periods, spatial resolution, or accuracy standards (Poulter et al., 2015). For example, the Copernicus Global Land Cover layers have a higher resolution (100 m), but only cover 2015–2019, and the Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Type Product (MCD12Q1) provides global land cover types at yearly intervals (2001–2016), but its resolution is only 500 m. In order to resolve these shortcomings, the European Space Agency (ESA) consistently maps moderate-resolution global land cover over time using surface reflectance from the Medium-Spectral Resolution Imaging Spectrometer (MERIS) aboard the Environmental Satellite (ENVISAT) and Systeme Probatoire d’Observation de la Terre (SPOT) 4 and 5. To validate the data, Bontemps et al. (2013) used the GlobCover 2009 verification database and found that the overall accuracy of the ESA’s Climate Change Initiative Land Cover (CCI-LC) maps was 75.4%. Reinhart et al. (2021) evaluated the consistency between datasets through proportional area comparison, and found that the overall accuracy of the ESA CCI-LC with Coordination and Information on the Environment (CORINE) Land Cover (CLC) reached ∼76%. Although the dataset divided the world into 37 land use types, to match Costanza’s biomes and consider the actual conditions present within the study area, we classified 29 land use types in the study area into seven major LULC types: cropland, forestland, grassland, wetland, urban, bare land, and water (including snow and ice) (Table S1).

Digital elevation data, obtained from the National Aeronautics and Space Administration (NASA) Shuttle Radar Topographic Mission (SRTM) at a spatial resolution of 90 m, were utilized to extract the stream network and contour lines of the Aral Sea Basin (Table 1). They were acquired with the same sensor in a single mission and produced with a single technique: Synthetic Aperture Radar (SAR) interferometry. The absolute and relative vertical accuracies of the Digital Elevation Model (DEM) were ±16 m and ±6 m, respectively (Rabus et al., 2003).

In addition, we also collected nine variables, including topographical, population, and land cover parameters, to simulate the contribution of various factors to land expansion (Table 1). The precipitation and temperature data for current climate conditions were extracted from the WorldClim global normal climate (the average for 1970–2000) database with a 30-arc-second (∼1 km) resolution (Fick & Hijmans, 2017). Population data were obtained from the WorldPop program with a resolution of 30-arc-seconds (∼1 km). The road information was extracted from the Socioeconomic Data and Applications Center (SEDAC) managed by the NASA Earth Science Data and Information System (ESDIS) project (Sherbinin et al., 2002). The water surface and slope data were extracted from the ESA CCI-LC maps and DEM, respectively. ArcGIS’s spatial analysis function was used to calculate the Euclidean distance between the water body and the roads, and all variables were resampled to a 1 km resolution. The data for all environmental variables were then resampled to a 1 km × 1 km resolution to match the land use simulation model.

Model projection of future LULC changes

In order to improve the geographic cellular automata model in the conversion rule mining strategy and landscape dynamic change simulation strategies, Liang et al. (2021) developed the patch-generating land use simulation (PLUS) model, which was based on a comprehensive analysis of the influence of man-made and natural factors.

The model proposes a rule mining framework based on the Land Expansion Analysis Strategy (LEAS), which comprehensively considers factors that affect land use changes, such as population, economy, traffic topography, and climate background, so that it can better simulate the LULC pattern with scattered patches (Fig. 2; Li et al., 2021).

The land use data from the two periods, 2015 and 2020, were overlaid together and the cells were extracted with the changed states from the later date, which represented the change regions for each land use type. Sampling points were randomly selected and divided into subsets according to their land use types, which were then analyzed separately using Random Forest Classification (RFC) algorithms to explore the relationship between the growth of each land use type and multiple driving factors (elevation, slope, temperature, precipitation, population, distance to water, distance to highway, distance to road, and distance to railway) (Fig. 3). On this basis, we proposed a cellular automata (CA) model based on multi-type Random Patch Seeds (CARS) to better simulate the fine-scale patch growth of various LULC types.

A multi-type random patch seeding mechanism based on threshold descent was used in the PLUS model in order to predict the patch evolution for various land-use type scenarios (Shi et al., 2021). By applying a Monte Carlo approach, the probability surface $P_{i,k}^{d=1}$ for each land-use type was determined when the neighborhood effect of land use k was 0, ${\displaystyle O{P}_{i,k}^{d=1,t}=\left\{\begin{array}{cc}{P}_{i,k}^{=1}\times \left(r\times {\mu }_k\right)\times D_k^t\hfill & \mathrm{if}{\Omega }_{i,k}^t=0\mathrm{and}r<P_{i,k}^{d=1}\hfill \\ P_{i,k}^{d=1}\times {\Omega }_{i,k}^t\times D_k^t\hfill & \text{all others}\hfill \end{array}\right.}$

where r is a random value ranging from 0 to 1; μ_k is the threshold used to generate the new land use patches for land use type k, which was determined by the model users. Seeds can create new forms of land use and expand into a series of new patches. If the new land-use type won in a round, a declining threshold τ was used to determine the candidate land-use type c chosen by the roulette wheel as follows: ${\displaystyle \begin{array}{cc}\hfill \hfill & \hfill \mathrm{If}\sum _{k=1}^N\left|G_c^{t-1}\right|-\sum _{k=1}^N\left|G_c^t\right|<Step\text{Then},l=l+1\hfill \\ \hfill \hfill & \hfill \left\{\begin{array}{cc}\text{Change}\hfill & P_{i,c}^{d=1}>\tau \mathrm{and}T{M}_{k,c}=1\hfill \\ \mathrm{No}\text{change}\hfill & P_{i,c}^{d=1}\le \tau \mathrm{or}T{M}_{k,c}=0\hfill \end{array}\right.\tau ={\delta }^l\times r1\hfill \end{array}}$

where Step is the step size of the PLUS model used to approximate the land use demand; δ is the decay factor of the decreasing threshold τ, which ranges from 0 to 1 and is set according to the expert; r1 is the normally distributed stochastic value with a mean of 1, ranging from 0 to 2; and l is the number of decay steps. TM_k,c is the transition matrix that defines whether land use type k is allowed to convert to type c (Verburg & Overmars, 2009). This study’s neighborhood size was 3. We set the declining threshold τ and decay factor δ as 0.7 and 0.2, respectively.

To ensure the reliability of the simulation results, we used the Kappa coefficient and Figure of Merit (FoM) to measure the goodness of fit for land use simulation. The calculation formula can be described as follows (Mitsova, Shuster & Wang, 2011): ${\displaystyle \text{kappa}=\frac{p_0-p_c}{1-p_c}}$

where kappa is the simulation accuracy index, p_c is the expected simulation accuracy in a random state, and p₀ is the actual simulation accuracy. ${\displaystyle \mathrm{FoM}=\mathrm{B}/\left(\mathrm{A}+\mathrm{B}+\mathrm{C}+\mathrm{D}\right)}$

where A refers to the error area in which the observed change is simulated as persistence, B refers to the correct area in which the observed change is simulated as change, C refers to the error area in which the observed change is simulated as change in the wrong classification, and D refers to the error area in which the observed persistence is simulated as change.

Previous research has demonstrated that the PLUS model can achieve good results in predicting land use changes (Liang et al., 2021) and that its simulation accuracy is better than that of the Future Land Use Simulation (FLUS) model, the Conversion of Land Use and its Effects (CLUE) model, and the Artificial Neural Networks–Cellular Automata (ANN-CA) (Li et al., 2021). In this study, we used the transition probability matrix from 2005 to 2010 to develop the 2015 LULC map (Table S2). The forecast data in 2015 was compared with the actual respective data to verify the model’s accuracy. The simulation performance evaluation showed that the kappa coefficient was 0.92 and the FoM was 0.24, which suggested that the simulation results were accurate and reliable. Then, we used the 2015 and 2020 maps and the FLUS model in the transition probability matrix to predict the 2025 LULC maps (Table S3).

Assessment of ESVs

In this study, nine ESVs, proposed by Xie et al. (2008), were calculated using the BTM. As the Aral Sea Basin is located in Central Asia, and 83.2% of the area is shared by Turkmenistan, Kazakhstan, Tajikistan, Uzbekistan and Kyrgyzstan, so the equivalent value coefficients of each ES referred to those proposed by Li et al. (2019; Table 2), which were derived from the value coefficient matrix proposed by Costanza et al. (Costanza et al., 2014). We used the assessment model proposed by Costanza et al. (Costanza et al., 1997) as a basis to evaluate the ESVs of the Aral Sea Basin. ${\displaystyle ES{V}_t=\sum A_k\times V_k}$ ${\displaystyle V_k=\sum _{i=1}^{n}ES{V}_{ki}}$

where ESV_t is the total ESVs, A_k is the area of the LULC type k, V_k is the ESV of LULC type k, and ESV_ki is the i type of ESV of LULC type k.

The changes in ESVs in the study area were calculated using the following equation derived from a study by Song & Deng (2017). ${\displaystyle C{R}_i=\frac{ES{V}_{end}-ES{V}_{start}}{ES{V}_{start}}\times 100\text{\%}}$

where CR_i is the change in the ESVs in grid i, ESV_end is the ESVs at the end of the study period for grid i, and ESV_start is the ESVs at the beginning of the study period.

Elasticity of ESV changes in relation to LULC

In economics, elasticity is a measure of how responsive one economic variable is to changes in another (Song & Deng, 2017). In this study, we aimed to investigate the responses of ESVs to LULC. We used a sensitivity analysis to evaluate changes in the ESVs in response to 50% ESV coefficient adjustments for each LULC type (Kindu et al., 2016). The formula for calculating the sensitivity coefficient refers to the standard economic concept of elasticity: ${\displaystyle \mathrm{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 CS is the coefficient of sensitivity; i and j represent the initial and adjusted values, respectively; and k is the LULC category. If CS was greater than one, then the estimated ecosystem value was considered elastic relative to that coefficient. However, if CS was less than one, then the estimated ecosystem value was considered inelastic and the results were reliable even if the value coefficient had a relatively low accuracy (Kindu et al., 2016; Kreuter et al., 2001).

Correlation analysis

Correlation analysis can be used to quantitatively describe the degree of linear correlation and clarify the direction of correlation between two variables (Liu, Wu & Chen, 2018). The larger the value, the stronger the correlation and the smaller the value, the weaker the correlation. When the result is greater than 0, this suggests that as one variable increases, the other also increases, making it a coordinated relationship. When the result is less than 0, there is a trade-off relationship. The formula is as follows: ${\displaystyle R_{xy}=\frac{\sum _{i=1}^n\left(x_i-\bar{x}\right)\left(y_i-\bar{y}\right)}{\sqrt{\sum _{i=1}^n{\left(x_i-\bar{x}\right)}^2}\sqrt{\sum _{i=1}^n{\left(y_i-\bar{y}\right)}^2}}}$

where R_xy is the correlation coefficient, n is the number of samples, x_iy_i are the ith value of x and y respectively, and $\bar{x}$ and $\bar{y}$ are the average values of x and y, respectively.

LULC change detection analysis

Figure 4 shows the spatial distribution of land use in the Aral Sea basin from 1995 to 2025. From 1995 to 2025, the largest cover type in the study area is bare land, accounting for approximately 45.3% of the total study area, followed by grassland and cropland, which account for 33.5% and 17.3%, respectively (Table 3). Among all types of land use, bare land increased the most. The total area increased by 111.44 × 10⁴ ha, and urban land increased by 546.54%. The average annual urban land notably increased by 18.2% from 1995 to 2025. The area of cultivated land increased by 15.2 × 10⁴ ha, but it tended to be stable from 2005 to 2015. It is worth noting that water bodies experienced the largest reduction in area, decreasing by 42.42% from 1995 to 2025, with an average annual reduction of 1.62 × 10⁴ ha. The forestland area decreased by 9.45 × 10⁴ ha. The wetland area has not significantly fluctuated in recent years.

Using the geospatial analysis software ArcGIS, combined with a digital elevation model map, we extracted the natural river network structure in the basin and contour lines within the basin with a distance of 100 m. Through our analysis, we found that the 1,000-meter contour line acted as the dividing line between the river and the alluvial plain of the mountain pass, and the dissipation area of the river was completely below the 100-meter contour line. Therefore, the 1,000-meter and 100-meter contour lines were used as the dividing lines of the upper, middle, and downstream parts of the river, and the land use changes in the upper, middle, and lower reaches were counted separately (Fig. 5).

The area of urban land in the upstream, midstream, and downstream watersheds is predicted to increase rapidly by 622.39%, 536.99%, and 608.98%, respectively (Table S4). Correspondingly, the area of cropland and grassland is also increasing. However, the forestland, wetland, bare land, and water bodies show obvious spatial differentiation. The forestland and water body areas in the upstream and downstream have shown an overall decreasing trend, while the middle regions are increasing. The wetland and bare land areas are constantly decreasing in the upstream area, but are increasing in the downstream reaches. According to the conversion relationship between the different land use types, the upstream forestland and water bodies have been transformed into grassland and cultivated land. In the middle reaches, excluding the shrinking bare land and wetlands, the other land use types have increased to varying degrees. The areas covered by water bodies in the downstream areas have been greatly reduced, and most have been converted into bare land. Grassland, urban land, cropland, and wetland areas have also expanded to varying degrees, while forestland has slightly decreased in area.

Changes in total ESVs

According to our estimation, the total ESVs of the Aral Sea basin in 1995 was approximately 492.55 billion USD (Table 4). Grasslands constituted the highest contribution at 48.08%, followed by cropland and water body areas (33.79% and 16.39%, respectively) (Fig. 6). Due to the observed LULC changes, the regional ESVs decreased by 16.55 billion USD from1995 to 2005, which can be mainly attributed to the decreased ESVs associated with water bodies that counteracted the ESV gains from cropland, grassland, and urban areas. The regional ESVs further decreased by 9.36 billion USD from 2005 to 2015. Similarly, the decrease in the water ESVs was the main reason for the decrease in the regional ESVs. However, the regional ESVs are predicted to increase by 3.03 billion USD from 2015 to 2025. This is mainly because the ESVs of the water bodies will no longer decline, and the ESVs of the grassland and urban areas show an upward trend, which will overcompensate for the ESV losses associated with cropland.

We further analyzed the ESV changes in the upper, middle, and lower reaches of the Aral Sea Basin (Table S5). From 1995 to 2025, the upstream ESVs are predicted to fluctuate the least, with an overall reduction of 1.9 billion USD. The ESVs in the midstream will continue to increase, with a total increase of 6.82 billion USD and, except for an observed slight decrease in the wetland ESVs, the ESVs of the other land use types will increase. Among the increased ESVs, the grassland ESV increase, reaching 43.47%, accounted for the largest proportion, followed by urban land and cultivated land at 34.41% and 15.03%, respectively (Table S5). The downstream ESVs showed a continuous decrease, with a total reduction of 27.79 billion USD. Despite the reduction in ESVs associated with water bodies (−29.79 billion USD) and forestland (−0.04 billion USD) areas in the downstream regions, the ESVs of other land use types increased to varying degrees.

Changes in ecosystem service function values

Table 5 shows the individual ecosystem function values (ESV_ki) of different land use types. The most important ecosystem functions in the Aral Sea Basin were biodiversity, food production, and water regulation, which contributed 35.9%, 28.12%, and 15.53%, respectively, of the total ESVs in 1995; 37.45%, 29.29%, and 12.78% in 2005; and 38.27%, 29.89%, and 11.2% in 2015. They are expected to contribute 38.21%, 29.75%, and 11.21%, respectively, of the total ESVs in 2025. The relative contributions of biodiversity and food production to regional ESVs have been increasing and overcompensated for losses in water regulation. Most of the ESV_ki were projected to increase between 1995 and 2025, except for water regulation, waste treatment, and culture and tourism, which decreased by 31.21%, 9.96%, and 5.17%, respectively. Across all ecosystem functions, water regulation had the highest change rate (Fig. 7).

Ecosystem sensitivity analysis

Although the wetland and bare land areas had the highest service value coefficients, the CS of the wetland and bare land areas was 0 because the wetland area was the smallest and the value of the bare land was zero. Meanwhile, the CS for the cropland increased from 0.34 to 0.36, and that of the grassland increased from 0.48 in 1995 to 0.51 in 2025 (Table 6). Compared with the grassland and cropland, the CS for the water body areas decreased from 0.16 to 0.1. In this study, all CS values were far smaller than “1”, indicating that the total estimated ecosystem value was inelastic with respect to the ecosystem value coefficients.

Analysis of the synergy degree of ecosystems trade-offs

After calculating the value of each ecosystem service in the Aral Sea Basin using the correlation analysis model, we obtained the correlations of nine types of ecosystem services (Table S6). When the result was positive, it indicated the presence of a synergistic relationship between two ecosystem services, i.e., two targeted ecosystem services had the same upward or downward trend in the same period of time, and an increase in one service promoted an increase in the other. When the result was negative, there was a trade-off relationship between the two ecosystem services, i.e., an increase in one ecosystem service led to a decrease in the other.

Among the 81 values representing the nine ecosystem services in the Aral Sea Basin, 36 were smaller than 0 and 45 were greater than 0 (Table S6), which suggested that a synergistic relationship was dominant (55.56%). All of the trade-off relationships in the Aral Sea Basin were related to water regulation, waste treatment and prescription, and culture and tourism, but there was a strong synergy among all of them (Fig. 8). Additionally, strong synergy (R > 0.9; p < 0.01) was identified pairwise between food production-raw material, climate regulation-gas regulation, soil formation retention-raw material, waste treatment-water regulation, biodiversity-food production, biodiversity-gas regulation, and biodiversity-climate regulation (Table S6). Significant negative correlations (R < −0.9; p < 0.01) were also found between pairwise water regulation-gas, water regulation-climate, waste treatment-gas, waste treatment-climate, biodiversity-water, and biodiversity-water regulation.

Impact of LULC change on ecosystem services in the Aral Sea Basin

Land use changes are caused and accelerated by many factors, including rapid urban expansion and cultivated land reclamation initiated by population expansion (Nahuelhual et al., 2014). In arid regions, water is the main factor affecting surface land use changes, and the use and distribution of water resources also directly lead to changes in ecosystem services (Kreuter et al., 2001). By analyzing the land use changes and calculating the ESVs of the Aral Sea Basin, we found significant urban area expansion (+599.22%) and the largest area increase from 1995–2025 can be attributed to grassland areas, which increased the ESVs of the Aral Sea Basin by 4.4 billion USD and 4.04 billion USD, respectively. Correspondingly, the area covered by water bodies was predicted to shrink the most (−39.9%), by 255.55 × 10⁴ ha from 1995 to 2025, causing a loss of 31.97 billion USD from the ESVs. The forestland area also shrank (0.09 × 10⁴ ha), resulting in an ESV loss of 2.73 million USD. Additionally, the cropland ecosystem service value increased by 0.68 billion USD from 1995 to 2025, which was mainly caused by the expansion of cropland areas (Table 4). Changes in the land use types of other areas were small.

To determine the partitioning and usage of water resources in the Aral Sea Basin, we calculated the changes in land use and ESVs in the upper, middle, and lower reaches of the Aral Sea Basin. The results showed that in the upper reaches of the Aral Sea Basin, where the LULC changed significantly, included urban land (+622.39%), waterbodies (−10.93%), and cropland (−1.83%) (Table S4). However, due to the small proportion in the total land area of the upper reaches, the ESVs of the upper reaches only slightly decreased (−1.9 billion USD) (Table S5). From 1995 to 2025, the ESVs of the middle reaches of the Aral Sea Basin will continue to increase. Except for a small decrease in the wetland ESVs, the ESVs of the other land use types increased to varying degrees. Urban, grassland, and cropland areas increased the most, by 3.78 billion USD, 1.79 billion USD, and 1.07 billion USD, respectively (Table S5). In this region, the bare land area decreased by 121.58 × 10⁴ ha, and the wetlands decreased by only 0.03 × 10⁴ ha. The fastest growth rates were found in urban (+536.99%) and forestlands (+24.57%) (Table S4). Unlike the upstream and midstream ESVs, the downstream ESVs experienced a drastic reduction. The ESVs of the water body areas experienced the largest reduction, accounting for 61.37% of the total ESVs reduction in the region, with a total reduction of 29.79 billion USD. However, the ESVs of the other land use types increased to varying degrees, including cropland (+0.46 billion USD), grassland (+1.18 billion USD), and urban land (+0.39 billion USD) (Table S5).

In comparison, the cropland, grassland, and urban land areas in the Aral Sea basin continued to increase, especially in the middle reaches, and accounted for 70.18%, 44.34%, and 86.02% of the total increase, respectively. With the continuous expansion of urban land, cropland, and grassland areas, water consumption continued to increase. The forestland and water body areas in the upstream and downstream reaches decreased, but they increased in the middle reaches. We believe that the midstream area intercepted a large amount of water from the Syr Darya and Amu Darya that allowed it to reclaim the bare land and develop cropland, grassland, and forestland areas through higher water consumption and building reservoirs, which has intensified the shrinkage of the Aral Sea Basin.

More cropland areas can produce more food to feed more people and bring considerable economic benefits. However, the Aral Sea Basin has an arid-semiarid climate, and the development of cropland consumes many water resources. Although more cropland brings considerable economic benefits, it can also result in the loss of natural ecosystem services (Li et al., 2019; Shi et al., 2020). In the Aral Sea Basin area, the original water body became bare land, which directly led to a rapid reduction in ESVs in this area. In the middle reaches of the region, the ESVs increased due to the expansion of cropland, grassland, and urban areas, but it could not compensate for the loss of downstream ESVs. At the same time, we projected that the expansion of agricultural land largely will lead to an increase in service functions related to food production (+1.2 billion USD), biodiversity (+2.5 billion USD), climate regulation (+0.7 billion USD), and soil formation (+0.1 billion USD) between 1995 and 2025 in the Aral Sea Basin (Table 5). However, the results show that agriculture and urban expansion have a negative impact on the management of other important ecosystem services, such as water regulation (−23.7 billion USD), waste treatment (−2.2 billion USD) and recreation, and culture and tourism (−1.6 billion USD).

Additionally, LULC changes are ecological processes dominated by temperature and precipitation, and climate change can also cause changes in land use types, which indirectly lead to changes in ESVs (Zhou et al., 2015). For example, the warming trends observed in Central Asia initially enhanced vegetative greenness before 1991, resulting in the enhancement of climate regulation in this region (Zhou et al., 2015). In recent years, however, precipitation in the Aral Sea basin has slightly decreased, the temperature has increased significantly, and the increased surface evaporation has put more pressure on the water resources in the region (Jiang et al., 2017), which will further affect LULC changes and the resulting ecosystem services.

Analyzing the cooperative relationship of ecosystem services

The Aral Sea Basin is located in an inland arid area with diverse landform types, scarce rivers, and a unique relationship between various ecosystems. According to our results, there is a dominant synergistic relationship between the various ecosystem services in the Aral Sea Basin. A trade-off relationship exists between the ecosystem services of water regulation, waste treatment and recreation, culture and tourism, among others. However, these three ecosystem services (water regulation, waste treatment and recreation, culture and tourism) are synergistic in pairs.

In arid areas, water is a scarce resource that has unparalleled advantages in cultural tourism services. Additionally, water has a purification function and water in the Aral Sea Basin dominates almost all ecological processes. Since almost all of the nine ecosystem services (except for waste treatment and recreation, and culture and tourism) require water, they all have a trade-off relationship with water regulation. There are similar reasons behind the exceptions of waste treatment and recreation, and culture and tourism.

Because most of the land in the Aral Sea Basin is bare desert, increasing food production can increase the coverage of surface vegetation. Therefore, there is a synergistic relationship between raw material production, gas regulation, soil-formation retention, and biodiversity. The raw material designation generally refers to areas with a higher proportion of surface area covered by forests with strong soil retention capabilities. Surface vegetation cover has a strong effect on climate regulation, especially in arid areas. Areas with strong climate regulation functions also have high vegetation coverage and diverse plant species, so there is a strong synergistic relationship between climate regulation and biodiversity. In the same way, we also observed a strong synergy between climatic and atmospheric regulation and biodiversity.

Protection policies for the Aral Sea basin

The reduction and irrational distribution of water resources in the Aral Sea Basin are the main factors decreasing the regional ESVs. Therefore, the protection and suitable use of water resources should be the focus of conservation efforts. Increases in population growth and social development will encourage urban land areas to continue to expand, and urban expansion will take up a large amount of arable land. To maintain agricultural production, more arable land can be developed from bare land and grassland areas, which consume more water resources. However, the Aral Sea Basin encompasses seven countries, and the coordinated use of water resources is particularly important. We suggest that the Aral Sea Conservation Organization distribute water resources from the Syr Darya and Amu Darya basins in a fair and reasonable way in order to prevent overexploitation.

Additionally, environmental protection policies for each country in this area should mainly focus on coordinating the relationship between cultivated and urban land areas, encouraging the return of unfavorable cultivated land to grassland, and preventing the further conversion of grassland or bare land to cultivated land in order to prevent the consumption of more water. First, permanently cultivated land should be designated as such in order to avoid the encroachment of built-up areas. This designation should not change under any circumstances, especially in areas with fertile soils. Otherwise, more arable land would need to be cultivated to produce the same yield, which would require more water consumption. Second, greater fees should be imposed for converting cultivated land into built-up areas. These fees should be used to subsidize farmers so that they may install more drip irrigation systems. Finally, to encourage the integration of rural farmland, certain types of crops should be planted according to their distance from a water source. Plants that consume a large amount of water, such as cotton and rice, should be cultivated on arable land that is close to rivers in order to prevent water loss during transportation.

Limitations and areas of future research

In this study, the calculation of ESVs depended on changes in land use types, and the accuracy of determining land use types using remote sensing methods was affected by many factors, such as classification methods and data resolution (Li et al., 2019). These factors may be inaccurate during land use classification, which may lead to the underestimation or overestimation of ESVs (Sexton et al., 2015). During the evaluation process, we divided 37 land use types into seven categories (Table S1), among which shrubs were classified as grasslands, which may have caused an overestimation of ESVs.

When simulating land use, we chose the newly-proposed PLUS model, which takes into account many factors affecting land use distribution and uses Random Forest learning algorithms to analyze the impact of different influencing factors on land use. The reliability of the model was verified. As indicated by the results of the model validation and comparison, the simulation results obtained higher simulation accuracy (Kappa: 0.92; FoM: 0.24) than other land use simulation models, which was similar to the results of Liang et al. (2021). However, the PLUS model cannot simulate the evolution of land use composition and mixed structure at the sub-cell scale. Additionally, the premise of our prediction was that the region’s land use policies will not change (Zhang et al., 2020). These factors and assumptions may have caused deviations between the forecast results and the actual future conditions (López-Marrero et al., 2011).

Indeed, when addressing land-use conversion, we only assessed the responses of ecosystem services to changes in the area constituting a given land-use type. There was no discussion of the micromechanisms behind any changes in ecosystem services as part of the land use conversion process (Song & Deng, 2017), and no analysis of the socioeconomic response to land use changes was performed. Water from the Aral Sea is used to reclaim arable land to irrigate crops (considered a reasonable use), which can produce more food to feed more people. Future studies should explore which approach provides more social and economic benefits when protecting the Aral Sea.