Watershed Drought and Ecosystem Services: Spatiotemporal Characteristics and Gray Relational Analysis

: This study explored the spatiotemporal characteristics of drought and ecosystem services (using soil conservation services as an example) in the YanHe Watershed, which is a typical water basin in the Loess Plateau of China, experiencing soil erosion. Herein, soil conservation was simulated using the Soil and Water Assessment Tool (SWAT), and the relationship between drought, soil conservation services, and meteorological, vegetation, and other factors since the implementation of the ‘Grain for Green’ Project (GFGP) in 1999, were analyzed using the gray relational analysis (GRA) method. The results showed that: (1) The vegetation cover of the Watershed has increased signiﬁcantly, and evapotranspiration (ET) increased by 14.35 mm · a − 1 , thereby increasing water consumption by 8.997 × 10 8 m 3 · a − 1 (compared to 2000). (2) Drought affected 63.86% of the watershed area, gradually worsening from south to north; it decreased in certain middle areas but increased in the humid areas on the southern edge. (3) The watershed soil conservation services, measured by the soil conservation modulus (SCM), increased steadily from 116.87 t · ha − 1 · a − 1 in 2000 to 412.58 t · ha − 1 · a − 1 in 2015, at a multi-year average of 235.69 t · ha − 1 · a − 1 , and indicated great spatial variations, with a large variation in the downstream and small variations in the upstream and midstream areas. (4) Integrating normalized difference vegetation index ( NDVI ) data into SWAT model improved the model simulation accuracy; during the calibration period, the coefﬁcient of determination (R 2 ) increased from 0.63 to 0.76 and Nash–Sutcliffe efﬁciency (NSE) from 0.46 to 0.51; and during the validation period, the R 2 increased from 0.82 to 0.93 and the NSE from 0.57 to 0.61. (5) The GRA can be applied to gray control systems, such as the ecosystem; herein, vegetation cover and drought primarily affected ET and soil conservation services. The analysis results showed that vegetation restoration enhanced the soil conservation services, but increased ET and aggravated drought to a certain extent. This study analyzed the spatiotemporal variations in vegetation coverage and the response of ET to vegetation restoration in the YanHe Watershed, to verify the signiﬁcant role of vegetation restoration in restraining soil erosion and evaluate the extent of water resource consumption due to ET in the semi-arid and semi-humid Loess-area basin during the GFGP period. Thus, this approach may effectively provide a scientiﬁc basis for evaluating the ecological effects of the GFGP and formulating policies to identify the impact of human ecological restoration on ecosystem services.


