2.2.3. Terrestrial Water Storage

As a key variable in hydrological cycle, terrestrial water storage value was obtained from the GRACE satellite and designed mainly to observe gravity field changes with time [35]. Monthly terrestrial water storage data (Unit: cm) from 2002 to 2014 was produced by the Jet Propulsion Laboratory (JPL) with a spatial resolution of 0.5 degrees (https://grace.jpl.nasa.gov).

JPL data were the results computed from GRACE Level-1 data by mascon method [36]. In the process of calculation, the parameters required in the model, e.g., C20 term, geocentric first-order term value and post-ice rebound correction coefficien were calculated by relevant methods [37] and models [38,39]. Furthermore, based on constraints, the area of the world is divided into 4551 spherical caps with equal areas for calculation to reduce the measurement error. In order to improve spatial resolution of the data, CLM 4.0 hydrological model and coastline resolution refinement (CRI) filtering method are used to recover the signal of the solution from the mascon model and separate data from the land and sea for producing spatial distribution grid with 0.5◦ resolution [40]. While previous studies showed that GRACE satellite data still have some shortcomings, such as lower spatial and temporal resolutions, these data are still widely used because they can simulate water resources, including groundwater, for different land types [35].

#### 2.2.4. Auxiliary Data from Global Land Data Assimilation System (GLDAS)

Due to the lack of continuous monitoring of hydrological data, GLDAS-NOAH data was used to analyze hydrological factor data [41]. Soil moisture, snow water equivalent and total canopy water storage data in this paper were selected from NOAH data products in GLDAS (https: //giovanni.gsfc.nasa.gov and https://disc.gsfc.nasa.gov/). GLDAS-NOAH data time series was from January 2002 to December 2014, with monthly temporal resolution and 0.25◦ spatial resolution. In order to be consistent with GRACE satellite data for calculating groundwater change, ArcGIS software was used to resample and transform its spatial scale to 0.5◦ . The above data units are kg/m<sup>2</sup> , and unit for soil moisture data was depth, and it was recorded at depths of 10 cm, 40 cm, 100 cm, and 200 cm, respectively.

#### 2.2.5. Normalized Difference Vegetation index (NDVI)

Third-generation global inventory modeling and Mapping Research (GIMMS ndvi3g) NDVI data set from NASA Goddard Space Center (https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v1/), which was used in this paper. The changes in vegetation were modified by NDVI data, which had spatial resolution of 0.083◦ in 15-day intervals from 2002 to 2014. The changing trend of vegetation was analyzed through the seasonal data (winter: Dec, Jan-Feb; spring: Mar-May; summer: Jun-Aug; autumn: Sep-Nov; and growing season: May-Oct).

#### *2.3. Methods*

#### 2.3.1. Determination of Terrestrial Water Storage Change (TWSC)

Data of GRACE satellite's monthly gravity model reflects the components related to the earth's static structure, which is the difference between the cell's monthly water storage and the multi-year average of the cell's water storage [42,43]. Thus, the TWSC value was determined through JPL data for a total of 153 months from April 2002 to December 2014 and subtract the average in the period between 2002 and 2010. The numerical value was used to reflect the TWSC, and the positive and negative signs are used to reflect the direction of change, representing the accumulation or loss of TWSC, respectively. However, during the commissioning and operation phase of the GRACE satellite, there were problems such as sensor performance degradation and insufficient energy supply, which resulted in poor quality of the observation data of the GRACE satellite during these two periods [44,45]. The data in thirteen months were not available during the study period, June 2002, July 2002, June 2003, January 2011, June 2011, May 2012, October 2012, March 2013, August 2013, September 2013, February 2014, July 2014 and December 2014. In this paper, in order to maintain the average seasonal cycle well, interpolation, which was the average of the values for each cell from the months either side of the missing data, was used to fill in missing data [46].

#### 2.3.2. Estimation of Groundwater Change (GWC)

Spatiotemporal changes of the GWC were obtained from the GRACE (TWSC) and GLDAS data among, monthly data represented by the monthly average values [35]. The GWC had a significant correlation coefficient with the in situ observed groundwater changes, which can be used to characterize regional groundwater conditions [45]. The water balance equation was expressed as follows:

$$\text{GWC} = \text{TWSC} - \text{GLDAS}\_{\text{(SMC} + \text{SWEC} + \text{TCWSC)}} \tag{1}$$

where GWC is groundwater storage change, TWSC is terrestrial water storage change, SMC is soil moisture change, SWEC is snow water equivalent change, TCWSC is total canopy water storage change, and GLDAS(SM+SWE+TCWS) is the sum of SMC, SWEC and TCWSC. The above data units are cm.

#### 2.3.3. Analysis of NDVI

#### Maximum-Value Composite

Changes in vegetation was analyzed through the method of maximum-value composite [47]. Based on the pixel-by-pixel data of the NDVI image from January to December each year, maximum value of a pixel was determined by the calculation, and the MVC image was then generated.

#### Analysis of Spatial Trend

The spatial trend of NDVI was analyzed through the unitary linear regression method, in which time and effect factors were independent and dependent variables, respectively. Slope of the straight line from linear regression was applied to illustrate the spatial trend of NDVI [48].

$$slope = \frac{n \times \sum\_{i=1}^{n} i \times a\_i - \left(\sum\_{i=1}^{n} i\right) \left(\sum\_{i=1}^{n} a\_i\right)}{n \times \sum\_{i=1}^{n} i^2 - \left(\sum\_{i=1}^{n} i\right)^2} \tag{2}$$

where slope is the linear tendency index, *a<sup>i</sup>* is the annual NDVI in each grid, *n* = 13 is the number of years, and *i* is the year, e.g., 2003 was the 1st year, 2004 was 2nd year, etc. Value of slope > 0 represents increasing trend; and value of slope < 0 represents decreasing trend.

#### Analysis of Hurst Index

The hurst index is an effective method for quantitatively representing the long-range correlation of time series, and it has been widely used in hydrology, economics, climatology, geology, and geochemistry [49,50]. Its basic principle is to define the mean sequence for a time sequence {*NDVI*(*t*)}, *t* = 1, 2, . . . , *n*:

$$\overline{NDVI\_{\mathbb{T}}} = \frac{1}{\pi} \sum\_{t=1}^{\tau} \text{NDVI\_{\mathbb{T}}} \; \text{\; \; \tau = 1,2,\cdots,n} \tag{3}$$

$$X\_{(t,\tau)} = \sum\_{t=1}^{\tau} \left( \text{NDVI}\_t - \overline{\text{NDVI}\_{\tau}} \right) \mathbf{1} \le t \le \tau \tag{4}$$

$$R\_{\tau} = \max\_{1 \le t \le \tau} \mathcal{X}\_{(t,\tau)} - \min\_{1 \le t \le \tau} \mathcal{X}\_{(t,\tau)} \; \pi = 1, 2, \dots, n \tag{5}$$

$$S\_{\pi} = \left[\frac{1}{\pi} \sum\_{t=1}^{\pi} (\text{NDVI}\_{t} - \text{NDVI}\_{t})^{2}\right]^{\frac{1}{2}} \pi = 1,2,\cdots,n\tag{6}$$

where τ is the number of elements, *t* is the time step the year, *NDVI*<sup>τ</sup> is the time series of *NDVI*, *X*(*t*, τ) is the cumulative deviation, *R*τ is the extreme deviation sequence, and *S*τ is the standard deviation.

Taking the ratio of *R*(τ) and *S*(τ), we arrive at the following:

$$
\log\left(\frac{R}{S}\right)\_n = H \times \log(n)\tag{7}
$$

The *H* value determines whether the *NDVI* sequence is completely random or persistent. There are three indications according to the value of *H* index. Value of 0.5 < *H* < 1 indicates that *NDVI* time series is a continuous sequence, i.e., the change in the future would maintain the same trend with the past change trend, and the closer the *H* is to 1, the stronger the persistence. Value of *H* = 0.5 indicates that time series is a random sequence and there would not be long-term correlation. Value of 0 < *H* < 0.5 indicates that future change trend would be opposite to the past change trend, and the closer the *H* is to 0, the stronger the anti-persistence.

#### Analysis of Trend Test

The *F* test was used for the trend significance test. Significance test only represented the confidence level of the changing trend, regardless of change speed. Statistic calculation formula is as follows:

$$F = \mathcal{U} \times \frac{n-2}{Q} \tag{8}$$

where *U* is the error sum of the squares, *Q* is the regression sum of squares, and *n* is the number of years. According to the test results, the trend was divided into five levels, i.e., extremely significant decrease (slope < 0, *P* < 0.01), significant decrease (slope < 0, 0.01 < *P* < 0.05), non-significant change (*P* > 0.05), significant increase (slope > 0, 0.01 < *P* < 0.05), and extremely significant increase (slope > 0, *P* < 0.01).

#### **3. Results**

#### *3.1. Changes of Hydrological Factors over Time*

The observations in Figure 2 showed that hydrological factors changed annually. In general, total annual precipitation (P) and evapotranspiration (ET) increased by 0.05 cm/yr. and 0.01 cm/yr., respectively, on average from 2002 to 2014. P showed a slightly greater rate of increase after 2008 than before. During 2002 and 2014, mean annual temperature (T) decreased on average by −0.05 ◦C/yr. Terrestrial water storage change (TWSC) increased from 2002 to 2005 and decreased significantly from 2006 to the beginning of 2012, which led to a value that was approximately 6 cm less than the mean of the whole obtained TWSC series (Figure 2b). TWSC increased again after 2012 and reached a value equivalent to the average of the whole series' mean. Monthly soil moisture change (SMC), snow water equivalent change (SWEC), total canopy water storage change (TCWSC) obtained from GLDAS, fluctuated from 2002 to 2014 (Figure 2c). SMC in the study area had obvious seasonal variations and ranged from −5.24~8.26 cm, and this parameter was sensitive to changes in regional water resources. There were significant differences in the time series of SWEC, which generally reached a peak in December. In the process of freezing and thawing in spring, the value decreased slowly until the end of April. The maximum value of SWEC was 1.36 cm in December 2012. The TCWSC value reached the maximum in July or August in summer and showed an upward trend in the study area, with a range from −0.01 to 0.01 mm. Analysis of TWSC and auxiliary data-Global Land Data Assimilation System (GLDAS = SMC + SWEC + TCWSC) showed that temporal patterns groundwater change (GWC) could be described based on time series analysis (Figure 2d). GWC in the study time showed a clear decreasing trend at rate of −0.43 cm/yr., and especially experienced a significant decrease from 2007 to 2012 at a rate of −0.99 cm/yr. Furthermore, TWSC lagged behind other factors at the time scale, e.g., P, SMC and so on in Figure 2.

*Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 8 of 19

**Figure 2.** Changes of hydrometeorological factors in XRB from 2002~2014 (**a**) P, T & ET; (**b**) TWSC; (**c**) SWEC, TCWSC & SMC; (**d**) GWC. **Figure 2.** Changes of hydrometeorological factors in XRB from 2002~2014 (**a**) P, T & ET; (**b**) TWSC; (**c**) SWEC, TCWSC & SMC; (**d**) GWC.

In addition to annual changes, hydrological factors changed seasonally as well in Figure 3. P, ET and T values showed similar seasonal change patterns, with a high peak values occurring in summer and a low peak value occurring in winter. TWSC and GLDAS increased from spring to summer and then decreased from autumn to winter. A similar change trend was observed for TWSC and GLDAS with P, T and ET, indicating that there was an interaction between these factors. However, a lag effect was also observed, which required further quantitative analysis. The inter-annual trend of GWC showed a continuous downward trend, gradually decreasing from the maximum value in spring to the minimum value in winter, indicating that the annual groundwater storage gradually decreased, especially in summer where this phenomenon was more obvious. In addition to annual changes, hydrological factors changed seasonally as well in Figure 3. P, ET and T values showed similar seasonal change patterns, with a high peak values occurring in summer and a low peak value occurring in winter. TWSC and GLDAS increased from spring to summer and then decreased from autumn to winter. A similar change trend was observed for TWSC and GLDAS with P, T and ET, indicating that there was an interaction between these factors. However, a lag effect was also observed, which required further quantitative analysis. The inter-annual trend of GWC showed a continuous downward trend, gradually decreasing from the maximum value in spring to the minimum value in winter, indicating that the annual groundwater storage gradually decreased, especially in summer where this phenomenon was more obvious.

**Figure 3.** Seasonal P, T, ET, TWSC, GLDAS(SWEC&SMC&TCWSC) and GWC from 2002–2014. **Figure 3.** Seasonal P, T, ET, TWSC, GLDAS(SWEC&SMC&TCWSC) and GWC from 2002–2014. **Figure 3.** Seasonal P, T, ET, TWSC, GLDAS(SWEC&SMC&TCWSC) and GWC from 2002–2014.

#### *3.2. Spatial Distribution of Annual Hydrological Factors 3.2. Spatial Distribution of Annual Hydrological Factors 3.2. Spatial Distribution of Annual Hydrological Factors*

#### 3.2.1. Situ Observation of Hydrological Factors 3.2.1. Situ Observation of Hydrological Factors

3.2.1. Situ Observation of Hydrological Factors The in-situ observations of hydrological factors between 2002 and 2014 showed that the mean annual P increased from northwestern to southeastern areas, with an average value of 37.2 cm, while the mean annual T decreased from southwestern to northeastern areas, with an average value of 7.5 °C (Figure 4a,b). This result indicated that hydrological and meteorological factors in the study area have significant spatial heterogeneity, while spatial distributions between P and T were different. The in-situ observations of hydrological factors between 2002 and 2014 showed that the mean annual P increased from northwestern to southeastern areas, with an average value of 37.2 cm, while the mean annual T decreased from southwestern to northeastern areas, with an average value of 7.5 ◦C (Figure 4a,b). This result indicated that hydrological and meteorological factors in the study area have significant spatial heterogeneity, while spatial distributions between P and T were different. The in-situ observations of hydrological factors between 2002 and 2014 showed that the mean annual P increased from northwestern to southeastern areas, with an average value of 37.2 cm, while the mean annual T decreased from southwestern to northeastern areas, with an average value of 7.5 °C (Figure 4a,b). This result indicated that hydrological and meteorological factors in the study area have significant spatial heterogeneity, while spatial distributions between P and T were different.

**Figure 4.** Spatial distribution of P (**a**) and T (**b**) from 2002~2014 (Grid cell size: 1 km). **Figure 4.** Spatial distribution of P (**a**) and T (**b**) from 2002~2014 (Grid cell size: 1 km). **Figure 4.** Spatial distribution of P (**a**) and T (**b**) from 2002~2014 (Grid cell size: 1 km).

#### 3.2.2. Hydrological Factors from Satellite Data 3.2.2. Hydrological Factors from Satellite Data 3.2.2. Hydrological Factors from Satellite Data

TWSC data observed from the satellite in the study area showed a decreasing trend from northeast to southwest on the whole in Figure 5. The northern and central region were in a state of accumulation, while the southern region was in a state of deficit with a relatively obvious decreasing trend (Figure 5a). Based on GLDAS data, SMC, SWEC and TCWS showed equivalent water increases, as shown in Figure 5b. Spatial distribution of GLDAS was similar to that of TWSC, indicating that TWSC data observed from the satellite in the study area showed a decreasing trend from northeast to southwest on the whole in Figure 5. The northern and central region were in a state of accumulation, while the southern region was in a state of deficit with a relatively obvious decreasing trend (Figure 5a). Based on GLDAS data, SMC, SWEC and TCWS showed equivalent water increases, as shown in Figure 5b. Spatial distribution of GLDAS was similar to that of TWSC, indicating that TWSC data observed from the satellite in the study area showed a decreasing trend from northeast to southwest on the whole in Figure 5. The northern and central region were in a state of accumulation, while the southern region was in a state of deficit with a relatively obvious decreasing trend (Figure 5a). Based on GLDAS data, SMC, SWEC and TCWS showed equivalent water increases, as shown in Figure 5b. Spatial distribution of GLDAS was similar to that of TWSC, indicating that there was an

there was an interaction between these parameters. The GWC value showed a slight surplus state in

0.5°).

3.3.1. Temporal Variation

interaction between these parameters. The GWC value showed a slight surplus state in most area of XRB, which ranged from −1.2 to 1.0 cm (Figure 5c). *Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 10 of 19

**Figure 5.** Spatial distribution of TWSC (**a**), GLDAS (**b**) and GWC (**c**) from 2002~2014 (Grid cell size: **Figure 5.** Spatial distribution of TWSC (**a**), GLDAS (**b**) and GWC (**c**) from 2002~2014 (Grid cell size: 0.5◦ ).

*3.3. Spatiotemporal Variations of Normalized Difference Vegetation Index (NDVI)* 

#### *3.3. Spatiotemporal Variations of Normalized Di*ff*erence Vegetation Index (NDVI)*

#### 3.3.1. Temporal Variation

Similar to other areas in northern China, NDVI had apparent seasonal variation, and showed an upward trend during 2002 and 2014 in Figure 6. The NDVI increased from May to the most vigorous growing season in July and August, and it then decreased to the minimum level in November and remained almost unchanged until the following April, with a range of the NDVI of 0.16~0.19 (Figure 6a). In addition, NDVI showed a single peak change during the year in Figure 6a as well as hydrological factors in the study area in Figure 2a, meaning that there were relationships between them. Furthermore, the annual value of NDVI increased slightly at a rate of 0.001 yr−<sup>1</sup> from 2002 to 2014 in Figure 6b. The analysis of the cumulative anomaly method showed that NDVI decreased from 2002 to 2009 and then increased again (Figure 6b). This change trend was also consistent with the change trend of hydrological factors in the study area (Figure 2b,d), indicating that water factors in semiarid areas, especially groundwater, may be the main factors affecting vegetation changes. *Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 11 of 19 remained almost unchanged until the following April, with a range of the NDVI of 0.16~0.19 (Figure 6a). In addition, NDVI showed a single peak change during the year in Figure 6a as well as hydrological factors in the study area in Figure 2a, meaning that there were relationships between them. Furthermore, the annual value of NDVI increased slightly at a rate of 0.001 yr−1 from 2002 to 2014 in Figure 6b. The analysis of the cumulative anomaly method showed that NDVI decreased from 2002 to 2009 and then increased again (Figure 6b). This change trend was also consistent with the change trend of hydrological factors in the study area (Figure 2b,d), indicating that water factors in semiarid areas, especially groundwater, may be the main factors affecting vegetation changes.

**Figure 6.** Variations (**a**) and curve of accumulation (**b**) of NDVI with time. **Figure 6.** Variations (**a**) and curve of accumulation (**b**) of NDVI with time.

#### 3.3.2. Spatial Variations 3.3.2. Spatial Variations

Figure 7 showed that the spatial distribution of NDVI was heterogeneous in the area. NDVI value was lower in the southwestern area, with a range of 0.0~0.33, than in the central belt from west to east (Figure 7a). NDVI was greater along the rivers, which ranged from 0.63~1.00, than other parts, which ranged from 0.0~0.63. Furthermore, the annual value of NDVI decreased over 63.2% of the area at a rate of −0.008~0.0 a−1, while 36.8% of the area showed an increase in the NDVI at rate of 0~0.02 a−1 (Figure 7b). F test results suggested that significant variance area accounted for 89.3% of the whole area and was mostly distributed in areas with a relatively lower NDVI value (Figure 7c). Significant variance was not observed in cultivated areas along rivers. The Hurst index (0.06~0.97) indicated that in 55.7% of the area, particularly in areas with a decreasing NDVI trend, the variance trend was continuously maintained for a long time while the remaining 44.3% of the area would show fluctuations (Figure 7d). Figure 7 showed that the spatial distribution of NDVI was heterogeneous in the area. NDVI value was lower in the southwestern area, with a range of 0.0~0.33, than in the central belt from west to east (Figure 7a). NDVI was greater along the rivers, which ranged from 0.63~1.00, than other parts, which ranged from 0.0~0.63. Furthermore, the annual value of NDVI decreased over 63.2% of the area at a rate of <sup>−</sup>0.008~0.0 a−<sup>1</sup> , while 36.8% of the area showed an increase in the NDVI at rate of 0~0.02 a−<sup>1</sup> (Figure 7b). F test results suggested that significant variance area accounted for 89.3% of the whole area and was mostly distributed in areas with a relatively lower NDVI value (Figure 7c). Significant variance was not observed in cultivated areas along rivers. The Hurst index (0.06~0.97) indicated that in 55.7% of the area, particularly in areas with a decreasing NDVI trend, the variance trend was continuously maintained for a long time while the remaining 44.3% of the area would show fluctuations (Figure 7d).

*Remote Sens.* **2020**, *12*, x FOR PEER REVIEW 12 of 19

**Figure 7.** Spatial distribution of annual NDVI value (**a**), annual variation trend (**b**), F test results (**c**), and Hurst index (**d**) from 2002~2014 (Grid cell size: 8 km). **Figure 7.** Spatial distribution of annual NDVI value (**a**), annual variation trend (**b**), F test results (**c**),and Hurst index (**d**) from 2002~2014 (Grid cell size: 8 km).

#### *3.4. Relationships between Hydrological Factors 3.4. Relationships between Hydrological Factors*

The value of P was highly related to ET, with a correlation coefficient of 0.685 (*P* < 0.01). TWSC had a significant correlation (*R* = 0.680, *P* < 0.05) with GWC, as compared to GLDAS in Table 2, indicating that TWSC mainly consisted of GWC in the area. GLDAS had a significant and positive correlation with P (*R* = 0.634, *P* < 0.05) and ET (*R* = 0.686, *P* < 0.01), meaning that a higher P rate corresponded to more SMC. Significant positive correlation was observed between T and GWC (*R* = 0.680, *P* < 0.05), which indicated that water demand of plants increased as T increased in XRB and groundwater was the main source of water in semiarid area and changed significantly. In summary, there are complex correlations among regional hydrological elements. The value of P was highly related to ET, with a correlation coefficient of 0.685 (*P* < 0.01). TWSC had a significant correlation (*R* = 0.680, *P* < 0.05) with GWC, as compared to GLDAS in Table 2, indicating that TWSC mainly consisted of GWC in the area. GLDAS had a significant and positive correlation with P (*R* = 0.634, *P* < 0.05) and ET (*R* = 0.686, *P* < 0.01), meaning that a higher P rate corresponded to more SMC. Significant positive correlation was observed between T and GWC (*R* = 0.680, *P* < 0.05), which indicated that water demand of plants increased as T increased in XRB and groundwater was the main source of water in semiarid area and changed significantly. In summary, there are complex correlations among regional hydrological elements.


**Table 2.** Correlation between hydrological factors at annual scale. **Table 2.** Correlation between hydrological factors at annual scale.

GWC 1 \*\*. Correlation is significant at the 0.01 level (2-tailed). \*. Correlation is significant at the 0.05 level (2-tailed). \*\*. Correlation is significant at the 0.01 level (2-tailed). \*. Correlation is significant at the 0.05 level (2-tailed).

P is not the only one of the most important factors underlying the interactions among the hydrological cycle but also the key variable in the water resources in semiarid areas. It accounted for approximately 75% of total water resource recharge in normal years and approximately 57% of recharge in dry years in study area [11]. A lag time of P infiltration is observed based on GLDAS and P is not the only one of the most important factors underlying the interactions among the hydrological cycle but also the key variable in the water resources in semiarid areas. It accounted for approximately 75% of total water resource recharge in normal years and approximately 57% of recharge in dry years in study area [11]. A lag time of P infiltration is observed based on GLDAS and GWC. Lag time is usually represented in hydrological forecasts based on P in the hydrological

equilibrium [51]. Table 3 showed their relationships between TWSC, GLDAS and GWC and P and ET, which were analyzed to determine lag times. The results showed that TWSC and GLDAS had significant correlation (*P* < 0.05) with P at different lag periods, i.e., a one-month time lag occurred (*R* = 0.312, *P* < 0.05). Moreover, these parameters were significantly and positively correlated with the ET in the current month (*R* = 0.276, *P* < 0.05) and one month prior (*R* = 0.276, *P* < 0.05) but not positively correlated (*p* > 0.01) with the ET two to three months later. This finding suggested that the P and ET in the current month and one-month prior could significantly affect the TWSC and GLDAS. Meanwhile, P and ET did not have correlation (*P* > 0.05) with GWC, which might be related to the considerable depth of groundwater table due to overexploitation, long time period for infiltration to recharge groundwater or lack of infiltration recharge of groundwater due to ET. As a result of this, P and ET were increasingly less sensitive to GWC.


**Table 3.** Correlation between hydrological factors at monthly scale in a different lag period.

\*\*. Correlation is significant at the 0.01 level (2-tailed). \*. Correlation is significant at the 0.05 level (2-tailed).

## *3.5. Relationships between Hydrological Factors and NDVI*

To obtain better understanding of relationships between hydrological factors and NDVI, correlation analysis was carried out. As shown in Table 4, the results showed that ET and GLDAS were significantly and positively correlated with NDVI (*R* = 0.747 and 0.704; *P* < 0.01) in the growing season (May-Oct) of vegetation. Correlation coefficient between NDVI and GWC in the growing season was significantly negatively correlated (*R* = −0.64, *P* < 0.05), indicating that groundwater represents an important water source for vegetation growth in the XRB. T and TWSC had no effect on the growth of vegetation in each season, and their correlation coefficients did not pass the significance test. In summary, ET, GLDAS and GWC could be the major hydrological parameters that affect vegetation dynamics in the growing season.

**Table 4.** Correlation between NDVI and hydrological factors in different seasons.


\*\*. Correlation is significant at the 0.01 level (2-tailed). \*. Correlation is significant at the 0.05 level (2-tailed).

#### **4. Discussion**

#### *4.1. Dynamics of Hydroecological Elements in Semiarid Area*

Water for agricultural irrigation relies heavily on groundwater due to shortages of surface water in study area [47], and reports have indicated that more than 80% of the water use was from groundwater [11]. Spatiotemporal variation and distribution characteristics of the hydrological factors were determined by multiple satellite data, which have been widely used as effective approaches for detecting the hydrological dynamics in ecotones in semiarid areas. Combined with the analysis results of climate change in the study area, spatial distribution of precipitation (P) and temperature (T) was uneven, as shown in Figure 4 and presented a decreasing and increasing trend from the southeast to the northwest, respectively. The results showed that, on the spatial scale, terrestrial water storage change (TWSC) results based on JPL satellite data and auxiliary data-Global Land Data Assimilation System (GLDAS = SMC + SWEC + TCWSC) values were also closely related to P and T. Application of the multi-satellite data would represent an effective approach to monitoring and modeling the dynamics of groundwater change (GWC). Our result confirms that the environment changes have occurred in the semiarid area and that they can be related to dynamics of hydrogeological elements were easily observable and measurable through satellite images. Through correlation analysis, there are complex relationships among regional hydrological elements in Table 2. Furthermore, soil moisture change (SMC) and snow water equivalent change (SWEC) obtained from the GLDAS data played an important role in the hydrological cycle in the area. The SWEC not only affected the surface runoff, groundwater recharge, and SMC in the spring [52–54], but also indicated the snow cover in the winter and the accumulated air temperature in the spring (Figure 2). SMC could even directly reflect the vegetation growth in growing season, and predict the GWC [55,56]. SMC was higher for irrigated land than grassland in the vegetation growth season (Figure 5). In particular, SMC and GWC directly affected the local hydrologic cycle, which could be used to optimize the management and utilization of water resources and improve vegetation growth. However, compared to other areas [57], the response between P and GWC was increasingly insensitive in the study area with a lag period of one month. There are many reasons for this phenomenon, but anthropic factors [58], e.g., irrigation and overexploitation of groundwater, have thickened the vadose zone of soil, and the relationship between P and groundwater, which is is becoming more and more complex.

While the application of the satellite data can significantly improve the level of water resources assessment in this area, there are still some deficiencies, such as time scale (monthly), spatial resolution (0.5 degree) and data accuracy verification, which have a significant impact on the research results of small and medium-sized watersheds, e.g., the correlation strength between them and other environmental molecules is very poor, only the correlation can be considered. Therefore, how to combine the limited water resources monitoring information with new methods, new technologies and new achievements to carry out effective regional water resources management has become a practical problem demand of relevant departments of water resource management.

#### *4.2. Impacts of Hydrological Factors on Vegetation*

Hydrological factors played a controlling role in terrestrial ecosystems [59,60]. Field survey in July of 2015 and August of 2016 and normalized difference vegetation index (NDVI) value showed that the northern part of the study was dominated by grassland with medium vegetation coverage, the middle part was dominated by cultivated land with high vegetation coverage, and the southern part was dominated by sandy land with low vegetation coverage. The results in the paper showed that ET, GLDAS and GWC could be the major hydrological parameters that affect the vegetation dynamics in the growing season. The northeastern XRB included grassland with some small rivers with high terrain and vegetation coverage, and it was affected by snow melt water in the spring as well as rainfall. As a result, GLDAS and TWSC showed a cumulative trend with higher values in the northern piedmont plain than in other areas, indicating that the piedmont plain had basically maintained its original ecology and was rarely affected by human activities. Change of GWC in the central part of the study area, which is a flat, wide area of cultivated land, was relatively reduced as shown in Figure 5c. Therefore, irrigation had a significant impact on regional SMC and even affected regional water resource reserves [61], which indicated that the GWC might be affected by human activities, such as the expansion of cultivated land area [62] and overexploitation and utilization of groundwater [63]. TWSC in the southwest of the study area showed a decreasing trend (Figure 5a) [64]. Major land use types in the area consisted of typical steppe, meadow steppe, and cropland in Figure 1, and vegetation growth in these areas was dependent on ET (*R* = 0.747, *P* < 0.05) and GLDAS (*R* = 0.704, *P* < 0.05) which was mainly consisted of SMC. However, in summer and the vegetation growing season, P and NDVI were positively correlated, although the occurrence of this phenomenon is mainly related to the serious desertification of the local surface and the inability of the soil to store water correlation coefficient failed the significance test (Table 3) which indicated that lower P did not have a significant effect on the vegetation in the semiarid area [65]. NDVI, as an intuitive indicator of vegetation growth, plays a fundamental role in reflecting the characteristics of vegetation growth and distribution. However, NDVI is hardly enough to explain vegetation changes, such as vegetation height, density and so on. Furthermore, grassland will often give a stronger NDVI value than forest but the effect on ET and GWC are very different. In future, more indexes related to water should be used in the study of regional hydroecological processes.

#### **5. Conclusions**

Using multiple satellite data, the spatiotemporal changes of hydrological elements in semiarid areas from 2002 to 2014 and their effects on vegetation were analyzed in this paper. The main conclusions are as follows.

(1) Analysis of spatiotemporal characteristics of hydrological elements showed that the hydrological process in the study area has changed significantly. Annual variations of precipitation (P), evapotranspiration (ET) and temperature (T) were all in the form of a single peak, and the interannual variation law was slightly different, with rates of 0.05 cm/yr., 0.01 cm/yr. and −0.05◦C/yr., respectively. Terrestrial water storage change (TWSC) showed a fluctuating trend, with initial increase, then decrease, and finally an increase, which was consistent with the change trend of soil moisture change (SMC). Groundwater change (GWC) showed a decreasing trend at a rate of −0.17 cm/yr.

(2) Complex correlations occurred among regional hydrological elements. P and ET were significantly correlated (*R* = 0.685, *P* < 0.01) and had a greater impact on the GLDAS (*R* = 0.634, *P* < 0.05 and *R* = 0.686, *P* < 0.01) than on the TWSC and GWC. GWC is an important component of the TWSC in the region (*R* = 0.680, *P* < 0.05), and it was more sensitive to the T response (*R* = 0.570, *P* < 0.05). Furthermore, P would lead to greater TWSC and GLDAS values when P preceded the TWSC by one month, whereas smaller changes would be observed when P preceded these parameters by two months. The time lag of the GLDAS that was influenced by P was more obvious than that of the TWSC. The TWSC and GLDAS were significantly and positively correlated with the ET in the current month and one month prior and not positively correlated with the ET two to three months later. Due to overexploitation, P and ET did not have any effect on the GWC.

(3) Normalized difference vegetation index (NDVI) had obvious seasonal variations and showed an upward trend at a rate of 0.001 yr−<sup>1</sup> during 2002 and 2014. Spatial distribution of NDVI was heterogeneous in study area. NDVI decreased by 63.2% of the area at rate of <sup>−</sup>0.008~0.0 yr−<sup>1</sup> . The area showing significant variance accounted for 89.3% of the whole area, and 55.7% of the area would maintain the variance trend continuously for a long time, with these areas mainly showing decreasing NDVI change trends. However, other 44.3% area would show fluctuations. Hydrological factors play a controlling role on terrestrial ecosystems. In growing season, ET, GLDAS and GWC were the parameters that limited vegetation growth, and they were more important than other factors in XRB.

Application of satellite data could significantly improve the water assessment capability in semiarid areas and could be used for regional water resource and eco-environment management in semiarid areas. Hydrological factors, such as TWSC, SM and GWC, spatiotemporal dynamic and their correlations were successfully determined and analyzed in this study. Impacts of hydrological change on NDVI were also identified based on the analysis. These results will help to understand regional hydroecological processes, and also provide a scientific basis for local environmental management.

**Author Contributions:** Methodology, S.Z.; software, B.S.; data curation, Y.Z.; writing—original draft preparation, Y.Z.; writing—review and editing, P.L. All authors have read and agreed to the published version of the manuscript. **Funding:** This research was funded by National Key R&D Program of China (2018YFE0103800), Ministry of Water Resources Public Welfare Special Scientific Research Project—Semi-arid Zone Water Cycle and Water Ecological Security Key Technology Research (number 201501031), National Natural Science Foundation of China (51779118, 51669021, 51869020), Fundamental Research Funds for the Central Universities, CHD (300102299302, 300102299102, 300102299104), One Hundred Talent Plan of Shaanxi Province, International Collaborative Research of Disaster Prevention Research Institute of Kyoto University(2019W-02), Excellent Projects for Science and Technology Activities of Overseas Staff in Shaanxi Province (2018038), Yan'an University Project (YDBK2017-19, YDBK2019-35, YDBK2019-36), and Yan'an Science and Technology Bureau (project 2018KS-02).

**Acknowledgments:** The authors would like to thank the Editors and the anonymous reviewers for their crucial comments, which improved the quality of this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
