Modeling and analyzing supply-demand relationships of water resources in Xinjiang from a perspective of ecosystem services

LI Feng, LI Yaoming, ZHOU Xuewen, YIN Zun, LIU Tie, XIN Qinchuan*

1 State Key Laboratory of Desert and Oasis Ecology, Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences,Urumqi 830011, China;

2 Research Center for Ecology and Environment of Central Asia, Chinese Academy of Sciences, Urumqi 830011, China;

3 University of Chinese Academy of Sciences, Beijing 100049, China;

4 School of Atmospheric Sciences, Sun Yat-sen University, Guangzhou 510275, China;

5 Geophysical Fluid Dynamics Laboratory, Princeton University, Princeton, New Jersey, NJ 08542, USA

Abstract: Water shortage is one bottleneck that limits economic and social developments in arid and semi-arid areas. As the impacts of climate change and human disturbance intensify across time,uncertainties in both water resource supplies and demands increase in arid and semi-arid areas. Taking a typical arid region in China, Xinjiang Uygur Autonomous Region, as an example, water yield depth (WYD)and water utilization depth (WUD) from 2002 to 2018 were simulated using the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model and socioeconomic data. The supply-demand relationships of water resources were analyzed using the ecosystem service indices including water supply-demand difference (WSDD) and water supply rate (WSR). The internal factors in changes of WYD and WUD were explored using the controlled variable method. The results show that the supplydemand relationships of water resources in Xinjiang were in a slight deficit, but the deficit was alleviated due to increased precipitation and decreased WUD of irrigation. WYD generally experienced an increasing trend, and significant increase mainly occurred in the oasis areas surrounding both the Junggar Basin and Tarim Basin. WUD had a downward trend with a decline of 20.70%, especially in oasis areas.Water resources in most areas of Xinjiang were fully utilized and the utilization efficiency of water resources increased. The water yield module in the InVEST model was calibrated and validated using gauging station data in Xinjiang, and the result shows that the use of satellite-based water storage data helped to decrease the bias error of the InVEST model by 0.69×108 m3. This study analyzed water resource supplies and demands from a perspective of ecosystem services, which expanded the scope of the application of ecosystem services and increased the research perspective of water resource evaluation.The results could provide guidance for water resource management such as spatial allocation and structural optimization of water resources in arid and semi-arid areas.

Keywords: ecosystem services; water resources; climate change; human activities; arid and semi-arid areas; InVEST model; Xinjiang

1 Introduction

Ecosystem services refer to various benefits that humans obtain from natural ecosystems directly or indirectly (Millennium Ecosystem Assessment, 2005). According to the report of Millennium Ecosystem Assessment, ecosystem services include a broad range of services such as supply,regulation, cultural, and support services. Water resource supplies belong to the supply services and have an important role in the evaluation systems of ecosystem services (Costanza et al., 1998;Wang et al., 2019; Ouyang et al., 2020). The number of studies on water resource supplies ranked top among supply services in researches on global watershed ecosystem services (Kaval, 2019)and the number of studies on water yield ranked the first among 130 articles on ecosystem services, as summarized by Ochoa and Urbina-Cardona (2017).

Water supplies are particularly important for arid and semi-arid areas and limit economic and social developments in these areas. Arid and semi-arid areas account for approximately 46.00% of the global land area, where nearly one-third of the global population live in (Lal, 2019). Water shortage brings ecological and social problems to arid and semi-arid areas, including drinking water reduction, food shortage, frequent sand storms, expanded desertification, land salinization,and biodiversity loss (Yu et al., 2019). Both climate change and human activities could intensify uncertainties in water resource supplies and demands in arid and semi-arid areas. The unique land surface environment and radiative transfer processes in arid and semi-arid areas cause temperature to increase more rapidly in these areas than in other regions, leading to frequent heat waves, increased evaporation, and glacier melting, as documented by studies using climate models (Huang et al., 2019; Yu et al., 2019). A series of unsustainable development practices,such as excessive cultivation, population overload, and rapid industrialization and urbanization,also contribute uncertainties to water resource demands.

Assessing the sustainable use of water resources in arid and semi-arid areas has become a hot topic in the context of climate change and human activities. The ecological footprint model, the fuzzy evaluation method, and the system dynamics model have been widely used in assessments on the sustainable use of water resources. The ecological footprint model originally proposed by William and Wackernagel (William, 1992; Wackernage and Rees, 1997; Seidl and Tisdell, 1999)was further evolved to the ecological footprint model of water resources, which was used to evaluate the sustainable use of water resources (Li et al., 2012). Novoa et al. (2019) used the ecological footprint model of water resources to assess the water footprint of agricultural production in a watershed in Chile and simulated the impact of inter-annual climate change on the water footprint. Jia et al. (2018) applied the fuzzy evaluation method to assess the carrying capacity of water resources for the Heihe River Basin in an arid region of Northwest China.Dawadi and Ahmad (2013) used a system dynamics method to simulate the changing process of water resources in the Las Vegas Valley in the semi-arid region of southern Nevada, USA and evaluated the impacts of climate change, population growth, and conservation policies on the water resource balances.

The supply-demand relationships in ecosystem services that reflect the coupling correlations between natural ecosystems and human society are now the research frontier in the field of ecosystem services (Yahdjian et al., 2015). Introducing the supply-demand relationships of ecosystem services into the assessments on the sustainable use of water resources is meaningful as compared with the abovementioned approaches. First, it enriches the assessment methods on sustainable use of water resources and likely circumvents problems in previous methods. For example, some concepts in the ecological footprint model of water resources, such as water ecological footprint and water carrying capacity, are artificially constructed, making it difficult to evaluate the accuracies of the models. By comparison, the water yield depth (WYD) results of the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model can be directly validated using gauging station data. Second, it expands the research direction of ecosystem services. Scholars have carried out researches on ecosystem services in arid lands, mostly focusing on ecosystem service tradeoffs (Kang et al., 2020), human well-being (Wei et al., 2018), and ecosystem service flow (Chen et al., 2020). However, there is also a need to assess and understand the sustainable use of water resources in arid and semi-arid areas based on ecosystem services.