Introduction
Ecosystem services (ESs) are defined as "the natural environmental conditions and effects that are formed and maintained by ecosystems and ecological processes and upon which people rely for existence"; they are the benefits that humans gain directly or indirectly from the ecosystem [1]. Based on their contribution to human beings, ESs can be divided into four categories: Supply, regulation, support, and cultural services. Soil conservation services refer to the ability of the ecosystem to control and prevent soil erosion, and store and maintain sediments. As an important regulating service, it helps prevent regional land degradation and reduce the risk of flooding [2]. With the launch of the Millennium Ecosystem Assessment (MEA), particularly the 2003 report, the study of ESs has become a research hotspot in the field of global change. The Loess Plateau in China is located in an ecologically sensitive area, with the coexistence of drought and soil erosion [3,4]. Drought can lead to vegetation degradation [5], reduction of net primary productivity of vegetation [6], land desertification [7], and soil erosion [8] in the natural ecosystem, and affect the supply of water resources [9], agricultural production [10], and energy and power supply [11] in the socioeconomic system. Soil erosion not only causes land degradation [12] and soil fertility decline [13] and affects agricultural production and food security [14], but also causes eutrophication of water bodies [15], siltation of reservoir channels [16], aggravation of flood disasters [17], and other problems in downstream areas. Because of these ecological problems, research on ESs in ecologically fragile areas has become a focus in academic research.
Since the 1970s, the Chinese government has adopted several soil and water conservation measures, including the construction of terraces and silt dams, to mitigate severe soil erosion [18]. In 1999, the government proposed and implemented a series of large-scale ecological modification projects, such as the 'Grain for Green' Project (GFGP) [19]. Relevant studies [20][21][22] have shown that the GFGP has significantly increased vegetation coverage, alleviated soil erosion, improved the ecological environment, and enhanced functioning of ecological services in the Loess Plateau. In the global terrestrial ecosystem, vegetation restoration is a primary method to prevent soil erosion, and the study of its impact on ESs has received increasing focus. Fu et al. [23] reported that vegetation has significantly improved in typical regions of the GFGP on the Loess Plateau, and soil conservation, hydrologic regulation, and carbon sequestration services have increased by varying degrees. Wang et al. [24] studied Yulin City in particular and found that soil erosion followed a downward trend with increasing normalized difference vegetation index (NDVI) values, and its temporal and spatial variations were primarily related to rainfall, slope, and land use. Kervroëdan et al. [25] revealed that vegetation-controlled runoff and erosion in the Loess Belt in Europe, by increasing hydraulic roughness and sediment retention; the ability of vegetation to control soil erosion increased with an increase in stem density, diameter, and leaf area. Gashaw et al. [26] reported that the reclamation of cultivated land and degradation of vegetation in the Andassa Watershed in Ethiopia increased food production and pollination, but degraded several other ESs to varying levels, indicating that vegetation has an irreplaceable role in regional ecosystems. Most studies on soil erosion control have analyzed the effect of vegetation macroscopically, with few studies examining the spatial and temporal distribution characteristics and influencing factors of vegetation restoration on soil conservation services.
Tobler's First Law of Geography [27] states that "everything is related to everything else". Vegetation restoration enhances soil conservation and reduces sediment load. However, vegetation growth increases regional evapotranspiration (ET) by changing surface roughness, surface reflectivity, interception of precipitation, and transpiration [28], which further aggravates regional water shortages [29] and the degree of drought [30]. Similar results have been reported by other researchers [31,32]. Drought is a major meteorological disaster affecting socioeconomic and natural ecosystems. Drought indices are the primary tools to identify and monitor drought and investigate drought characteristics [33]. Commonly used drought indices include the standardized precipitation index (SPI), stan-dardized precipitation evaporation index (SPEI), evaporative stress index (ESI), and Palmer drought severity index (PDSI). Although the SPI index is widely accepted by climatologists and hydrologists, it is based only on accumulated precipitation and does not consider the impact of water demand on drought. To overcome this, the SPEI was developed by Vicente-Serrano et al. [34] to incorporate potential evapotranspiration (PET). Potop et al. [35] showed that the SPEI might better characterize drought in areas with significantly warm climates; however, it requires long time-series data. Most studies calculate the SPI and SPEI using data obtained from meteorological stations. Moreover, due to the inconsistent research scales, the drought rules observed at the point scale are not scientifically reliable for evaluating vegetation restoration and policymaking at the regional scale [36]. Remote sensing data provide multi-resolution (spectral and spatial) and abundant information [37]; hence, the construction of drought indices based on remote sensing data has gradually become the primary method for evaluating drought. The ESI is the ratio of actual to PET and can monitor drought caused by multiple factors, such as decreased precipitation and increased radiation [38]. The ESI has been widely used in several drought monitoring studies [39,40]. Yoon et al. [40] identified the advantages of ESI over other satellite-based drought assessment indices and explored its applicability as a new indicator of agricultural drought monitoring. Liu et al. [41] showed that the ESI indicated drought in arid river basins of Northwest China better than RDI. Therefore, in this study, the ESI was constructed to monitor drought in the YanHe Watershed.
The impact of vegetation restoration on spatiotemporal pattern changes of drought and soil conservation requires urgent investigation because of the adverse effects of drought on ecology and socio-economy, and the significant impact of vegetation restoration in enhancing soil conservation services [42]. However, the ecosystem is a gray system with complex response mechanisms between drought, soil conservation, and their respective influencing factors. Herein, we introduced the GRA method to evaluate the effect of vegetation restoration on drought and supply of ecosystem services. The ecosystem service cascade [43] indicates that the process by which humans obtain ESs from natural ecosystems [44] and is a gradual transformation of ecosystem structure and processes. Therefore, this study analyzed the process of increasing water consumption through vegetation restoration using remote sensing data, and quantified soil conservation services using the Soil and Water Assessment Tool (SWAT) model. The cascade model of processeffect-service was adopted to analyze drought and soil conservation services in the YanHe Watershed.
This study aimed to: (1) Understand the spatiotemporal variations of vegetation coverage and the response of ET to vegetation restoration in the YanHe Watershed since 2000; (2) explore the spatiotemporal characteristics and evolutionary trends of drought in the basin; and (3) evaluate spatiotemporal variations and the influencing factors of soil conservation services. The results of this study were used to analyze the excessive consumption of water resources caused by vegetation restoration in the YanHe Watershed during 2000-2015, and the significant role of vegetation restoration in restraining soil erosion. The research results are conducive to comprehensively evaluate the ecological effects of the GFGP, thereby provide a reference for ecological control on the Loess Plateau and guide the construction of ecological civilization.

