Correlation assessment of NDVI and land use dynamics with water resources for the southern margin of Mu Us Sandy Land, China

To prevent desertification, countries worldwide have made diversified efforts, and vegetation restoration has been demonstrated to be an effective approach. However, in regard to sandy land with limited water resources, measures such as revegetation may lead to an increased drought risk. Despite confirmed sand utilization achievements, many controversies remain regarding the advantages of desert greening, especially considering water scarcity. Therefore, the long-run and causal relationships between sandy land, water consumption, and vegetation coverage are necessarily explored. Choosing the southern margin of the Mu Us Sandy Land as the study area, this study explored the interactions between sandy land, water consumption, and normalized difference vegetation index (NDVI) of 2000–2018 with the vector autoregression (VAR) model approach. In the study area, various revegetation projects have been implemented, resulting in a notable reduction in the sandy land area. In addition, the NDVI increased from 0.196 in 2000 to 0.371 in 2018, an increase of 89.3%. The results indicated that there exist long-term stable equilibrium and causal relationships existed between water consumption and sandy land and NDVI. NDVI enhancement is relatively the direct factor that causes the increase of water consumption. It could be inferred that the implemented revegetation measures may rely on a large water consumption amount, which may further aggravate water shortages and ecological damage issues. More scientific and stronger effective water resource management measures should be locally implemented to achieve a balance between water resources and revegetation.