Xinjiang Uygur Autonomous Region is the largest region in China and covers large arid and semi-arid areas in inland Asia. Accompanying with intensified climate change and human activities, uncertainties regarding both water supplies and demands affect water uses in Xinjiang and could exert considerable impacts on the sustainable development of the community and society. In recent years, the climate in Xinjiang has been warmer and wetter, with temperature increases at the rate of 0.34°C/10a, much higher than the global average (Shi et al., 2007; Chen et al., 2015; Huang et al., 2019). The changing climate in Xinjiang has also led to changes in land surface processes, such as increasing glacier melting (Piao et al., 2010) and evapotranspiration(Chen et al., 2015), and hence caused uncertainties to water supplies. In addition, Xinjiang has experienced rapid socioeconomic development since 2000. In this context, some scholars have carried out researches on supplies or demands of water resources in Xinjiang. In Xinjiang,increased water yield has been modeled and observed in Altay Mountains (Fu et al., 2017),Bosten Lake Basin (Yang et al., 2020), and nine high-alpine catchments (Luo et al., 2019). These studies revealed that increased precipitation and glacial meltwater both contributed to the increase in water yield. On the contrary, decreased water yield has been modeled and observed in the Manas River Basin (Xu et al., 2020) and Tarim River Basin (Wang et al., 2020). These studies explored that human activities may have led to the reduction in water yield. As for water utilization, it was considered in some studies that water utilization had an increasing trend in Xinjiang due to urbanization and industrialization (Lei et al., 2012; Zhang et al., 2020), but there are still some other studies that revealed that water-saving agricultural and industrial water recycling may led to the decrease in water utilization (Chen and Shi, 2016; Zhou et al., 2020). In the context of climate change and human activities, how do the changes of the supply-demand relationships of water resources in Xinjiang? This question needs to be answered urgently.

The goals of our study are to: (1) account for glacial meltwater to improve the performance of the InVEST model in arid and semi-arid areas; (2) analyze the spatiotemporal variation of water supplies and demands in Xinjiang as well as their driving factors; and (3) assess the sustainable use of water resources in Xinjiang based on ecosystem services. Specifically, the InVEST model was used to calculate WYD and social historical data were used to calculate water utilization depth (WUD) during the period of 2002–2018. The water supply rate (WSR) and water supply-demand difference (WSDD) were used to quantify the water utilization efficiency and the budget of water resources, respectively. The internal factors in changes of WYD and WUD were explored using the controlled variable method.

2 Study area and materials

2.1 Study area