Overview of the Study Area
The YanHe River is a primary tributary in the middle reaches of the Yellow River; it originates from the Baiyu Mountains in Jingbian County, Shaanxi Province, and flows through the Zhidan, Ansai, and Baota counties [45]. After 286.9 km, it flows into the Yellow River in Yanchang County. The YanHe Watershed (36 • 27 -37 • 58 N, 108 • 41 -110 • 29 E) covers 7725 km 2 and has an altitude of 461-1787 m (average elevation 1212 m) (Figure 1), with an average slope of 4.3% and gully density of 4.7 km/km 2 . The regional terrain is high in the northwest and low in the southeast, with forms including the hilly and gully regions (lower reaches), valley and plain regions (upper reaches), and residual loess platform regions (middle reaches). The study area lies in a warm temperate zone that has a continental semi-arid monsoon climate, with an average annual precipitation of 509 mm, approximately 78% of which occurs between June and October as rainstorms. The average annual temperature is 8.8-10.2 • C, average annual evaporation 1000 mm, and annual total radiation 125-134 kcal/cm 2 [46]. The main soil type is loose and soft loessial soil with poor erosion resistance. The YanHe Watershed is a part of the Loess Hilly-Gully Region, with fragmented landforms and crisscross gullies. Soil erosion is most severe in the YanHe Watershed area of the Loess Plateau, occurring in approximately 80% of the total area [47]. The GFGP, by implementing the largest set of revegetation measures in the YanHe Watershed, has changed the vegetation cover since 1999. Vegetation types vary with environmental gradients; steppes and desert steppes occur in the northwest, temperate forest-steppes in the middle, and warm temperate-deciduous broadleaf and mixed coniferous-deciduous broadleaf forests in the southeast of the watershed. 110°29′ E) covers 7725 km 2 and has an altitude of 461-1787 m (average elevation 1212 m) ( Figure 1), with an average slope of 4.3% and gully density of 4.7 km/km 2 . The regional terrain is high in the northwest and low in the southeast, with forms including the hilly and gully regions (lower reaches), valley and plain regions (upper reaches), and residual loess platform regions (middle reaches). The study area lies in a warm temperate zone that has a continental semi-arid monsoon climate, with an average annual precipitation of 509 mm, approximately 78% of which occurs between June and October as rainstorms. The average annual temperature is 8.8-10.2 °C, average annual evaporation 1000 mm, and annual total radiation 125-134 kcal/cm 2 [46]. The main soil type is loose and soft loessial soil with poor erosion resistance. The YanHe Watershed is a part of the Loess Hilly-Gully Region, with fragmented landforms and crisscross gullies. Soil erosion is most severe in the YanHe Watershed area of the Loess Plateau, occurring in approximately 80% of the total area [47]. The GFGP, by implementing the largest set of revegetation measures in the YanHe Watershed, has changed the vegetation cover since 1999. Vegetation types vary with environmental gradients; steppes and desert steppes occur in the northwest, temperate forest-steppes in the middle, and warm temperate-deciduous broadleaf and mixed coniferous-deciduous broadleaf forests in the southeast of the watershed.

Data Sources and Preprocessing
This study utilized multiple datasets, including the digital elevation model (DEM), meteorological, hydrological, soil, land use, vegetation, and ET data. (1) DEM data, with a resolution of 30 m, were downloaded from the Geospatial Data Cloud

Data Sources and Preprocessing
This study utilized multiple datasets, including the digital elevation model (DEM), meteorological, hydrological, soil, land use, vegetation, and ET data.

Gray Relational Analysis
Gray relational analysis (GRA) is a multivariate statistical analysis method that assesses the degree of correlation among factors according to their similarity or difference in development trends [48]. As the similarity among the geometric shape of the time series increases, the correlation among these factors becomes more significant. Based on this principle, the major and minor factors that affect the system's development can be determined. According to Gray System Theory, ecosystems with complicated factor relationships and unclear internal principles are gray control systems. Accounting for the interference of unknown and undetermined factors has increased our understanding of ESs [49]. Gray relational analysis is also suitable for small data samples that do not meet standard distribution laws. While the number of GRA calculations is relatively small, the results are consistent with the qualitative analysis results. Therefore, GRA can be used to analyze the impact of vegetation restoration on drought and supply of ESs. Evaluations using this method are usually based on the gray relational coefficient, degree, and sequence. The specific steps are as follows: 1.
The evaluation criterion (reference sequence) and evaluation object (comparison sequence) were determined.

2.
The difference in absolute values was reduced by unifying the data into an approximate range, focusing on data changes and trends.

3.
The gray relational coefficient, ξ i (k), was calculated as where ξ i (k) is the correlation coefficient between reference sequence y(k) and several comparison sequences x 1 , x 2 , K, and x n at each moment; ρ is the distinguishing coefficient; usually ρ = 0.5.

4.
The gray relational degree, r i , was obtained as a quantitative expression of the degree of association between the reference and comparison sequences as where r i is the gray correlation degree.

5.
A correlation sequence of the evaluation subject was established according to the calculated value of the gray relational degree.

Drought Index Based on Remotely Sensed Data
The evaporative stress index (ESI) was used to assess the degree of drought, and the data source was the MOD16 product, which has a good simulation effect on land ET. The ESI considers the role of ET in reflecting the response of vegetation and soil moisture to drought, as the ratio of ET to PET reflects water supply and demand and responds rapidly to water consumption by plants. Anderson et al. [50] standardized the ratio of actual ET to PET (f PET ), which was termed the evaporative stress index (ESI), calculated as follows: where N is the length of the data sequence, f PET the mean f (i) PET , and α fPET the standard deviation f (i) PET . A positive ESI indicates that the area is wet, while a negative ESI indicates that the area is dry; the specific correspondence is shown in Table 1.