Introduction
Desertification is one of the most serious ecological problems globally. Drylands cover approximately 41% of the Earth land surface and are home to more than 38% of the total global population of 6.5 billion people (Reynolds et al. 2007). To address the issue of desertification, governments worldwide have implemented various measures, such as sand fixing, revegetation (Amiraslani and Dragovich 2011), windbreaks installation, and natural sealing (Zhu et al. 1984;Bonkoungou 1996). China, which is one of the countries most seriously affected by desertification, has implemented large-scale conservation programs, including the Key Shelterbelt Construction Program, Wildlife Conservation and Nature Reserve Development Program, Forest Eco-Compensation Program, and Natural Forest Conservation Program (Liu et al. 2008;Yin and Yin 2010). The implementation of these projects has achieved remarkable results. According to the NASA satellite data, the global green area increased by 5% from 2000 to 2017 (equivalent to the entire Amazon Responsible Editor: Philippe Garrigues * Siyuan Liu liusiyuan22@126.com * Hongzhen Liu liuhzh@yrec.cn rainforest). Notably, China has led the growth of global greening and contributes 25% of the net growth in the global vegetation area (Chen et al. 2019). The Mu Us Sandy Land, which is located in the typical agro-pastoral transitional zone, is one of the most sensitive, vulnerable, and severely degraded areas in northern China. This desert area contains abundant water resources and a high grass coverage, but with increasing population and human activities, excessive reclamation and grazing have led to water resource depletion and continuous desertification (Runnstrom 2003;Huang et al. 2009;Bai and Cui 2019). However, in recent years, the trend of desertification has been reversed. The desertification land area of the Mu Us Sandy Land decreased at an average rate of 62.37 km 2 /year from 1990 to 2017 (the desertification area in 2017 amounted to 1684.09 km 2 ) (Han et al. 2019). Over the past three decades, the interannual normalized difference vegetation index (NDVI) in the Mu Us Desert has mostly exhibited increasing trends due to human activity impacts (Karnieli et al. 2014;Wu et al. 2014;An et al. 2014;Li et al. 2016). The area with the most dynamic NDVI trends is distributed in Shaanxi Province (Wu et al. 2002;Liu et al. 2020). This area occupies one-third of the area of the Mu Us Sandy Land and is referred to as the agropastoral ecotone of the Northern Shaanxi Province (ANS). According to data released by the Shaanxi Provincial Forestry Bureau on April 22, 2020 (Yulin Municipal People's Government 2020), the desertification land control rate in Yulin has reached 93.24%, which indicates that all mobile sand dunes have been fixed, thus suggesting that a real realization of "sand degradation and green coverage increase".
Due to the distinctiveness of the study location and the desertification prevention achievements in recent years, the ANS has become a research area of heightened interest for scholars. Numerous studies have explored the effectiveness of the aforementioned projects. Studies have assessed desertification causes (Wu and Ci 1999;Wu 2001), monitored and evaluated desertification (Wu and Ci 2002;Zhang et al. 2003), examined land use and land cover (LULC) changes (Zhang et al. 2012;Meng et al. 2012;Li et al. 2013), explored vegetation change trends (Ju et al. 2008;Karnieli et al. 2014), and conducted ecological risk and ecological vulnerability assessments in desert areas (Meng et al. 2015). However, the transformation process of a natural desert ecosystem is a complex coupling process, and the resultant impacts should be further investigated.
In previous studies, both field investigations and statistical data analysis have revealed that, owing to the common influence of growing human activities and still fragile eco-environments, ecological restoration in the Mu Us Sand Land is experiencing increasing challenges due to water scarcity (Li et al. 2017). A previous study indicated that although revegetation yielded positive effects on ecological systems in the Mu Us Sandy Land, revegetation caused an evapotranspiration increase of 51 mm per summer, thus suggesting that revegetation exacerbated water shortage problems in the aforementioned region (Zheng et al. 2020). In water-limited arid areas, a reduction in water resources due to vegetation commonly results in a reduction in the amount of water available for human activities. Consequently, the competition for water between the environment and humans is exacerbated, which results in ecological degradation (Quiggin 2001;Wu et al. 2002). Field investigations and statistical data analysis have demonstrated that water-scarcity-related challenges are associated with ecological restoration efforts in the Mu Us Desert due to growing human activities and fragile regional eco-environment (Li et al. 2016(Li et al. , 2017. In terms of sandy land containing limited water resources and highly susceptible to drought, measures, such as vegetation planting, may lead to an increase in water consumption and new desertification crisis. Despite the reported sand utilization achievements, many controversies remain regarding the advantages of desert greening measures, especially considering the scarcity of water resources in these areas. Therefore, whether the sandy land degradation and vegetation growth caused by human activity may lead to a substantial increase in water resource consumption must be considered. To address this issue, this study adopts the ANS as the study area and aims (1) to analyze the change in land use in the study area and clarify the conversion relationship between sandy land and other land types, (2) to explore NDVI changes across the ANS, and (3) to assess the dynamic relationships among sandy land, water resources, and NDVI based on the vector autoregression (VAR) model. This study is expected to provide references for ecological improvement and sustainable development in water-deficient areas.

Description of the study area
In this study, ANS is taken as the study area, which occupies one-third of the area of the Mu Us Sandy Land. The ANS is located in Northwest China (ranging from 107°35′E to 111°29′E and from 37°35′N to 39°02′N) and covers an approximate area of 33,992 km 2 (Figure 1). This area comprises six counties, namely Yuyang, Shenmu, Fugu, Hengshan, Dingbian, and Jingbian. More than 44.2% of the study area is classified as desertified land. The annual mean temperature is approximately 6-8.5°C, with monthly mean temperatures of 22°C in July and −11°C in January. The considered region experiences a typical continental semiarid climate, and the water resources here are very limited. The per capita water resource level reaches only 9.97 million m 3 , which is lower than the water resource standard proposed by the United Nations Population Action Organization in 1993 for severely water-scarce countries (10 million m 3 /per capita) (Wu and Ci 2002). The average annual rainfall of ANS is 401 mm. And the spatial distribution of rainfall is very uneven, ranging from 440 mm in the southeast to 250 mm in the northwest. A total proportion ranging from 60 to 80% of the annual precipitation is received from June to August. Based on the historical rainfall in the study area, the precipitation in this area showed a slight growth trend in 2000-2018, with an average annual growth of about 4.3 mm/a and a total growth rate of about 19.38% in 19 years. However, the increase rate of water consumption in this area is about 44.68% correspondingly. For the study area which suffers water resources depletion, the rapid growth of water consumption is unfavorable for the sustainable development of the region.

Data sources
The land use type reflects the achievements of humans in regard to land use and transformation as well as the various forms and uses (functions) of land. In this study, land use data from 1990 and 2000-2018 were adopted. The NDVI is a quantitative parameter of the vegetation coverage and reflects the ecological environmental quality. Monthly average NDVI data were obtained from the spatial distribution dataset of the monthly NDVI in China (2000China ( -2018. Land use and NDVI data could be accessed from the Resource and Environment Data Cloud Platform, Chinese Academy of Sciences (http:// www.resdc.cn/), and the NDVI dataset named "China monthly vegetation index (NDVI) spatial distribution dataset" (https://www.resdc.cn/DOI/doi.aspx?DOIid=50 ). Water resource data could be accessed from the Water Resources Bulletin of Shaanxi Province and Water Resources Statistics Bulletin of Yulin city, which are available on the websites of relevant government departments (http://slt.shaanxi.gov.cn/ zfxxgk/fdzdgknr/zdgz/szygb/ and http://slj.yl.gov.cn/index. html, respectively).

Analysis of the land use and NDVI changes
To better understand the land cover and NDVI change trends in the ANS, the evolution characteristics of land cover and NDVI changes were carried out. With the ArcGIS platform, land cover maps and NDVI change maps pertaining to 2000-2018 were overlaid, and key information on areas exhibiting major NDVI changes and the land cover type across this area was extracted.
Additionally, a transition matrix, a classic method for LUCC detection, was employed to explore the dynamics of LUCC in the ANS. In regard to land cover changes from 2000 to 2018, the typical year (2000, 2005, 2010, 2015, and 2018) is selected as the time node in order to better show the conversion of land use types. As the large-scale expansion of Fig. 1 Location of the study area sandy land use scale in the study area mostly started since 2000, to achieve a better comparison, data of 1990 was chosen for further comparison. In regard to each pair of compared data sets, an extended transition matrix was constructed. Maps reflecting the initial and subsequent times were overlaid in ArcGIS to produce a matrix providing the LULC area according to the categorical transition between two time points. Due to the large amount of data, to better reveal the flow direction of the considered data, a Sankey chart was selected to visualize the transfer matrix data.

A VAR model approach
The VAR model is a widely adopted econometric technique for multivariate time-seriesmodeling and was proposed by Sims (1980) to analyze the mutual influence relationship among macroeconomic variables. The VAR model exhibits very attractive features and has proven to be a valuable tool for the analysis of dynamic transmission mechanisms among time-series processes (Mcmillin 1991;Lu 2001). In each model equation, the endogenous variables return all lagged endogenous variable items to estimate the dynamic relationships occurring among all the endogenous variables. Although the VAR model was originally developed to analyze the mutual influence relationship among macroeconomic variables, the model relies on multiple simultaneous equations, which are not based on economic theory and can be applied to several other fields. For example, the VAR model was applied to analyze and forecast the mutual relationship between variables in a system, determine leading and potential factors, and quantify the influence of these factors (Kumar et al. 2009;Adenomon et al. 2013;Wu et al. 2018).
In addition to forecasting, VAR model equations have been adopted to simulate the effects of sudden changes (impulses) in one variable on other variables. These effects are quantified through the impulse response function (IRF), which enables estimation of the timescale over which changes in water consumption, sandy land area, and NDVI affect each other. Another application of the VAR model is variance decomposition considering the forecast error variance. Variance decomposition indicates the contribution (as a percentage) of each model variable to the forecast error variance. The implementation of the VAR model is described in the following text.

Augmented Dickey-Fuller test
Standard Granger causality tests were conducted on the above stationary time series. The unit roots of X t were determined to confirm the stationary properties of each variable. The aforementioned unit roots were obtained with the augmented Dickey-Fuller(ADF) test Fuller 1979, 1981). Regarding time series X t , the ADF relationship is expressed as follows: where △ is the difference operator; k is the autoregressive lag length, which should be sufficiently large to eliminate any possible serial correlation in β i , and α and β are the coefficients of interest. When the aforementioned variables were found to be nonstationary, ADF tests were repeated targeting the first and second differences.

Johansen multivariate cointegration test
The Johansen multivariate cointegration test (Johansen and Juselius 1990) can be represented as follows: where Z t is a 3 × 1 vector of the variables lnU t , lnY t , and lnE t , μ is a 3 × 1 vector of constant terms; Γ and Π denote 3 × 3 matrices of the coefficients; and ε t is a 3 × 1 vector of white noise error terms. The Johansen multivariate cointegration test is based on maximum likelihood estimation and trace statistics (λ trace ), where the λ trace statistic tests the null hypothesis of no cointegration (rank (r) = 0) against an alternative hypothesis of cointegration (r > 0).

Selection of the lag length
Most VAR models are estimated considering symmetric lags, i.e., the same lag length is considered for all variables in all model equations. The lag length is mostly selected based on criteria such as the final prediction error (FPE) and Akaike information criterion (AIC). These information criteria are statistical model fit measures (Gao 2006). They quantify the relative goodness of fit of various previously derived statistical models for a given sample of data.

Impulse response function (IRF)
Because the VAR model is a nontheoretically model, direct determination of the relationship among the model variables based on parameter matrices is difficult. The aforementioned relationships can be obtained by analyzing the dynamic effect on the system when an error term changes or when the model experiences shocks. Therefore, IRFs have been proposed as tools for interpreting VAR models. The impulse technique of response analysis involves representing the reaction of each variable to a shock in each equation of the system. IRFs can be applied to reflect the effect of a random distribution term on the current and future values of endogenous variables. Thus, the dynamic effect of random disturbances on endogenous variables can be characterized. The obtained characterization reflects the influences of a random disturbance on the other variables in accordance with the VAR model as well as the dynamic process of feedback to itself. IRFs with strong time characteristics can indicate the degree of response in regard to any new information generated by any system variable. The VAR model can be formulated as a triangular moving average process, which can be written in the vector form as follows: where the φ s element in row i and column j indicates the effect of a one-unit increase in innovation of the jth variable at time t(ε ij ) on the value of the ith variable at time t + s(x i, t + s ) considering all other innovations as remaining constant at all times. Bootstrapped confidence intervals are commonly generated around IRFs.

Variance decomposition
Variance decomposition examines short-run dynamic interactions between model variables. Variance decomposition is generally employed to describe relativity effects, and it indicates the contribution of each impact factor to the forecast error variance in the VAR system. Therefore, it is necessary to analyze variance decomposition results to trace any shocks to the system. Variance decomposition provides information regarding the relative importance of each random innovation in affecting the variables of the VAR model.

Results and discussion
Identifying the dynamic evolution in land use As the rapid utilization of sandy land in the ANS mainly began after 2000, LULC data for 1990 were selected here for comparison to 2000-2018 data to evaluate the land use conditions before and after sandy land utilization. As shown in Figure 2, in 1990, the main land use categories in the ANS were cultivated land, grassland, and sandy land. The other land types occupied a limited area. The area of sandy land has decreased greatly since 1990, and sandy land has been increasingly occupied by arable land and grassland. This trend has remained the same after 2000. Moreover, the area of urban and rural construction land has rapidly increased, especially after 2000, which indicates that human activities have increased notably in the study area. These human activities are concentrated in the marginal zone of the sandy land distribution, which further verifies that the scale of sand land utilization measures is gradually expanding.
The transition matrix, which is a classic tool for LULC change detection, was employed in this study to explore the LULC changes in the ANS. An extended transition matrix was determined for each pair of compared datasets. Maps reflecting the initial and subsequent periods were overlaid to generate a matrix indicating the categorical transition in LULC occurring over time (Pontius et al. 2004;Manandhar et al. 2010;Li et al. 2010). To better display transfer data, a transfer matrix of the land use was obtained in ArcGIS software and visualized as a Sankey diagram (Figure 3) with the data flow direction. The Sankey diagram, which is also referred to as the Sankey energy balance diagram, is a specific type of flow chart in which the width of the extended branch corresponds to the size of the data flow. This diagram was divided into five periods from 1990 to 2018. Figure 3 indicates the transfer degreed of the different types of land use during each period. The highest frequency of data conversion was obtained for the 2015-2018 period, followed by the 1990-2000 period.
It should be noted that the highest conversion frequency between the various types of land use occurred from 2015-2018 (mainly grassland, farmland, and sandy land). However, a certain transfer balance existed, and the final change proportion was not significant. The final change proportions for grassland, farmland, and sandy land were 1.3%, −2.7%, and −8%, respectively. Of all the land types, only the area of sandy land area continually decreased. Since 1990, the sandy land area in the ANS decreased from 5909.2 to 4114.8 km 2 , for a total decrease of 30.4%. The main transferred-out sandy land type was grassland before 2010 and farmland after 2010. Moreover, construction land exhibited the highest increase in area, from 89.5 km 2 in 1990 to 1016.4 km 2 in 2018 (an 11.4 times increase).
As mentioned above, the main transferred-out sandy land type was grassland before 2010 and farmland after 2010. In line with the climatic conditions and rural living habits in the ANS, numerous pastures and considerable cultivated land are distributed in this area. Irrespective of whether sandy land was transferred into grassland or farmland, the local agricultural water consumption and total water consumption levels were predicted to increase. To address this issue, it is necessary to explore the relationship between sandy land area change and the total water consumption and agricultural water consumption levels.

Spatial and temporal evolution characteristics of the NDVI of 2000-2018
Remote sensing data of the monthly average NDVI from 2000 to 2018 were obtained. The NDVI exhibited a fluctuating trend over 19 years but an increasing overall trend. The NDVI increased from 0.196 in 2000 to 0.371 in 2018. The total increase in the NDVI reached 89.3%, at an average annual increase of 4.7%. After 2000, the highest obtained NDVI value was 0.594. This value was obtained in 2013. The annual maximum NDVI value was consistent with the overall trend of the annual average NDVI value. The annual maximum and annual average NDVI values exhibited growth trends. The annual average NDVI value exhibited a relatively stable growth, with a slope of 0.0092. The maximum NDVI value exhibited a higher volatility and a more significant growth trend than did the annual average NDVI value.
The slope of the trend of the maximum NDVI value was 0.0132. The average annual NDVI value of the entire ecotone reached 0.279.
In terms of the NDVI distribution within year, the vegetation coverage increased significantly after May. The vegetation coverage from May to October was significantly higher than that during other periods of the year. In addition, the average monthly NDVI reached a maximum value of 0.473 in August, which was consistent with the growth period of vegetation and crops across the ANS. Subsequently, the vegetation gradually withered and the NDVI declined. The lowest NDVI value of 0.160 was observed in February.
Due to differences in natural conditions and geographical characteristics, the vegetation coverage exhibited an uneven distribution in space and time (Figure 4). Regions with high NDVI values were concentrated in the southeastern ANS. The vegetation coverage in the northwest was relatively low, which is consistent with the distribution of loess and desert from southeast to northwest. Hence, the vegetation growth was poor in regions close to sand-and water-deficient areas. Moreover, the NDVI in the middle of the ANS was lower than that in the other regions. Ultimately, the evolution of the NDVI indicated that the vegetation coverage in the entire ecotone gradually improved, which facilitates ecological improvement of the ecotone.

Analysis of stationarity and cointegration testing
Stationarity of standardized time series of the total water consumption (TC), agricultural water consumption (AC), sandy area (SA), and NDVI (NV) of the study area was examined. To avoid data fluctuation and eliminate possible heteroscedasticity, the data were first subjected to logarithmic processing. The obtained data sequence was then subjected to ADF testing. According to the ADF test results (Table 1), the data sequence remained stable if the probability value was lower than 0.05 (95% confidence interval) and the ADF statistic was lower than the critical values at the 1%, 5%, and 10% significance levels. According to Table 1, all data sequences were nonstationary in the original sequence state. The probability values of log(TC), log(AC), and log(NV) were lower than 0.05 after the first-order difference. Hence, these sequences reached the stable state. Only the log(SA) sequence reached the stable state after the second-order difference. Considering that only same order stationarity may prevent pseudo regression, the second-order difference was applied to all data sequences to obtain a stationary sequence.
After differential processing to obtain stable time series, data sequences may lack long-term information, which is essential for problem analysis. The cointegration test is performed to assess whether the causal relationship described by a regression equation encompasses pseudo regression, i.e., to determine whether a stable relationship exists among the variables. In the Johansen cointegration test, maximum likelihood estimation is performed to analyze the cointegration relationships among multiple variables. Table 2, when the number of cointegration equations was 0, the trace statistics and maximum eigenvalue statistics were 107.8306 and 55.35650, respectively, which are greater than the critical value at a significance level of 5%. Therefore, the original hypothesis of "zero cointegration relationships" is rejected. The results of the trace statistics test and Max-Eigen statistics test indicated that cointegration relationships, namely, long-term stable equilibrium relationships, existed among the sequences. These relationships were not affected by short-term fluctuations.

Construction of the VAR model
An important aspect of the VAR model is the determination of a suitable lag length. If the lag length is too short, the dynamic relationship between variables cannot be fully revealed.  Moreover, if the lag length is too long, the degrees of freedom can decrease, and the validity of the estimated model parameters may be affected. Therefore, a balance must be determined between the lag period and degrees of freedom. Different information criteria, including the final prediction error (FPE), Akaike information criterion (AIC), Schwartz information criterion (SC), Hannan-Quinn(HQ) criterion, and likelihood ratio (LR) test, were adopted to calculate the lag length of the VAR model. The test results (Table 3) indicate that the maximum lag recommended by these information criteria is 2.
The VAR model equation for runoff obtained in E-views software can be written as follows: Analysis of the impulse response and relative contribution The relationships among the variables were investigated by determining whether shocks in one variable were transmitted to the other variables. The aforementioned investigation was performed through impulse response analysis within the context of the VAR model. The results in Figure 5 reveal that one of the four variables (△ ′ log(TC), △ ′ log(AC), △ ′ log(SA), and △ ′ log(NV)) imposed a significant initial effect on the other variables. The initial effects of △ ′ log(SA) and △ ′ log(NV) on the other variables were marginally stronger and yielded longer lag lengths than those of △ ′ log(SA) and △ ′ log(NV), respectively, on the other variables ( Figure 5). This finding indicates that the relationships of the sandy land area and vegetation coverage with water consumption were significant. With decreasing sandy land conversion into vegetation cover, the increase in vegetation coverage is the direct factor causing an increase in water consumption. The change in sandy land area is an indirect factor. These results reveal that the response of water consumption to the change in NDVI is relatively stronger.
Furthermore, △ ′ log(TC), △ ′ log(AC), and △ ′ log(SA) responded to a positive △ ′ log(NV) shock beginning with the horizontal line and after one step for each point of the sample. Similarly, the responses of △ ′ log(TC) and △ ′ log(AC) to a negative △ ′ log(SA) shock exhibited a lag of one period. To a certain extent, water consumption and sandy land changes can be considered to respond to vegetation growth in the second year of growth. In other words, vegetation impacted water consumption and sandy land after a certain period of growth and accumulation, which exhibits lag characteristics.
The responses of sandy land to the other variables were highly significant. No signs of decline were observed over 10 periods, which verified that the responses of sandy land to the other variables, especially water consumption (△ ′ log(TC)and△ ′ log(AC)), were sensitive and exhibited long-term stability. The response of △ ′ log(SA) to △ ′ log(NV) exhibited hysteresis loops and an increasing trend over time, which may suggest a cumulative effect. These results indicated that the changes in the sandy land area may be negatively affected by the NDVI. Thus, the longer the NDVI continues to increase, the greater its effect on sandy land area transformation. Variance decomposition was employed to identify the proportion of the variance in one variable caused by the innovations in other variables in the vector autoregressive (VAR) system context. The variance decomposition results obtained in this study are listed in Table 4. Relatively, the NDVI yields a strong impact on the total and agricultural water consumption levels, which are approximately 8% and 9% respectively, and reveals an increasing trend over time. This also confirms the results of IRF, which indicate that the increase in vegetation coverage directly affects water consumption.
It should be noted that the sandy land area and NDVI are closely related to agricultural water consumption, as indicated in Table 4 and gradually increase over time. The arable land and grassland accounted for 78.7% of the total area of ANS, and that of sandy land is 12.3%. According to the previous analysis of the transformation of land use types, the land use type of sandy land utilization mainly turns into farmland or grassland. It is proposed that the process of sandy land transformation into the other land use types significantly impacts agricultural water consumption.

Summary of the results
From the above, the results indicate that the following: to 0.371 in 2018 at a rate of 89.3%. However, its spatial distribution was uneven. Especially along the margin of the sandy land distribution, the intensity of human disturbances was high, while the vegetation coverage was low. Regarding further sand control measures, the ecological vulnerability should be scientifically considered.
(3) There occurred long-term stable equilibrium and causal relationships among water consumption, sandy land and NDVI. Relatively, the NDVI exerted a notable impact on water consumption and exhibited an increasing trend over time. The increase of NDVI may be the direct factor causing the increase in water consumption, and the change in the sandy land area may be an indirect factor.
It could be concluded that the development and utilization measures of sandy land in the ANS area may involve a large amount of water consumption. Thus, vegetation cover growth is based on investment in water resources. This phenomenon may serve as a reminder that while various measures are implemented to combat desertification, the increase in water consumption should be considered, which may further impact areas with on water shortages. Therefore, the scientific management of water resources is an effective strategy for the transformation of sandy land and accomplishment of vegetation growth. To achieve sustainable development of the local natural system, a balance must be realized among water resources, sandy land area, and vegetation coverage.
Author contribution MZ and YW designed the concept and developed the methodology in this article and selected the case study, together with the other coauthors; SL prepared the initial draft manuscript and analyzed the results under supervision of YW and PZ.; HL and RL made the final editing and revision work. All authors read and approved the manuscript for submission.

Declarations
Ethics approval and consent to participate Not applicable.
Consent for publication Not applicable.