Xinjiang Uygur Autonomous Region, located in Northwest China, covers an area of approximately 1.66×106km2and nearly one-six of the land area in China. Xinjiang has a typical geographic feature named as ''three mountains and two basins'', where Altay, Tianshan, and Kunlun mountains are distributed from north to south and both the Junggar Basin and Tarim Basin are situated in the gaps among the mountains (Fig. 1). As located in inland and surrounded by high mountains, the climate for most areas in Xinjiang is dry with low precipitation. The annual precipitation across the entire region is only 130.00 mm (Zhang et al., 2015). There are only a few wetlands located in the Ili River Valley and the mountainous areas. Numerous glaciers are distributed in high-altitude mountains, which contain rich freshwater that supplies most of the rivers (Chen et al., 2015; Zhang et al., 2021). In Xinjiang, although oases account for approximately 5.00% of the entire region (Cai et al., 2021), they support major population in Xinjiang. The potential evapotranspiration is high across the entire region and could reach over 1000.00 mm/a due to long sunshine duration and dry climate (Dong et al., 2020). The unique climatic and geographic conditions of Xinjiang make most rivers belonging to the inland river system. Most of the 570 rivers in Xinjiang have small runoff, and only 18 rivers of them have an annual runoff greater than 1.0×108m3(Ministry of Water Resources of the People's Republic of China, 2018).

Fig. 1 Digital elevation model (DEM) and rivers in Xinjiang (a) and land cover types and main cities in Xinjiang (b). Note that the figures are based on the standard map (新S(2021)023) of the Map Service System(https://xinjiang.tianditu.gov.cn/main/bzdt.html) marked by the Xinjiang Uygur Autonomous Region Platform for Common Geospatial Information Services, and the standard map has not been modified. Altay, Tacheng,Karamay, Shihezi, Changji, Urumqi, Turpan, Hami, Bole, Yining, Korla, Aksu, Kashgar, Artux, and Hotan are the capitals of Altay Prefecture, Tacheng Prefecture, Karamay City, Shihezi City, Changji Hui Autonomous Prefecture, Urumqi City, Turpan City, Hami Prefecture, Bortala Mongolian Autonomous Prefecture, Ili Kazak Autonomous Prefecture, Bayangol Mongolian Autonomous Prefecture, Aksu Prefecture, Kashgar Prefecture,Kizilsu Kirgiz Autonomous Prefecture, and Hotan Prefecture, respectively.

The economy and technology in Xinjiang have witnessed rapid development in recent years.The annual growth rate of gross domestic product (GDP) in Xinjiang remained above 10.00%,higher than the average in China, for most years in the past two decades (Bao and Zou, 2017).The population in Xinjiang increased from 18.5×106in 2000 to 24.9×106in 2018 (Statistical Bureau of Xinjiang Uygur Autonomous Region, 2019). Irrigation water accounted for more than 90.00% of the total water uses in Xinjiang (Wang et al., 2019). While cropping areas in Xinjiang increased from 3.38×106hm2in 2000 to 6.25×106hm2in 2018 (Statistical Bureau of Xinjiang Uygur Autonomous Region, 2019), developments and promotions of water-saving irrigation technology helped to increase the water utilization efficiency of irrigation from 0.38 in 1990(Shen et al., 2013) to 0.52 in 2015 (Wang et al., 2019).

2.2 Data acquisition

The details on all used data are shown in Table 1. To address the differences in the spatial resolution and the projection system of the data, we resampled all data in this study to a spatial resolution of 500 m under the projection of GCS_Krasovsky_1940 Albers. We prepared the model input data from 2000 to 2019 based on the official InVEST document and optimized the model parameters to ensure the model performance.

Tropical rainfall measuring mission (TRMM) data serve as the precipitation data in this study.Because TRMM satellites have high accuracy and high time resolution, they can be used in research areas with sparse meteorological stations (Luo et al., 2017; Wang et al., 2018). This study used the TRMM 3B43 v7 data product with a spatial resolution of 0.25°×0.25°. The Climatic Research Unit gridded Time Series (CRU TS) v4 dataset serves as the potential evapotranspiration data (Harris et al., 2020). The potential evapotranspiration data were calculated using the Penman-Monteith formula recommended by Food and Agricultural Organization (FAO), with a spatial resolution of 0.5°×0.5° (Dong et al., 2020; Zhou et al., 2020).The Gravity Recovery and Climate Experiment (GRACE) satellite was jointly developed by the German Space Agency (DLR) and the National Aeronautics and Space Administration (NASA).The GRACE satellite provides useful observations that reflect changes in glaciers and groundwater storage as GRACE has a low satellite orbit and is highly sensitive to the Earth's gravity field (Jin et al., 2016; Rodell et al., 2018; Beveridge et al., 2018; Meng et al., 2019). If the equivalent water column height of GRACE satellite data decreases, then the terrestrial water storage decreases correspondingly, especially on the glacier melting and groundwater extraction.

Table 1 Details of the used data

For the soil data, Chinese soil texture data including the composition of sand, clay, and silt were downloaded from the Chinese soil texture spatial distribution data of the Resources and Environmental Science and Data Center (Table 1), and soil organic carbon data (Fischer et al.,2000) in the Harmonized World Soil Database (HWSD) Chinese soil dataset (v1.1) were downloaded from the National Tibetan Plateau Data Center (Table 1). Further, we calculated the plant available water capacity (PAWC) using soil texture data mentioned above according to the formula proposed by Gupta and Larson (1979).

The European Space Agency (ESA) CCI v2.0.7 and v2.1 land cover datasets were used as the crop (irrigation) distribution data, and the distribution map of irrigated cultivated land in Xinjiang was obtained by extracting the land cover type numbered 20 with the spatial resolution of 300 m.The GDP spatial distribution data were downloaded from the Resources and Environmental Science Data Center of Chinese Academy of Sciences with the spatial resolution of 1 km (Table 1). Multiple factors such as land use type, night light brightness, and residential density that are closely related to human economic activities are considered comprehensively. Landscan data were used to be the population distribution density data with the spatial resolution of 1 km. Spatial data and multivariate asymmetric modeling methods were used to classify the census data in the administrative areas. The land use types related to ecological water were extracted from the land use data of ESA. The data of agricultural comprehensive water consumption, industrial water consumption per 10,000 CNY, domestic water per capita, and ecological water utilization volume were all from the Xinjiang Water Resources Bulletin (Water Resources Department of Xinjiang Uygur Autonomous Region, 2002–2018).

3 Methods

3.1 Water yield depth (WYD)

The InVEST model jointly developed by Stanford University, World Wide Fund for Nature, and the Nature Conservancy is widely used in the field of ecosystem services (Ochoa and Urbina-Cardon, 2017; Agudelo et al., 2020) as it contains modules that can quantify a variety of ecosystem service functions, covering terrestrial, freshwater, and marine ecosystems (Polasky et al., 2011). The InVEST model requires fewer parameters and data, in which the data availability is high. The model is usually less prone to parameter estimation errors and error propagation and is easier to interpret and understand (Tallis and Polasky, 2009). The InVEST model uses the Budyko curve in the water yield module. The Budyko curve is widely applied in studies in the UK(Redhead et al., 2016) and Australia (Donohue et al., 2012) as well as global studies (Luo et al.,2020). We briefly describe the modeling of WYD using the InVEST model as follows.

The water yield module in the InVEST model follows the water-heat coupling balance principle described by the Budyko curve. The model simplifies the confluence process and does not distinguish surface runoff and subsurface flow in the raster form (Lang et al., 2017). The formula for modeling WYD is shown in Equation 1 (Sharp et al., 2015):

whereRbdenotes the Budyko drying index; andωdenotes the ratio of plant available water depth to annual precipitation.

Both parametersRbandωin the Budyko curve can be calculated using Equations 3 and 4,respectively, as follows (Tallis et al., 2011):

where ET0denotes the potential evapotranspiration (mm);Kdenotes the evapotranspiration coefficient of plants;Zdenotes the seasonal distribution of precipitation ranging from 1 to 30; and AWC denotes the available water capacity reflecting the amount of water that plants need from soil to grow (%). AWC can be determined based on parameters such as soil properties and root depth.

The modeled results were moved to average for every four years (i.e., 2000–2003, 2004–2007,2008–2011, 2012–2015, and 2016–2019) such that we were able to investigate long-term changes in modeled results given that there was high fluctuation from year to year.

To understand the performance of the InVEST model, we used the runoff data from 14 outfall gauging stations in Xinjiang (including seven in northern Xinjiang and seven in southern Xinjiang) as the validation data for the InVEST model. The upstream catchment area of each gauging station was derived using the soil and water assessment tool (SWAT) model. We then calculated the runoff of each catchment by multiplying the simulated WYD and the corresponding catchment area and validated the modeled runoff with the measured runoff from the corresponding gauging station from 2006 to 2011. The observed values in 2007 were used as the calibrated dataset to optimize theZparameter.

Glacial meltwater is an important recharge for inland rivers in Northwest China and accounts for 9.00%–43.00% of the total runoff in inland river basins (Guo et al., 2016). It is necessary to account for water supplies from precipitation and glacial meltwater in the study of water resources in Xinjiang, while previous studies failed to consider glacial meltwater when applied the InVEST model (Wang et al., 2019; Li et al., 2020; Yang et al., 2020). In this paper, we combined glacial meltwater data from GRACE and precipitation data from TRMM as total water yield sources with the hope of improving the performance of the InVEST model according to water balance (Eq. 5). Specifically speaking, the changes in terrestrial water storage (TWSC) for the year were obtained by subtracting the GRACE data at the beginning of the year from the data at the end of the year (Eq. 6) (Yang et al., 2015; Meng et al., 2019), and then the positive values of the changes were removed. The reduction in water storage in the glacier distribution area that was extracted from the second glacier inventory data of China is considered to be the glacial meltwater (Wei et al., 2014; Guo et al., 2015).

where TWSCidenotes the changes in terrestrial water storage in theithyear (mm); and TWSi_endand TWSi_begindenote the terrestrial water storage at the end and beginning of the year (mm),respectively.

3.2 Water utilization depth (WUD)

In this study, WUD (mm) was used to characterize the potential demand for water resources and was modeled as the sum of irrigated water utilization depth (WUDirr, mm), industrial water utilization depth (WUDind, mm), domestic water utilization depth (WUDdom, mm), and ecological water utilization depth (WUDeco, mm):

To calculate WUDirr, we extracted the distribution map of irrigated land in Xinjiang from the land cover data and then assigned agricultural comprehensive water utilization depth for different cities and regions to the distribution map of irrigated croplands.

The WUDindwas calculated as follows:

whereαdenotes the proportion of industrial GDP to total GDP (%); GDP denotes the total GDP per square kilometer of land (10,000 CNY/km2); andwinddenotes the industrial water utilization per 10,000 GDP (m3/10,000 CNY).

The WUDdomwas calculated as follows:

where POP denotes the population per square kilometer of land (person/km2); andwdomdenotes the domestic water utilization per capita (m3/capita).

Ecological water utilization refers to the amount of water needed and used for ecosystem maintenance. In this study, ecological water utilization was accounted for different ecosystem types, including forests, grassland, water body (including wetland), and urban land (Li et al.,2009; Shang et al., 2013). The WUDecowas calculated as follows:

where EW denotes the ecological water utilization of each administrative region (m3); and A denotes the area of land use type (km2).

3.3 Water supply-demand relationships

Compared with the previous methods like the ecological footprint model, the perspective of ecosystem services allows for the improvement of verifiability on the assessment of water sustainability and avoids empirical parameters. For example, WYD simulated by the InVEST model is more easily verified using gauging station data than water footprint and carrying capacity of water resources. The perspective of ecosystem services can avoid empirical parameters such as global average yield capacity, yield factors, and equilibrium factors. In addition, existing studies have used ecosystem services to explore ecosystem service tradeoffs(Kang et al., 2020), human well-being (Wei et al., 2018), and ecosystem service flow (Chen et al.,2020), but few applied ecosystem services to assess the sustainable use of water resources in arid areas. The supply-demand relationship indices have been applied in many regions such as Spain(Boithias et al., 2014), Taihu Lake in China (Li et al., 2016), and coastal regions in China (Zhai et al., 2020). These indices provided a clear picture of the supplies and demands for local ecosystem services, but few applied them to assess the sustainable use of water resources. In this study, we referenced and refined these indices based on the unique circumstances of Xinjiang to assess the sustainable use of water resources more accurately. We opted to use the difference index instead of the ratio index in existing studies because the arid natural environment of Xinjiang resulted in many zero values for WYD and WUD.

WSR and WSDD were used to account for the water supply-demand relationships. WSR can reflect the utilization efficiency of water resources and was calculated as follows:

where WSR denotes the water resource supply rate; Spdenotes the potential water supply depth(mm); and Sadenotes the actual water supply depth (mm). In this study, Sprefers to the WYD modeled by the InVEST model and Sarefers to the actual water supply depth in each region as published in Xinjiang Water Resources Bulletin (Water Resources Department of Xinjiang Uygur Autonomous Region, 2002–2018). WSR indicates that water resources are not fully utilized if it is greater than 1, and otherwise it indicates that water resources are fully utilized.

WSDD that reflects the balance of water resource supplies and demands was calculated as follows:

Values of WSDD indicates that water resources are in surplus if it is much greater than zero, in balance if it is close to zero, or in deficit if it is much less than zero.

3.4 Spatial-temporal analysis and variable analysis

Cold and hot spot analysis is used for understanding the spatial patterns and distributions of water supplies and demands in our study. The Gi*statistical method is able to reflect the spatial aggregation of points or patches with high (or low) values (Ord and Getis, 1992):

Gi*can be statistically tested based on theZscore and theP-value.Zscore is the standard deviation of Gi*to measure the aggregation degree.P-value is the significance level to determine whether spatial aggregation is produced. IfZscore is positive, it is a hot spot area that clusters with high values, and ifZscore is negative, it is a cold spot area that clusters with low values. TheP-value is indicative to the significant level of the cold or hot spots (Guerri et al., 2021). The cold and hot spot analysis was used to analyze the aggregations of high (or low) values about WYD,WUD, WSR, and WSDD.

Both the Sen trend analysis and Mann-Kendall (MK) test were used for time series analysis.The Sen trend analysis is a robust non-parametric statistical method widely used to derive the trend of time series data. It is insensitive to outlier data and errors but does not reflect the significance level of the trend. The Sen trend analysis is derived as follows:

The MK test is a non-parametric statistical test method that is able to conduct a significance test for the trend of the time series and is not sensitive to outliers:

whereZmis the standardized test statistic, obeying the standard normal distribution;Smis the test statistic; sign is the sign function; andlis the number of data. If |Zm|>1.96, the data series has obvious trends at the significance level of 0.05, otherwise, the data series has insignificant trends at the significance level of 0.01 (Schr?ter and Remme, 2016). The Sen trend analysis and MK test were used to analyze temporal trends about WYD, WUD, WSR, and WSDD.

Different variables that are input to the model play varied roles in the model. The controlled variable method was used to investigate the contribution of each input variable to the modeled result. Taking the water yield module of the InVEST model as an example, variables such as precipitation, potential evapotranspiration, and land cover data are key inputs to the model for predicting WYD. To obtain the contribution rate of precipitation to WYD, we treated the precipitation data in 2002 as the baseline and used them to simulate land surface processes in 2006, 2010, 2014, and 2018, while keeping the other variables the same. The contribution rate of the variable to the output results can be calculated as follows:

wherePijdenotes the contribution rate of variablejto the result in theithyear;Vactdenotes the result obtained by inputting all variables in theithyear according to the actual situation; andVijdenotes the result obtained by inputting the variablejin the baseline year and other variables in theithyear. To further understand the internal factors of changes in supplies and demands of water resources, we investigated the relative contributions of precipitation, potential evapotranspiration,and land cover to the modeled WYD, as well as the relative contributions of WUDirr, WUDind,WUDdom, and WUDecoto the total WUD.

4 Results

4.1 Model calibration and validation

The calibration experiment showed that the root mean square error (RMSE), mean absolute error(MAE), and Pearson's correlation coefficient were 5.22×108(±1.29×108) m3, 4.08×108m3, and 0.65, respectively. The validation experiment using remaining four-year data shows that the RMSE, MAE, and Pearson's correlation coefficient were 9.67×108(±2.82×108) m3, 6.87×108m3,and 0.72, respectively, when comparing modeled and observed runoff values.

Mean relative bias error between observed runoff data and modeled results without considering glacial meltwater was 1.96×108m3. By comparison, mean relative bias error between observed runoff data and modeled results considering glacial meltwater using the GRACE data decreased to 1.29×108m3(Fig. 2).

4.2 WYD and WUD

WYD in Xinjiang had a spatial distribution pattern close to precipitation and generally decreased from the north to the south and from the west to the east with a range from 0.00 to 800.00 mm(Fig. 3a1–a5). The areas with WYD above 50.00 mm were mainly distributed in northwestern Xinjiang and sparsely distributed at high altitudes in both the Tianshan Mountains and Kunlun Mountains. These areas, accounting for 5.26% of total area in Xinjiang, were often rich in water supplies and were shown as the hot spots above the 95.00% confidence level (Fig. 3b1–b5). The areas with WYD ranging from 0.10 to 50.00 mm were mainly distributed in both northern and southern regions of the Junggar Basin and in the oasis areas surrounding the Tarim Basin.Although there was runoff in these basin areas, WYD was relatively small as compared to those in wet areas. The areas with WYD ranging from 0.10 to 50.00 mm were mostly identified as the cold spots above the 95.00% confidence level, and about 0.13% of total area in Xinjiang belonged to these areas. The areas with WYD below 0.10 mm were mainly distributed in deserts which almost have no runoff, including the Taklimakan Desert, the Gurbantungut Desert, and the Kumtag Desert. Their total area accounted for more than 89.00% of the land area in Xinjiang,where deserts were widely distributed. Therefore, these areas were considered as insignificant areas in the cold and hot spot analysis as a result of non-aggregated effects.

Fig. 2 Mean relative bias error between observed runoff data and modeled results without considering glacial meltwater (red violin boxes) and with considering glacial meltwater (blue violin boxes). The hollow white dots represent the median. The upper and lower horizontal lines represent the maximum and minimum values,respectively. The bottom and top lines of the box represent the 25.00% and 75.00% quartiles, respectively. The width of the violin represents the frequency of the data.

In terms of temporal changes, WYD had a continuous upward trend (Fig. 4a), with increases from 8.51 mm in 2002 to 10.86 mm in 2018, rising by 27.61%. The increasing trend was significant (P<0.05) and the Pearson's correlation coefficient across time was 0.91. The Sen trend analysis showed that 70.88% of the study area that passed the trend analysis exhibited increasing trends of WYD. The Sen slope was high in some areas of northwestern Xinjiang. The MK test was significant in the north and south sides of the Junggar Basin (Fig. 4b), indicating that WYD was significantly increasing in these areas. There were some areas with decreased WYD, mainly concentrated in the east of both the Tianshan Mountains and Kunlun Mountains and Ili River Valley.

WUD was found high in oasis areas and cities and low in desert areas (Fig. 5a1–a5). WUD in Xinjiang mainly ranged from 0.00 to 4000.00 mm. The areas with WUD above 500.00 mm were mainly distributed in urban and irrigated croplands such as the agricultural areas near the Irtysh River, Ili River, and Tarim River, and in the oases near the northern slope of the Tianshan Mountains. These areas were often identified as the hot spots above the 95.00% confidence level,accounting for 6.35% of the total area in Xinjiang (Fig. 5b1–b5). The areas with WUD below 1.00 mm were widely distributed in the desert areas and not clustered. Many areas with WUD ranging from 1.00 to 100.00 mm were identified as the cold spots above the 95.00% confidence level and were mainly distributed in the mountainous areas of the Altay Mountains, Tianshan Mountains,and Kunlun Mountains. These areas totally accounted for 13.94% of the total area in Xinjiang.

My sledge! Don t forget my sledge! He thought of that first, and it was fastened to one of the great white birds who flew behind with the sledge on its back

In terms of temporal changes, averaged WUD had a downward trend and its correlation with time had a Pearson's correlation coefficient of –0.85. Averaged WUD descended from 67.50 to 53.53 mm in the study period, reducing by 20.70% (Fig. 6a). The Sen trend analysis showed that WUD in the hot spots mostly had a downward trend (Fig. 6b) and the MK test showed that the the north of the Urumqi City (provincial capital of Xinjiang) and the central region in the Changji Hui Autonomous Prefecture experienced significantly decreased WUD above the 95.00% confidence level. The areas that experienced decline in WUD accounted for 43.21% of the total area in Xinjiang, and the areas with significant decline in WUD were 1.80 times the areas with significant increase. The areas with increased WUD were mainly located in the Altay Mountains,Tianshan Mountains, and Kunlun Mountains. The MK test showed that the increasing trend of WUD was significant in these mountainous regions. WUD was also increased in the downtown area of Urumqi City.

4.3 WSR and WSDD

In Xinjiang, 93.39% of the areas had WSR less than 1.00; WSR in a few areas was greater than 1.00. Most of these areas were identified as the hot spots (Fig. 7), such as areas in the Altay Prefecture and northern Tacheng Prefecture, around the Tianshan Mountains, and in the west of the Kunlun Mountains. These areas were highly overlapped with the hot spots of WYD. Areas with WSR less than 1.00 were widely distributed in the study area and no obvious spatial aggregation could be found in these areas, so the WSR had little cold spots.

Fig. 3 Spatial distribution of water yield depth (WYD) in Xinjiang from 2002 to 2018 (a1–a5) and the corresponding cold and hot spot analysis results from 2002 to 2018 (b1–b5). Note that the figure is based on the standard map (新S(2021)023) of the Map Service System (https://xinjiang.tianditu.gov.cn/main/bzdt.html)marked by the Xinjiang Uygur Autonomous Region Platform for Common Geospatial Information Services, and the standard map has not been modified.

Fig. 4 Temporal analysis on averaged WYD in Xinjiang from 2002 to 2018 (a) and the Sen trend analysis and Mann-Kendall (MK) test of WYD in Xinjiang from 2002 to 2018. Note that the figure is based on the standard map (新S(2021)023) of the Map Service System (https://xinjiang. tianditu.gov.cn/ main/bzdt.html) marked by the Xinjiang Uygur Autonomous Region Platform for Common Geospatial Information Services, and the standard map has not been modified.

Fig. 5 Spatial distribution of water utilization depth (WUD) in Xinjiang from 2002 to 2018 (a1–a5) and the corresponding cold and hot spot analysis results from 2002 to 2018 (b1–b5). Note that the figure is based on the standard map (新S(2021)023) of the Map Service System (https://xinjiang.tianditu.gov.cn/main/bzdt.html)marked by the Xinjiang Uygur Autonomous Region Platform for Common Geospatial Information Services, and the standard map has not been modified.

Fig. 6 Temporal analysis on averaged WUD in Xinjiang from 2002 to 2018 (a) and the Sen trend analysis and MK test of WUD in Xinjiang from 2002 to 2018. Note that the figure is based on the standard map (新S(2021)023) of the Map Service System (https://xinjiang.tianditu. gov.cn/main/ bzdt.html) marked by the Xinjiang Uygur Autonomous Region Platform for Common Geospatial Information Services, and the base map has not been modified.

During the study period, averaged WSR in Xinjiang decreased from 0.30 to 0.22, descending by 26.67% (Fig. 8a) and its correlation with time reached a Pearson's correlation coefficient of–0.79. The areas with increased WSR were mainly distributed in the oases surrounding the Junggar Basin and Tarim Basin (Fig. 8b). The areas with decreased WSR included areas in the east of the Tianshan Mountains and in the northern Bayangol Mongolian Autonomous Prefecture,Ili Kazak Autonomous Prefecture, and northern Altay Prefecture.

In the study area, WSDD mostly ranged from –4000.00 to 800.00 mm (Fig. 9a1–a5) and averages of WSDD generally ranged from –59.00 to –42.00 mm. The areas with negative values of WSDD accounted for 82.40% of the total area in Xinjiang. Spatially, WSDD was often found in surplus in mountainous areas and in deficit in irrigation areas. Large areas with WSDD greater than 30.00 mm were identified as the hot spots and were located in high-altitude areas. The other hot spots were areas represented by Altay Prefecture, Tacheng Prefecture, and Bortala Mongolian Autonomous Prefecture (Fig. 9b1–b5). The areas with hot spots above 90.00% confidence level accounted for 27.81% of the total area in Xinjiang. The areas with WSDD less than –500.00 mm were often identified as cold spots, of which the spatial distribution was highly coincident with the hot spots of WUD. The areas with cold spots above 90.00% confidence level accounted for 7.73% of the total area in Xinjiang. The other areas where WSDD was close to 0.00 mm were considered as non-significant areas in the cold and hot spot analysis.

Averaged WSDD in Xinjiang increased from –59.10 mm in 2002 to –42.72 mm in 2018(P<0.05; Fig. 10a). The Sen trend analysis showed that WSDD in 64.65% of the study area that passed the trend analysis increased, especially in the northern Tacheng Prefecture, Bortala Mongolian Autonomous Prefecture, and in the west of the Tianshan Mountains and Kunlun Mountains. The MK test showed that WSDD in the rim of the Junggar Basin significantly increased with high confidence level (Fig. 10b). There were also some areas where WSDD decreased. For example, tensions in water supply-demand relationships occurred in the east of the Tianshan Mountains and Kunlun Mountains. Areas in the southern slope of the Tianshan Mountains and in the Ili River Valley also presented tensions in water supply-demand relationships.

4.4 Controlled variable analysis

Fig. 7 Spatial distribution of water supply rate (WSR) in Xinjiang from 2002 to 2018 (a1–a5) and the corresponding cold and hot spot analysis results from 2002 to 2018 (b1–b5). Note that the figure is based on the standard map (新S(2021)023) of the Map Service System (https://xinjiang.tianditu.gov.cn/main/bzdt.html)marked by the Xinjiang Uygur Autonomous Region Platform for Common Geospatial Information Services, and the standard map has not been modified.

Fig. 8 Temporal analysis on averaged WSR in Xinjiang from 2002 to 2018 (a) and the Sen trend analysis and MK test of WSR in Xinjiang from 2002 to 2018 (b). Note that the figure is based on the standard map (新S(2021)023) of the Map Service System (https://xinjiang. tianditu.gov.cn/main/ bzdt.html) marked by the Xinjiang Uygur Autonomous Region Platform for Common Geospatial Information Services, and the standard map has not been modified.

Fig. 9 Spatial distribution of water supply-demand difference (WSDD) in Xinjiang from 2002 to 2018 (a1–a5)and the corresponding cold and hot spot analysis results (b1–b5). Note that the figure is based on the standard map (新S(2021)023) of the Map Service System (https://xinjiang.tianditu.gov.cn/main/bzdt.html) marked by the Xinjiang Uygur Autonomous Region Platform for Common Geospatial Information Services, and the standard map has not been modified.

The controlled variable method was used to analyze the relative contribution of precipitation,potential evapotranspiration, and land cover to WYD, as well as the contributions of WUDirr,WUDind, WUDdom, and WUDecoto the total WUD (Fig. 11). The contribution rates of variables could reflect the importance of factors on the changes in WYD. Figure 11a showed that precipitation had large impacts on WYD and the contribution rate of precipitation to WYD reached 19.52% in 2018. The contribution rate of precipitation to WYD increased with increasing precipitation in Xinjiang over time (Fig. 11a and c). The contribution of potential evapotranspiration and land cover to WYD were both maintained at less than 6.00%. Potential evapotranspiration had negative impacts on WYD, and WYD decreased with increasing potential evapotranspiration (Fig.11a and d). Land cover changes in Xinjiang showed a positive contribution to WYD, and the contribution rate had a slight upward trend. In 2018, the contribution rate of land cover changes to WYD reached 5.52%. The land cover transfer matrix shown in Table 2 indicated that large areas of bare land were transformed into grassland, ice and snow, cropland, shrubland, and water body. The area of bare land transformed into grassland reached 40,470 km2, accounting for nearly 3.40% of the total bare land in 2002. Land cover changes in Xinjiang led to increases in water yields.

Fig. 10 Temporal analysis on averaged WSDD in Xinjiang from 2002 to 2018 (a) and the Sen trend analysis and MK analysis of WSDD in Xinjiang from 2002 to 2018 (b). Note that the figure is based on the standard map (新S(2021)023) of the Map Service System (https://xinjiang. tianditu.gov.cn/ main/bzdt.html) marked by the Xinjiang Uygur Autonomous Region Platform for Common Geospatial Information Services, and the standard map has not been modified.

Fig. 11 Contributions of precipitation, potential evapotranspiration (PET), and land cover to WYD (a),contribution of irrigated WUD (WUDirr), industrial WUD (WUDind), domestic WUD (WUDdom), and ecological WUD (WUDeco) to the total WUD (b), and changes in precipitation (c) and PET (d) in Xinjiang from 2002 to 2018

The results shown in Figure 11b indicated that the contribution of WUDirraccounted for a majority of nearly –30.00% in changes of the total WUD. The relative contribution of WUDirrwas negative, different from WUDindand WUDdom. The reduction in WUDirrhas become increasingly important to the changes in the total WUD over time. The relative contributions of WUDindand WUDdomto the changes in the total WUD were relatively low but positive, and both of them experienced an increasing trend. The contribution rate of WUDecowas maintained at a relatively unimportant level across years.

Table 2 Land cover transfer matrix in Xinjiang from 2002 to 2018 (unit: km2)

5 Discussion

The mean relative bias error of the InVEST model decreased by 34.00% when we accounted for glacial meltwater using the GRACE data. Our study found that considering glacial meltwater into the InVEST model using satellite-based water storage data could help to decrease bias errors in arid and semi-arid areas such as Xinjiang. Previous studies indicated that incorporating glacial meltwater into water supply could improve the performance of the InVEST model in areas with high altitudes and in areas with glaciers (Scordo et al., 2018; Benra et al., 2021). Although the Budyko curve did not account for detailed hydrological processes, the simple hydrological model was accurate enough in areas without precise input data. According to validation results, the accuracy of the InVEST model could meet the requirement of studies on ecosystem services at the annual scale (Pan et al., 2015; Scordo et al., 2018).

High WYD was mainly distributed in two major areas: northwestern Xinjiang and mountainous areas. Precipitation was the main source of WYD in northwestern Xinjiang, while glacial meltwater was the main source of WYD in mountainous areas. It's worth noting that decreased WYD mainly occurred in the east of both the Tianshan Mountains and Kunlun Mountains with fewer glaciers, possibly due to the arrival of the "bonus period". Global warming resulted in considerable glacier melting in the Tianshan Mountains and Kunlun Mountains, leading to increased runoff. However, it was an unsustainable way to increase water resources. After the"bonus period" of glacial melting, runoff will decrease significantly, and water resources will correspondingly decrease significantly. This is also a major problem for water resources in the future (Hartmann et al., 2016; Bolch et al., 2021).

Decreased WUD principally occurred in the oasis that contained vast irrigated croplands.Agricultural irrigation was the largest consumption of water resources in Xinjiang (Chen et al.,2016). As technology advances, water-saving irrigation technology is becoming more and more improved and is widely used in agriculture. By the end of 2016, water-saving agricultural area accounted for 74.59% of the total cultivated area in Xinjiang, which ranked the first in China(Cao et al., 2020). Water-saving agriculture has contributed greatly to the reduction of WUD.Increased WUD was mainly distributed in mountainous areas. Since 2000, ecological engineering projects have greatly improved vegetation cover in Xinjiang (Niu et al., 2019), and warmer and wetter climate also resulted in increased Normalized Difference Vegetation Index (NDVI) in Xinjiang (Xu et al., 2016; Luo et al., 2020). In a word, vegetation status was slightly improved,especially in mountains where vegetation was concentrated, which can increase the consumption of ecological water. It's worth noting that the hot spots of WYD and WUD mismatched spatially according to the spatial analysis. Mountainous areas were generally hot spots of water yield but cold spots of water utilization, whereas oases were hot spots of water utilization but cold spots of water yield.

Water resources in most of Xinjiang with WSR less than 1.00 were fully utilized and local water supplies needed new facilities. Water yields in northwestern Xinjiang were not fully used for water supply. The hot spots and increased areas of WSR were in step with WYD, so the spatial and temporal patterns of WSR were greatly affected by WYD. During the study period, water resources in 82.40% of Xinjiang were in deficit, and previous researches had similar conclusions.Specifically, current water resources cannot meet leap-forward development in Xinjiang (Dai et al., 2017), and water resources in Alar, a city located in southern Xinjiang, also cannot satisfy demand for economic development (Liu et al., 2019). Fortunately, WSDD in 64.65% of the region that passed the trend analysis exhibited an increasing trend, and the areas with increased WSDD were mainly distributed in northern Xinjiang, indicating that water deficit in some areas of Xinjiang has been alleviated. Exiting researches showed that the improvements of water-saving technology and ecological management has alleviated the water shortage in Xinjiang (Ling et al.,2019; Zhou et al., 2020).

According to the controlled variable method, precipitation dominated the changes in WYD, and the contribution rates of precipitation to WYD increased as precipitation rose. Previous studies also indicated that water yield prediction was more sensitive to precipitation than to potential evapotranspiration in Ethiopia (Sahle et al., 2019), and the similar conclusion was confirmed in UK (Redhead et al., 2016). The Budyko curve used in the InVEST model allocated precipitation into evapotranspiration and runoff according to the aridity index. Increase in potential evapotranspiration often led to decrease in runoff if precipitation remains unchanging. While local climate and environmental changes have led to increase in WYD in Xinjiang in the past years,people still need to pay attention to uncertainties in water resources due to global climate change.When the climate became warmer and wetter, increases in precipitation positively impacted on water yield, but increases in potential evapotranspiration had negative effects on water yield.Further research and continuous observations may be required to determine which of these two factors is stronger. As for WUD, controlled variable analysis illustrated that WUDirrdominated the changes in WUD, and the contribution rates of WUDirrto WUD decreased as the WUDirrdecreased. Previous studies also indicated that irrigation accounted for the largest proportion of total water consumption in Xinjiang, up to 97.00% (Chen et al., 2016), and the similar conclusion also appeared in the study of Guo et al (2015). In addition to the improvement of water-saving agriculture, increases in precipitation also mitigated water use pressures on agricultural irrigation.The contribution rates of WUDindand WUDdomwere relatively low but positive with an upward trend. Local population, economy, and urbanization rate were all increased in recent years and positively contributed to WUD (Tang et al., 2012; Gong et al., 2020). It was worth noting that upgrades in industrial structure and increases in the proportions of secondary and tertiary industries in Xinjiang led to higher WUDindand WUDdom. It is necessary for Xinjiang to control the total WUD through measures such as developing water-saving industries and advocating life-saving water conservation strategies. As WUDirrdominated changes in the total WUD in Xinjiang in the past years, it is important to continue prompt water-saving irrigation technologies and also control the expansion of cropland areas. Reasonable planning of the cropland layout and expansion is also needed (Lei et al., 2012; Dai et al., 2017).

This research focused on discussing the inter-annual spatial distributions and temporal changes of water resource supplies, demands, and their relationships. Due to climatic and topographical conditions in Xinjiang, intra-annual variation in water resources is also considerable. On the water supply side, northern Xinjiang, western Xinjiang, and high-altitude areas often experienced large amounts of snowfall in winter. Although snowfall in one year tended to accumulate in the form of solid water, it did not generate runoff in the corresponding year but gradually turned into runoff with increasing temperature in the following year. The transformation of snow cover to runoff had a lagging effect, meaning that most snowfall benefited water supply services in the following year. There is a need to refine the temporal scale to the monthly scale instead of the annual scale in future studies to account for the lagging effect. On the water demand side, the peak of WUDirrgenerally occurred from July to August (Shen et al., 2013) and coincided the period when temperature and evapotranspiration were high, resulting in seasonal water shortage of agricultural irrigation and decrease of water utilization efficiency. There is a need to study ways of storing water resources in surplus months for use in deficit months.

This research assessed the sustainable use of water resources from the perspective of ecosystem services. As ecosystem services were broad topics, the trade-off and synergy among service functions could be further discussed to better reflect the integrity of the ecosystems (Hong et al.,2020). With a better understanding of ecosystem services, research hot spots continued to emerge,including ecosystem service evaluations, ecosystem service flows, and ecological compensation(Yu and Bi, 2011). The method to study ecosystem services could extend from one single function to multiple functions, which could help deepen and broaden the sustainable assessment of water resources under the framework of ecosystem services.

6 Conclusions

This study evaluated the supply-demand relationships of water resources in Xinjiang during the period of 2002–2018 from the perspective of ecosystem services. The InVEST model was used to simulate WYD and socioeconomic data were used to model WUD. WSR and WSDD served as indices to account for the supply-demand relationships of water resources, and then the controlled variable method was used to analyze the internal factors of changes in the supplies and demands of water resources. The results showed that WYD experienced a significant increasing trend,especially in oasis areas. WUD had a downward trend, reducing from 67.50 to 53.53 mm. The decreased WUD were also mainly distributed in oasis areas. During the study period, water resources in Xinjiang were in a slight deficit but mitigated over time with climate change and human activities. Increased precipitation dominated increased WYD, and decreased WUDirrdominated decreased WUD. Water resources in nearly 93.39% of Xinjiang were fully utilized and the utilization rate of water resources increased during the study period. The use of satellite-based water storage data helped to decrease the bias error of the InVEST model. The contradiction of water resources in Xinjiang is easing, but it is still necessary to be on the alert for the worsening contradiction between supply and demand in the later stage of glacier melting.

From the perspective of ecosystem services, this study assessed the water supply-demand relationships and evaluated water resource sustainability in arid and semi-arid areas, which diversified the evaluation methods of water resource sustainability by addressing some issues in previous methods and also extended the application scenarios of researches on the supply-demand relationships of ecosystem services. Compared with the previous methods like the ecological footprint model, considering ecosystem services in the InVEST model can improve the verifiability on assessment of water sustainability and avoid empirical parameters. This research may provide a reference for evaluation on sustainable uses of water resources in the other regions in the world.

Acknowledgements

This research was supported by the National Natural Science Foundation of China (41875122), the Western Talents (2018XBYJRC004), the Guangdong Top Young Talents (2017TQ04Z359), the Introducing Talents to Western China Project of Chinese Academy of Sciences (Y932121), and the Natural Science Foundation of Guangdong Province, China (2021A1515011429). We thank the Water Resources Department of Xinjiang Uygur Autonomous Region, China for providing the Xinjiang Water Resources Bulletin for this study. We also acknowledge the Resource and Environmental Science Data Center of Chinese Academy of Sciences for providing the spatial distribution of GDP. Special thanks are also given to the anonymous reviewers for their constructive comments.