. Soil Conservation Services Assessment
The SWAT model was used to simulate soil conservation services in the YanHe Watershed. The SWAT model, with a strong physical mechanism [51], is a long time-series watershed-scale distributed hydrological model developed by the United States Department of Agriculture-Agricultural Research Center (USDA-ARC) based on the Simulator Water Resources in Rural Basins (SWRRB) model [52]. It has been widely used in the simulation and prediction of hydrological processes [53] and is highly applicable to the Loess Plateau [54,55]. Equation (4) calculates the actual soil erosion using the Modified Universal Soil Loss Equation (MUSLE). This study assumed the absence of vegetation cover when calculating the potential soil erosion. Equation (5) presents the calculation of the amount of soil conservation. sed = 11.8· Q sur f ·q peak ·area hru 0.56 where sed is the actual amount of soil erosion (metric tons), sed c the amount of soil conservation, Q surf the surface runoff volume (mm·ha −1 ), p each the peak runoff rate (m 3 ·s −1 ), area hru the area of the hydrological response unit (HRU) (ha), K USLE the USLE (universal soil loss equation) soil erodibility factor, C USLE the USLE cover and management factor, P USLE the USLE support practice factor, LS USLE the USLE topographic factor, and CFRG the coarse fragment factor. The calibration and validation procedures were based on the SUFI-2 algorithm of the SWAT-CUP [56]. The Nash-Sutcliffe efficiency (NSE) and the coefficient of determination (R 2 ) were used to evaluate the simulation accuracy of the SWAT model.
where n is the number of observations; Q oi and Q mi are observed and predicted values at time i, respectively; Q o and Q m are the means of the observed and predicted values, respectively; NSE is the performance of the model output relative to the mean of the observed data; and R 2 is the consistency of the trend between the observed and predicted values. Generally, the model performance was considered 'satisfactory', if R 2 and NSE values were >0.5 for a monthly time-step [57]. The ArcSWAT2012 software was used to construct the SWAT model of the YanHe Watershed. Flow direction and accumulation were calculated based on DEM and the river network of the basin was extracted with 10,000 ha as the minimum catchment area threshold. The whole basin was divided into 47 sub-basins [58]; by superimposing land use, soil, and slope data, the 47 sub-basins were further divided into 1511 HRUs [59]. Finally, attribute data, such as soil and weather databases, and weather generators, were imported into the SWAT model. The first year (2000) was used as a warm-up period to mitigate the initial conditions. The SWAT parameters were calibrated using monthly sediment data from the Ganguyi station for a period of eight years (2001-2008) and then validated using monthly sediment data for seven years (2009)(2010)(2011)(2012)(2013)(2014)(2015). The model output included monthly runoff, sediment, and other data for each HRU during the simulation period. Based on previous studies [60], the soil conservation modulus (SCM) was divided into five gradients: Lower SCM (0-250 t·ha −1 ), low SCM (250-500 t·ha −1 ), medium SCM (500-750 t·ha −1 ), high SCM (750-1000 t·ha −1 ), and higher SCM (1000-1500 t·ha −1 ).

Integration of Remotely Sensed NDVI into SWAT
Under natural conditions, DEM and soil data vary slowly [61]; hence, these were not considered in this study. In comparison, human activity and vegetation restoration affected land use changes in the YanHe Watershed, but only one period of land use data was included in the original SWAT model.
The SWAT model reads the annual average C factor (C USLE,aa ) of a particular land use type from the crop database (crop.dat); calculates the minimum C factor (C USLE,mn ) of this land use; and calculates C USLE based on C USLE,mn , and soil surface residue (rsd surf ). In the long-term series simulation, the original model does not consider the impact of land use transfer on vegetation coverage changes and has some other limitations. Therefore, the study proposed recalculating C USLE using the remotely sensed NDVI to reflect real vegetation cover changes.
MOD13Q1 is a long-term serial product of the MODIS series, which can express vegetation changes in agricultural land in different seasons and farming modes without additional measurements or special satellite-data processing [62]. Considering the spatial resolution of the land use, soil, and slope data and the smallest unit of HRU as 30 × 30 m, we generated a 30 × 30 m fishnet based on the DEM raster map. Using the "Extract Multi Values to Points" tool provided by ArcGIS, the NDVI values were then extracted for every point and the average NDVI of each HRU was calculated to obtain the monthly NDVI dataset.
Van der Knijff [63] established the relationship between NDVI and C USLE based on National Oceanic and Atmospheric Administration (NOAA) Advanced Very-High-Resolution Radiometer (AVHRR) data as follows: where α and β are constants. Van Leeuwen and Sammons [64] showed that when using MODIS-NDVI data, α and β should be 2.5 and 1, respectively. Therefore, we calculated the new C USLE value based on Equation (8) also modified and NDVI data were introduced. Figure 2 shows the flowchart of the C USLE modification in the SWAT model based on the NDVI data.
where α and β are constants. Van Leeuwen and Sammons [64] showed that when using MODIS-NDVI data, α and β should be 2.5 and 1, respectively. Therefore, we calculated the new CUSLE value based on Equation (8), for α = 2.5 and β = 1. The SWAT source code was also modified and NDVI data were introduced. Figure 2 shows the flowchart of the CUSLE modification in the SWAT model based on the NDVI data.

Mann-Kendall Test
The Mann-Kendall (MK) trend test is a nonparametric test method [65] widely used for trend analysis of vegetation and hydrological series. The MK test sample does not need to conform to a certain distribution and is not affected by outliers. The τ-value in the MK test reflects the correlation between the test sequence and time, and the p-value represents the significance level of the test results. The test results were divided into seven categories according to τ-and p-values to determine the trend change (Table 2).
Decreased very significantly

Mann-Kendall Test
The Mann-Kendall (MK) trend test is a nonparametric test method [65] widely used for trend analysis of vegetation and hydrological series. The MK test sample does not need to conform to a certain distribution and is not affected by outliers. The τ-value in the MK test reflects the correlation between the test sequence and time, and the p-value represents the significance level of the test results. The test results were divided into seven categories according to τand p-values to determine the trend change (Table 2).

. Changes in Vegetation Cover
The NDVI effectively monitors regional and global vegetation and is a major factor in the inversion of vegetation coverage. It is suitable for vegetation monitoring in the middle stage of vegetation development or in areas with medium vegetation cover, such as the YanHe Watershed. Several types of remote sensing techniques are utilized in such monitoring, among which MODIS NDVI has a strong temporal sequence and is available for free. Moreover, studies have shown that MODIS NDVI is more suitable for indicating vegetation conditions on the Loess Plateau; therefore, it was adopted in this study. Figure 3 shows in the inversion of vegetation coverage. It is suitable for vegetation monitoring in the middle stage of vegetation development or in areas with medium vegetation cover, such as the YanHe Watershed. Several types of remote sensing techniques are utilized in such monitoring, among which MODIS NDVI has a strong temporal sequence and is available for free. Moreover, studies have shown that MODIS NDVI is more suitable for indicating vegetation conditions on the Loess Plateau; therefore, it was adopted in this study.  The MK trend of NDVI represents the spatial distribution of vegetation change in the YanHe Watershed during 2000-2015 ( Figure 4). After commencing the GFGP, vegetation followed a significant growth trend in most parts of the YanHe Watershed, with extremely significant and significant increases accounting for 78.08% and 11.87% of the total area, respectively. Areas with noticeable vegetation growth were distributed primarily in the upstream and downstream areas of the YanHe River. The upper reaches of the river are characterized by complex topography, with steep slopes (>25°) accounting for more than 70% of the area [66], which has the most severe soil erosion in the whole The MK trend of NDVI represents the spatial distribution of vegetation change in the YanHe Watershed during 2000-2015 ( Figure 4). After commencing the GFGP, vegetation followed a significant growth trend in most parts of the YanHe Watershed, with extremely significant and significant increases accounting for 78.08% and 11.87% of the total area, respectively. Areas with noticeable vegetation growth were distributed primarily in the upstream and downstream areas of the YanHe River. The upper reaches of the river are characterized by complex topography, with steep slopes (>25 • ) accounting for more than 70% of the area [66], which has the most severe soil erosion in the whole basin and is a key area for the implementation of the GFGP. Thus, the vegetation cover in this area improved significantly. In the downstream area, agricultural activities were frequent; vegetation cover was poor at the start of the study; however, it was well restored after the implementation of the GFGP. In the middle reaches before the GFGP, the vegetation was relatively dense and had not changed significantly since 2000. The regions with significant vegetation degradation accounted for less than 1% of the total watershed area and were primarily distributed in the urban areas of Yan'an City. During urbanization, land use types are transformed from agricultural land and grassland to construction land, resulting in vegetation destruction and a declining NDVI trend. In general, after implementing the GFGP, the NDVI in the study area exhibited a steady increasing trend, indicating that the vegetation status of the YanHe Watershed, along with the ecological environment, had overall improved over time.
gions with significant vegetation degradation accounted for less than 1% of the total watershed area and were primarily distributed in the urban areas of Yan'an City. During urbanization, land use types are transformed from agricultural land and grassland to construction land, resulting in vegetation destruction and a declining NDVI trend. In general, after implementing the GFGP, the NDVI in the study area exhibited a steady increasing trend, indicating that the vegetation status of the YanHe Watershed, along with the ecological environment, had overall improved over time.    were slightly lower than those in the corresponding previous years. In general, the changing trend of ET in the YanHe Watershed was consistent with the overall change in vegetation, indicating that the GFGP affected regional vegetation growth and the hydrological cycle process and had a significant impact on terrestrial ET.

Changes in Evapotranspiration
The results in Figure 5b show the MK trend in the spatial distribution of ET in the YanHe Watershed from 2000 to 2015. ET in the watershed had increased significantly during the 15 years (2000-2015) (except for a few areas with significance <95%), with noticeable regional differences. Because of the fertile soil, adequate water sources, flat terrain, and frequent human activities, vegetation cover near the river was relatively stable since 2000, and ET increased slowly. As urbanization restricted the area around Yan'an New City, vegetation restoration was not apparent, and ET remained almost unchanged. The increase in ET in the northern part of the basin was slightly higher than that near the river. Agricultural land and grassland were widely distributed in this area, but an increase in ET was not evident due to limited water availability and heat conditions, and limited changes in vegetation coverage. Adequate water and heat conditions are suitable for vegetation growth; therefore, a medium increase in ET was widely distributed in this region, with scattered areas of significant ET increase. ISPRS Int. J. Geo-Inf. 2021, 10, x FOR PEER REVIEW 11 of 21 GFGP affected regional vegetation growth and the hydrological cycle process and had a significant impact on terrestrial ET. The results in Figure 5b show the MK trend in the spatial distribution of ET in the YanHe Watershed from 2000 to 2015. ET in the watershed had increased significantly during the 15 years (2000-2015) (except for a few areas with significance < 95%), with noticeable regional differences. Because of the fertile soil, adequate water sources, flat terrain, and frequent human activities, vegetation cover near the river was relatively stable since 2000, and ET increased slowly. As urbanization restricted the area around Yan'an New City, vegetation restoration was not apparent, and ET remained almost unchanged. The increase in ET in the northern part of the basin was slightly higher than that near the river. Agricultural land and grassland were widely distributed in this area, but an increase in ET was not evident due to limited water availability and heat conditions, and limited changes in vegetation coverage. Adequate water and heat conditions are suitable for vegetation growth; therefore, a medium increase in ET was widely distributed in this region, with scattered areas of significant ET increase.

Drought Assessment Based on ESI
The variation trends of different drought levels in the YanHe Watershed from 2000 to 2015 are shown in Figure 6a,b. The mild drought area declined significantly, reaching 107.4 km 2 •a −1 . From 2000 to 2007, the proportion of mild drought fluctuated, with no apparent increasing or decreasing trend during the early growth and development of the vegetation. The area of mild drought continued to decline during 2007-2014, rebounding slightly in 2015 due to the increasing impact of vegetation restoration on the regional water balance during the later stage of the GFGP. In contrast, the moderate and severe drought areas fluctuated upward, with rates of 38.58 km 2 •a −1 and 12.98 km 2 •a −1 , respectively. From 2000 to 2015, the YanHe Watershed experienced a significant decrease in mild drought and increases in no drought (normal level), moderate drought, and severe drought conditions. A significant regional difference was observed between dry and wet conditions, indicating a serious imbalance between water supply and demand, which should be focused upon by the government and relevant departments.

Drought Assessment Based on ESI
The variation trends of different drought levels in the YanHe Watershed from 2000 to 2015 are shown in Figure 6a,b. The mild drought area declined significantly, reaching 107.4 km 2 ·a −1 . From 2000 to 2007, the proportion of mild drought fluctuated, with no apparent increasing or decreasing trend during the early growth and development of the vegetation. The area of mild drought continued to decline during 2007-2014, rebounding slightly in 2015 due to the increasing impact of vegetation restoration on the regional water balance during the later stage of the GFGP. In contrast, the moderate and severe drought areas fluctuated upward, with rates of 38.58 km 2 ·a −1 and 12.98 km 2 ·a −1 , respectively. From 2000 to 2015, the YanHe Watershed experienced a significant decrease in mild drought and increases in no drought (normal level), moderate drought, and severe drought conditions. A significant regional difference was observed between dry and wet conditions, indicating a serious imbalance between water supply and demand, which should be focused upon by the government and relevant departments.
The spatial distribution map of the average ESI from 2000 to 2015 represents the degree of drought in the YanHe Watershed (Figure 6c). The arid region accounts for 63.86% of the total area, of which mild, moderate, and severe droughts accounted for 60.52%, 3.30%, and 0.04%, respectively. Comparing Figures 5b and 6c, the distribution pattern of the average ESI was found to be similar to the average ET over time, indicating that the degree and distribution of drought were significantly affected by ET. In Figure 6c, the positive ESI areas (green) are concentrated in the southern edge of the basin and scattered in the middle reaches, indicating a humid climate, and small impact of drought. The negative ESI area (orange-red) is widely distributed in the upper reaches of the YanHe River, indicating severe drought; it is also slightly distributed in the lower reaches of the basin, however, the degree of drought is low. This study is based on a pixel-by-pixel calculation of ESI based on remote-sensing MODIS products. By analyzing the spatial distribution of ESI results, the spatial variation characteristics of drought can be explored in detail, which help formulate the drought prevention measures. The spatial distribution map of the average ESI from 2000 to 2015 represents the degree of drought in the YanHe Watershed (Figure 6c). The arid region accounts for 63.86% of the total area, of which mild, moderate, and severe droughts accounted for 60.52%, 3.30%, and 0.04%, respectively. Comparing Figures 5b and 6c, the distribution pattern of the average ESI was found to be similar to the average ET over time, indicating that the degree and distribution of drought were significantly affected by ET. In Figure  6c, the positive ESI areas (green) are concentrated in the southern edge of the basin and scattered in the middle reaches, indicating a humid climate, and small impact of drought. The negative ESI area (orange-red) is widely distributed in the upper reaches of the YanHe River, indicating severe drought; it is also slightly distributed in the lower reaches of the basin, however, the degree of drought is low. This study is based on a pixel-by-pixel calculation of ESI based on remote-sensing MODIS products. By analyzing the spatial distribution of ESI results, the spatial variation characteristics of drought can be explored in detail, which help formulate the drought prevention measures. Figure 6d shows the MK trend in the spatial distribution of ESI in the YanHe Watershed from 2000 to 2015. Areas with significantly worse droughts were concentrated in the southwest of the basin. The main vegetation types in this area included forest and grassland with severe ET, which accelerated the loss of water resources leading to dry conditions. In summary, it was observed that vegetation restoration aggravates drought, particularly in relatively humid areas with high vegetation coverage, which may negatively affect the balance of surface water, groundwater, and soil water storage. A con-  Areas with significantly worse droughts were concentrated in the southwest of the basin. The main vegetation types in this area included forest and grassland with severe ET, which accelerated the loss of water resources leading to dry conditions. In summary, it was observed that vegetation restoration aggravates drought, particularly in relatively humid areas with high vegetation coverage, which may negatively affect the balance of surface water, groundwater, and soil water storage. A continuation of this phenomenon could adversely affect the quality of forestry ecological projects and the impact of the GFGP.

Evaluation and Calibration of the Modified SWAT
In this study, soil conservation in the YanHe Watershed was simulated using the modified SWAT model (Figure 7). The validation process involved input of the parameter values determined during the calibration process into the model. The simulation results for the YanHe Watershed improved during the calibration (R 2 and NSE increased from 0.63 to 0.76 and 0.46 to 0.51, respectively) and verification periods (R 2 and NSE increased from 0.82 to 0.93 and 0.57 to 0.61, respectively). After introducing the remote-sensing data, the SWAT sediment simulation increased the R 2 , NSE, and the model simulation accuracy. During the calibration and verification periods, the monthly simulated value of sediment was consistent with the observed value and the modified model simulation result was closer to the observed value with a higher degree of correlation than that of the original model. In the study period, R 2 > 0.5 and NSE > 0.5 were achieved in the simulation results. Therefore, the modified SWAT model would be suitable for simulating sediment transport in the YanHe Watershed. from 0.63 to 0.76 and 0.46 to 0.51, respectively) and verification periods (R 2 and NSE increased from 0.82 to 0.93 and 0.57 to 0.61, respectively). After introducing the remote-sensing data, the SWAT sediment simulation increased the R 2 , NSE, and the model simulation accuracy. During the calibration and verification periods, the monthly simulated value of sediment was consistent with the observed value and the modified model simulation result was closer to the observed value with a higher degree of correlation than that of the original model. In the study period, R 2 > 0.5 and NSE > 0.5 were achieved in the simulation results. Therefore, the modified SWAT model would be suitable for simulating sediment transport in the YanHe Watershed.

Coupling Analysis of Soil Conservation Services and Drought
Since 2000, soil conservation services and drought conditions in the YanHe Watershed have been constantly changing under the influence of vegetation restoration. Based on the monthly data of ESI, NDVI, precipitation (PCP), relative humidity (RH), wind speed (WIND), and temperature (TEM), the GRA method was adopted to evaluate the correlation between soil conservation services, drought, and other factors. The analysis results are presented in Table 3. The correlation sequence of soil conservation services and related factors (was RESI > RNVDI > RPCP > RRH > RWIND > RTEM, which showed that,  Figure 8b shows the distribution pattern of SCM levels in the YanHe Watershed from 2000 to 2015, with significant spatial variations. The lower SCM accounted for 55.08% of the total area, most of which was located on both banks of the river, and the rest scattered on agricultural land. The soil erodibility factor (K) around the YanHe River was large; the erosion was severe and not conducive to the implementation of engineering measures. The surface soil activity of agricultural land was frequent because of human activities. In spring and winter, agricultural land was mostly without vegetation cover, with limited soil conservation functions. The low SCM accounted for 29.67% of the total area, which was concentrated at the bottom of the slope, with a low slope and undulating topography. The primary land use type was grassland, and the topography and surface cover of this area were relatively stable. The medium and high SCM had strong soil conservation functions, but their proportions were relatively small (10.47% and 4.78%, respectively). Most of them were located in the middle and lower reaches of the YanHe River because of the flat terrain, thick soil layer, and adequate water and heat conditions. The soil was not easily eroded, and acceptable vegetation and maintenance measures had been implemented. In contrast, soil erosion in the upper reaches was severe and not conducive to vegetation restoration because of the deep valleys, broken ground surface, and steep slopes (>25 • ). A high SCM was distributed in certain individual years. To reflect the average level of the soil retention function in the entire watershed, the high SCM was smoothed when calculating the multi-year average SCM. Figure 8c shows the varying trend of the SCM in the YanHe Watershed. To reduce the impact of extreme climatic conditions, the variations in SCM during the late stage of the GFGP (2011-2015) over those during the early stage (2001)(2002)(2003)(2004)(2005), were used to analyze the spatial trend in SCM. Soil conservation services in most areas of the YanHe River Basin have been increasing. The middle and lower reaches of the YanHe River are suitable for vegetation growth and ecological engineering development, which are important supply areas for soil conservation services in the basin. The upstream was found to be the key area of successful GFGP implementation, where soil erosion has been curbed and soil conservation functions have been enhanced. However, an increase in the supply of soil conservation services near the river was not apparent.

Coupling Analysis of Soil Conservation Services and Drought
Since 2000, soil conservation services and drought conditions in the YanHe Watershed have been constantly changing under the influence of vegetation restoration. Based on the monthly data of ESI, NDVI, precipitation (PCP), relative humidity (RH), wind speed (WIND), and temperature (TEM), the GRA method was adopted to evaluate the correlation between soil conservation services, drought, and other factors. The analysis results are presented in Table 3. The correlation sequence of soil conservation services and related factors (was R ESI > R NVDI > R PCP > R RH > R WIND > R TEM , which showed that, among many factors, soil conservation services had the highest correlation with ESI, followed by NDVI and precipitation, and the least correlation with temperature. Based on the GRA, the correlation between the SCM and ESI in the YanHe Watershed was analyzed on a spatial scale; the distribution of the spatial correlation coefficients between them are shown in Figure 9. The areas with negatively correlated SCM and ESI (red) are primarily distributed in the middle and lower reaches of the YanHe River, and the enhancement of soil conservation services in these regions was accompanied by the intensification of drought. Areas with a highly negative correlation (<−0.5) were scattered in the regions with high vegetation coverage. Vegetation can reduce the effect of raindrops, thereby reducing rainwater splash erosion; reduce the peak of floods; delay surface runoff; and reduce surface erosion by water; which effectively curbs watershed soil erosion and enhances soil conservation functions. However, the plant canopy intercepts rainfall and reduces the amount of water that reaches the ground. Vegetation restoration increases transpiration and biomass, resulting in increased surface crust and reduced infiltration. These factors lead to a decrease in groundwater recharge, an increase in water consumption from deep groundwater, a decrease in soil water content, and an aggravation of regional drought. The areas where SCM and ESI are positively correlated (green) are scattered in the upper reaches of the YanHe River, occupying a small proportion. Drought in these areas has been alleviated by the enhancement of soil conservation services. Because of the differences in vegetation types, topography, climate, and other factors, the correlation between soil conservation services and drought exhibited spatial variations.

Discussion
This study focused on drought and soil conservation services to evaluate the ecological effects of vegetation restoration in the YanHe Watershed. Using remote sensing data, this study analyzed spatiotemporal variations in vegetation cover and ET in the watershed. To meet the requirements of the study period and analyze data on a spatial scale, this study utilized currently available MODIS products. The results of this study are consistent with other studies: The GRGP significantly increased vegetation cover in the YanHe River Basin, and the large-scale restoration of vegetation [67,68] has been accompanied by a significant increase in regional ET and water consumption [69]. This study used ESI to characterize the degree of drought in the basin and evaluate the response of drought to the increase in ET. Compared to SPI and SPEI, the ESI accounted for the role of ET in the evolution of drought, which can better reflect the response of vegetation and soil moisture to drought [70], although it did not consider the impact of meteorological factors (precipitation) on drought. Herein, soil conservation services in the YanHe Watershed were simulated, based on the modified SWAT model which was integrated with remote sensing data. The original SWAT model calculates the C-factor based on land use data, which does not accurately express the variation in vegetation cover. The study recalculated CUSLE using a remotely sensed NDVI value to better reflect real vegetation cover changes when simulating soil conservation. The limitations of this study include: (1) Because of limited availability of remote sensing data sources, the modified SWAT did not integrate NDVI data with high temporal and spatial resolutions;

Discussion
This study focused on drought and soil conservation services to evaluate the ecological effects of vegetation restoration in the YanHe Watershed. Using remote sensing data, this study analyzed spatiotemporal variations in vegetation cover and ET in the watershed. To meet the requirements of the study period and analyze data on a spatial scale, this study utilized currently available MODIS products. The results of this study are consistent with other studies: The GRGP significantly increased vegetation cover in the YanHe River Basin, and the large-scale restoration of vegetation [67,68] has been accompanied by a significant increase in regional ET and water consumption [69]. This study used ESI to characterize the degree of drought in the basin and evaluate the response of drought to the increase in ET. Compared to SPI and SPEI, the ESI accounted for the role of ET in the evolution of drought, which can better reflect the response of vegetation and soil moisture to drought [70], although it did not consider the impact of meteorological factors (precipitation) on drought. Herein, soil conservation services in the YanHe Watershed were simulated, based on the modified SWAT model which was integrated with remote sensing data. The original SWAT model calculates the C-factor based on land use data, which does not accurately express the variation in vegetation cover. The study recalculated C USLE using a remotely sensed NDVI value to better reflect real vegetation cover changes when simulating soil conservation. The limitations of this study include: (1) Because of limited availability of remote sensing data sources, the modified SWAT did not integrate NDVI data with high temporal and spatial resolutions; however, the simulation accuracy was higher than that of the original model; (2) there is a lack of research on value accounting of ESs, which plays a limited role in the ecological compensation field of the basin. Further research should be conducted to consider the impact of climate change and human activities, analyze the relationship between drought and ESs from the perspective of physical mechanisms, and evaluate the value of ESs to guide eco-compensation.

Conclusions
The main findings of the study are as follows: 1.
GRA has significant advantages in studying the factors affecting ESs. Soil conservation services in the YanHe Watershed are closely related to ESI and NDVI, followed by rainfall. Spatial variations exist in the correlations between soil conservation services and drought. In the middle and lower reaches of the YanHe River, drought in the regions covered by vegetation was aggravated with the enhancement of soil conservation services. The drought in parts of the mountainous area in the upper reaches of the YanHe River decreased with the enhancement of soil conservation services.

2.
From 2000 to 2015, with the implementation of the GFGP, vegetation status significantly improved, but plant water consumption increased with increasing ET. The average annual growth rate of NDVI in the entire basin was 0.0245 × 10a −1 , whereas the growth rate of ET was 14.35 mm·a −1 .

3.
Arid regions in the YanHe Watershed accounted for 63.86% of the total area, and the spatial distribution of drought increased from south to north. From 2000 to 2015, the overall distribution pattern of wet and dry conditions in the basin varied negligibly. Drought in parts of the midstream area of the YanHe River has been alleviated, while the humid regions have become drier.

5.
The integration of remote sensing NDVI into SWAT enabled the fusion of remote sensing and non-remote sensing data, which enhanced the reliability of the original model and improved simulation accuracy.
Considering the water conservation and supply of soil conservation services, vegetation restoration in the YanHe Watershed should incorporate low-water-consumption tree species and restrict density of afforestation. The study results will help in evaluating the ecological effects of the GFGP and provide valuable information for formulating ecological policies on the Loess Plateau thereby enabling sustainable development of the ecosystem.