How Do Trade-O ﬀ s and Synergies between Ecosystem Services Change in the Long Period? The Case Study of Uxin, Inner Mongolia, China

: Ecosystem services management should often expect to deal with non-linearities due to trade-o ﬀ s and synergies between ecosystem services (ES). Therefore, it is important to analyze long-term trends in ES development and utilization to understand their responses to climate change and intensiﬁcation of human activities. In this paper, the region of Uxin in Inner Mongolia, China, was chosen as a case study area to describe the spatial distribution and trends of 5 ES indicators. Changes in relationships between ES and driving forces of dynamics of ES relationships were analyzed for the period 1979–2016 using a stepwise regression. We found that: the magnitude and directions in ES relationships changed during this extended period; those changes are inﬂuenced by climate factors, land use change, technological progress, and population demonstrate that trade-o ﬀ s and synergies among ES are spatially heterogeneous.


Introduction
Ecosystem services (ES) refer to the conditions and utilities of the natural environment, as well as to the benefits people obtain from these ecosystems [1]. ES not only maintain the life cycle process, ensure biodiversity, purify the environment, and perform other functions and ecological processes humans rely on, but they also provide raw materials, food, fiber, clean water and recreation, which are necessary for human well-being and production [2,3]. Over the past decades, the average life span of humankind significantly improved, and poverty alleviated in some areas. At the same time, adequate knowledge of ES and understanding of scientific principles of effective management of ES are still lacking. For example, the excessive pursuit of provisioning services, approximately 60% (15 out of 24) of the ES are being degraded in the Millennium Ecosystem Assessment [4]. The degradation and loss of ES have a negative impact on human well-being, and directly threaten regional, national, and global ecological security [5].
The nonlinear relationship between different ES and complex patterns of their utilization by humans are common [6], which often makes ES dependent on each other and interact in multiple ways, displaying both trade-offs and synergies [7]. Trade-offs, or negative relationships, refer to an increase in the supply of one ES that leads to the reduction in the supply of another ES [8]. On the contrary, synergies, or positive relationships, mean that the supply of multiple ES increases or decreases simultaneously [5]. Due to the existence of such relationships, unexpected outcomes of ecosystem management practice are common. For example, the increase of provisioning services (grain, wood, etc.) may lead to the reduction of regulating services (nutrient cycling, soil conservation, etc.) [9]. As another example, an increase in soil retention would improve soil carbon sequestration capacity [10]. Some European countries tried to encourage farmers to adopt more environmentally friendly cultivation techniques by offering financial incentives. However, these measures had no effect on the conservation of endangered species [11]. Therefore, consideration of trade-offs and synergies between ES is crucial for landscape planning and land management. It will help avoid costly adverse effects, and promote multi-functionality [12].
Influenced by natural and anthropogenic drivers, relationships between ES change spatially and temporally. Studies have analyzed the common drivers that underlie changes in ES, including climate change, land use changes, as well as policy interventions [13][14][15]. For example, the transformation of cultivated land into shrub and grassland enhanced soil conservation (SC) and carbon sequestration capacity but reduced water yield (WY) in the Loess Plateau, China [16]. In northern Shaanxi, SC and net primary productivity (NPP) had variable synergistic relationships in different land use types [10]. Finally, in the Yanhe river, SC, water retention (WR), and WY changed at sub-basin scale after the implementation of the Grain to Green Program (GTGP) [17]. Moreover, ES relationships occur at a different spatial and temporal scale [8], which increases uncertainties to be managed [18]. Previous studies have documented trade-offs and synergies among multiple ES [19][20][21][22][23]. However, they were mainly focused on single time periods and did not consider temporal changes [24]. Few studies have recently focused on changes in those relationships at different temporal scales. For example, simulations of ES relationships for the period of 2001−2070 revealed changes in some of those relationships involving regulating services [24]. A case study of Switzerland focused on a ten-year period assessment of ES and revealed robustness of the determined relationships between regulating and cultural services [25]. Although quantitative testing of ES and their relationships across temporal scales have been conducted recently, more evidence is needed to understand their dynamics and elucidate mechanisms underlying changes in ES relationships in a longer term to sustainably manage multiple ES.
In this study, we attempted to: (1) describe the spatial distribution and spatial trends of 5 ES indicators (i.e., livestock breeding (SHEEP); grain production (GRAIN); NPP; sandstorm prevention (SP); and WR) from 1979 to 2016 in the region of Uxin; (2) quantify changes in ES relationships over the same period; and to (3) reveal the driving factors behind change in the ES and their relationships.

The Study Area
Uxin is located in the southeast part of the Ordos Plateau in Inner Mongolia, North China. It has an average altitude of 1300 m, and extends in latitude from 37 • 38 N to 39 • 23 N and in longitude from 108 • 17 E to 109 • 40 E. Uxin has a typical temperate continental climate, with a mean annual precipitation of about 350 mm, a mean annual evaporation of 2200 mm, and a mean annual temperature of 6.8 • C. Fixed and moving sand dunes cover the majority of its landscape. Aeolian sandy soils and kastanozems are the most common soil types. Shrubs and subshrubs are the dominant vegetation type (e.g., Caragana intermedia and Artemisia ordosica). As a typical agro-pastoral transitional zone of northern China, the main land use is livestock husbandry and farming ( Figure 1).

Quantification of Ecosystem Services
We investigated 5 ES indicators, including livestock breeding (SHEEP), grain production (GRAIN), net primary productivity (NPP), sandstorm prevention (SP), and water retention (WR). The selection criteria included the following. ES should: (1) be relevant to the well-being of local households and surrounding ecosystem; (2) has available remotely sensed data for spatially explicit monitoring of ES; and (3) require data that were available for long-term temporal ES relationship assessment. SHEEP and GRAIN are important to human well-being of local households. NPP reflects the growth status of vegetation. SP helps to reduce sandstorm disasters for Uxin and downwind areas (e.g. Xi'an city, Shaanxi Province, etc.). WR plays an important role in ecosystem vegetation recovery and household livelihoods. We mapped the annual value layer for each indicator from 1979 to 2016, which was then used to describe changes in ES. All data were resampled to a spatial resolution of 250 m (The results of 5 ES are available in Supplementary Materials). Table 1 shows the data used to assess the 5 ES indicators. Climate data, such as mean temperature, precipitation, and wind speed, were interpolated using ArcGIS 10.3. Based on data availability, we used land use and land cover (LULC) maps of the years 1978,1987,1992,1997,2002,2007,2012,2014, and 2016

Quantification of Ecosystem Services
We investigated 5 ES indicators, including livestock breeding (SHEEP), grain production (GRAIN), net primary productivity (NPP), sandstorm prevention (SP), and water retention (WR). The selection criteria included the following. ES should: (1) be relevant to the well-being of local households and surrounding ecosystem; (2) has available remotely sensed data for spatially explicit monitoring of ES; and (3) require data that were available for long-term temporal ES relationship assessment. SHEEP and GRAIN are important to human well-being of local households. NPP reflects the growth status of vegetation. SP helps to reduce sandstorm disasters for Uxin and downwind areas (e.g., Xi'an city, Shaanxi Province, etc.). WR plays an important role in ecosystem vegetation recovery and household livelihoods. We mapped the annual value layer for each indicator from 1979 to 2016, which was then used to describe changes in ES. All data were resampled to a spatial resolution of 250 m (The results of 5 ES are available in Supplementary Materials). Table 1 shows the data used to assess the 5 ES indicators. Climate data, such as mean temperature, precipitation, and wind speed, were interpolated using ArcGIS 10.3. Based on data availability, we used land use and land cover (LULC) maps of the years 1978, 1987, 1992, 1997, 2002, 2007, 2012, 2014, and 2016 to represent the conditions for the following intervals: 1979-1987, 1987-1992, 1992-1997, 1997-2002, 2002-2007, 2007-2012, 2012-2014, and 2014-2016.  [26]. In order to convert the aggregated, or non-spatial, statistical data into a spatial layer, the normalized difference vegetation index (NDVI) from remotely sensed satellite data was selected as a proxy to characterize the capacity of an ecosystem to provide forage. The density of sheep units is dependent on accumulated NDVI values as follows:

Data Sources
where SU i , a vector map with map unit attributes linked to livestock numbers in 6 counties of Uxin. NDVI i , sum was aggregated within the boundaries of each county for the year i. NDVI i, max is the maximum value component (MVC) for the year i. SU i, spati is the spatial layer of the livestock number for the year i. The method of preparing spatial data for the GRAIN was similar to the one used for the livestock number. In this paper, farmland is extracted according to the NDVI threshold of 0.34. With the ArcGIS Zonal Statistics tool [27], this statistic mean value is calculated for each cropland zone defined by a zone dataset, based on values from NDVI dataset (a value raster). This method is consistent with the method of spatial representation of livestock number.

Net Primary Productivity (NPP)
The Carnegie-Ames-Stanford Approach (CASA) model was used to calculate the NPP [28]. This model takes into account the effects of solar radiation, temperature, and water stress on NPP. where APAR(x, t) is the vegetation photosynthetic effective radiation (g C·m −2 ·month −1 ) absorbed by grid x at month t; and ε(x, t) is the actual utilization rate of light energy of vegetation (g C). APAR(x, t) was calculated using the following equation: where SOL(x, t) is the total amount of solar radiation (MJ·month −1 ) at grid x for a specific month t [29]; 0.5 refers to the proportion of the solar effective radiation (the wavelength is 0.4-0.7 µm) to the total solar radiation; FPAR(x, t) is the absorption ratio of vegetation photosynthetic effective radiation. The utilization rate of light energy by vegetation is mainly affected by temperature and water conditions and calculated as follows: where ε max represents the maximum light energy utilization (g C·MJ −1 ); T ε1 (x, t) and T ε2 (x, t) indicate the stress effects of the lowest and the highest temperatures on the utilization of light energy by vegetation; and W ε (x, t) is the influence coefficient of water stress, which was set to 0.542 gC·MJ −1 for the case study area [30].

Sandstorm Prevention (SP)
Based on the Revised Wind Erosion Equation (RWEQ) [31], sandstorm prevention (SP) can be estimated as the difference between potential wind erosion (S Lp ) and actual wind erosion (S L ). The formulas are as follows: where S L (t·km −2 ·a −1 ) is the actual soil loss caused by wind erosion; S Lp (t·km −2 ·a −1 ) is the potential soil loss; QMAX (kg·m −1 ) is the maximum transfer capacity; S (m) is the critical field length, defined as the distance at 63% of QMAX; z (m) is the maximum wind erosion distance, set to 50 m for the study area; WF (kg·m −1 ) is the climate factor, influenced by soil moisture, wind speed, and snow cover; EF is the soil erodibility factor; SCF represents soil crusting; C is the vegetation cover factor; and K' is the surface roughness factor. For the detailed formulas of RWEQ, see Jiang, et al. [32]. Considering factors such as topography, soil thickness, and permeability, WR was estimated based on the principle of water balance at watershed scale [33] as follows: where WR is the average water retention (mm); Velocity is the flow coefficient [29]; TI is the topographic index; K sat is the soil saturated hydraulic conductivity (cm·d −1 ) [34]; and WY is the water yield of the basin (mm), calculated by InVEST 3.1.0 [33]. It is calculated as follows: where Yield xj (mm) is the annual water yield for LULC type j in a given grid x; P x (mm) is the annual precipitation of grid cell x; and AET xj is the annual actual evapotranspiration (mm) of the LULC type j in the grid x. The AET xj was calculated using the following equations: where ω x is the parameter of natural climate and soil properties, defined as the ratio of annual water demand for plants and precipitation; PET xj is the annual potential evapotranspiration of LULC type j in grid x [35]; Z is the seasonal parameter for seasonal rainfall distribution, which was set to 11.54 for the study area [36]; AWC x (mm) is the available water for plants in grid x [37]. In addition, the InVEST water yield model also requires tabulated biophysical parameters including the maximum root depth of the vegetation of each LULC type and the evapotranspiration coefficient (K c ). The details are shown in Table 2.

Trend Analysis of ES and Their Relationships
An image trend analysis was conducted to analyze ES changes in the period 1979-2016 in Uxin using the following equation: where, y k represents a specific ES layer in year k (k = 1979, 1980 . . . , 2016, n = 38 years). x k is a grid layer; each grid cell is set to the same value for the year k; b is the trend slope layer; A grid cell with b < 0 indicates a decreasing trend in ES in the period analyzed, while a grid cell with b > 0 indicates an increasing trend. The Arcpy scripts for trend analysis is available in Supplementary Materials. ES relationships can be modeled by the Pearson correlation coefficient [18]. For some services, a positive correlation indicates a synergy, while a negative correlation represents a trade-off. It is calculated as follows: where, r is correlation coefficient layer of an ES pair. In spatial distribution characteristics of trade-offs analysis, x k and y k represent 5 ES layers in year k (k = 1979, 1980 . . . , 2016, n = 38 years). In temporal dynamics of ES trade-offs analysis, x k and y k represent k th pixel value of each ES pair in a specific year (e.g., 1979, 1980 . . . , 2016). The significance of b value and r were assessed by t-test using 2-tail t-distribution at 95% confidence interval.

Driving Forces of Ecosystem Services Interactions
Climate change, humans activities, and technological progress are main drivers of quantity and quality changes in ES [38]. We hypothesize 5 ES indicators are driven by climate factors and land use and indirectly influenced by technical progress and population growth, which might result in spatiotemporal changes in relationships of 5 ES. Land cover transformations driven by economy and population growth on one hand and technological progress, on the other, change the demand for natural resources and affect consumption patterns of different ES. To test this hypothesis, 7 independent variables were selected to perform a stepwise regression (Table 3), which minimizes the problem of multi-collinearity in variables. A max-min standardization was performed for all variables; the statistical analysis was performed using R 3.5.0 [39]. The full reproducible R code and data are available in Supplementary Materials.  ranged from 0.38 to 0.60 (with an area ratio of 1.37%). The slope of the eastern and northern regions was between 0.03 and 0.38 (with an area ratio of 9.07%).

GRAIN
During the 38-year period GRAIN also increased 7.54 times from 1.50 × 10 4 t to 1.28 × 10 5 t, with an average annual change of 2.98 × 10 3 t (see Figure 2-B1). The highest value of GRAIN is found in about 10% of the total area, mainly in the southeastern, eastern, and northern parts of Uxin. Areas with lowest GRAIN value are scattered across the Uxin region and occupy about 4% of the total area (see Figure 2-B2). GRAIN exhibited an increasing trend (see Figure 2-B3), especially in the southern region, where its slope ranged from 0.38 to 0.60 (with an area ratio of 1.37%). The slope of the eastern and northern regions was between 0.03 and 0.38 (with an area ratio of 9.07%).

Net Primary Productivity
Two trends of NPP can be identified over the past 38 years (see  (Figure 2-D2). SP showed a consistent increasing trend its distribution characteristics were similar to its spatial distribution (see Figure 2-D3).

Water Retention
Different from the study of similar region in Xilin River Basin [36], WR showed an increasing trend in the period 1979-2016, with its rate changing from 3.75 × 10 8 m 3 in 1979 to 7.65 × 10 8 m 3 in 2016 resulting in an annual change of 1.03 × 10 7 m 3 (see Figure 2-E1). The mean WR was greater than 15 mm in more than 90% of the total region area, mainly in the east and south of Uxin (see Figure 2-E2). About 24.16% of the total region area had mean WR in the range of 15-30mm; 49.54% had mean WR of 30-50 mm; and 16.82% had mean WR of 50-100 mm. Over 16.33% of Uxin, mainly the northeastern area, showed a slight increasing trend for WR. Only 1.88% of the total region, areas scattered the south and central part, showed a decreasing trend for WR (see Figure 2-E3). Water is the main limiting factor of vegetation growth and agricultural production in this area. With WR increased gradually, there will be enough water resources to meet the demand of agricultural expansion and population growth, and indirectly reducing vegetation degradation.

Temporal Dynamics of ES Relationships
Trends in ES relationships varied in the 38 years analyzed (Figure 3, the bottom-right triangle diagram matrix). Positive correlations between SHEEP-GRAIN, SHEEP-NPP, and GRAIN-NPP have been gradually strengthening, while negative correlations among some ES pairs, such as SHEEP-SP, GRAIN-SP, and NPP-SP, displayed a weakening trend during the study period. Yet some positive

Spatial Patterns of ES Relationships
According to the spatial patterns of relationships between various ES ( Figure 3, the top-left triangle diagram matrices), the pairs GRAIN-NPP, NPP-WR and SP-WR exhibited almost uniform positive correlation across the whole region, while SHEEP-GRAIN, SHEEP-NPP, SHEEP-SP and NPP-SP exhibited a spatial clustering in positive and negative correlations. Finally, no specific patterns are observable for other pairs of ES.

Driving Forces of ES Relationship Changes
We used the stepwise regression to analyze the driving forces of change in the 10 relationships among ES pairs (see Table 4). In the period 1979-2016, SHEEP-GRAIN showed a positive relationship with GRS_PRS and AGMACH, suggesting that both the increased grazing pressure and the investments in the mechanical power of agriculture and animal husbandry promoted an improvement of the synergies with SHEEP-GRAIN. SHEEP-NPP also showed a positive relationship with GRS_PRS and AGMACH, while the increase of POP tended to inhibit the increase of these two synergy relationships. SHEEP-SP was negatively correlated with WIN and AGMACH. SHEEP-WR had a negative correlation with AGMACH, indicating that the increase in investments in artificial grassland by herdsmen were accompanied by an increase in sheep number and by an excessive consumption of groundwater, which tended to decrease their synergy relationships. GRAIN-NPP had a positive correlation with TEM, whereby increases in cumulative temperature tended to increase the synergies with GRAIN-NPP. GRAIN-SP was positively correlated with TEM, while GRAIN-WR was negatively correlated with changes in AGMACH. NPP-SP had a significantly positive correlation with TEM. NPP-WR was negatively correlated with GRS_PRS. SP-WR was positively correlated with WIN, GRS_PRS, and POP, and was negatively correlated with AGMACH.

Temporal Dynamics of ES Relationships
The temporal dynamics of the trade-offs and synergies among ES reflected the responses to social and environmental changes [24]. Furthermore, this study also found that both trade-offs and synergies were not invariable in the long period. On the contrary, they were likely influenced by a variety of drivers, such as climate, human activities, and technological progress. According to the results of the stepwise regression (see Table 4), there are four possible types of relationship between driving forces and ES pairs (see Figure 4). The direct and indirect driving forces of the relationship between ecosystem service pairs are shown in Table 5. The four types of relationships between driving forces and ES pairs are briefly described below.
(i) Driving factors directly affect two ES, increasing or decreasing the provision of ES, and resulting in changes in trade-offs and synergies (see Figure 4a). In this study we showed the presence of a synergistic relationship between SHEEP and GRAIN (see Figure 2, bottom-right part), which followed an increasing trend over the 38 years analyzed. The stepwise regression analysis showed that SHEEP and GRAIN were positively affected by technical progress factors. The mechanical input in crop farming and cultivated land during the last three decades led to an increase in food production by farm families. Triggered by the prohibition of open grazing policy since the early 2000s [40], herdsmen families increased their investments in irrigation for artificial pasture and fenced grazing, resulting in an increase in sheep breeding. Water is the main limiting factor of vegetation growth and agricultural production in this area. This study showed that the synergies between GRAIN-WR were transformed into trade-offs in the 38 years analyzed (see Figure 3), and were negatively affected by technical progress factors (see Table 4). By the late 1990s, when the mechanical input in crop farming increased, the water consumption by intensive agriculture increased gradually, while the WR in the region did not increase significantly (see Figure 2). Therefore, the synergy between GRAIN-WR gradually changed into a trade-off. Like GRAIN-WR, SHEEP-WR was negatively correlated with AGMACH, suggesting that the mechanical input in artificial grassland directly increased SHEEP provision and an excessive consumption of groundwater, tending to decrease their synergy relationship. input in crop farming increased, the water consumption by intensive agriculture increased gradually, while the WR in the region did not increase significantly (see Figure 2). Therefore, the synergy between GRAIN-WR gradually changed into a trade-off. Like GRAIN-WR, SHEEP-WR was negatively correlated with AGMACH, suggesting that the mechanical input in artificial grassland directly increased SHEEP provision and an excessive consumption of groundwater, tending to decrease their synergy relationship.     (ii) Driving factors directly affect one ES and indirectly affect another ES, resulting in changes in trade-offs and synergies (see Figure 4b). SHEEP-NPP showed a negative relationship with POP. Over the 38 years investigated, population growth led to a direct increase in sheep supply and to an increase in water consumption, thereby indirectly inhibiting the increase of NPP. SHEEP-SP had a negative correlation with AGMACH (see Table 4), indicating that increased investments in artificial grassland by herdsmen were accompanied by an increase in sheep number and by an excessive consumption of groundwater, which prevented vegetation recovery, thereby indirectly decreasing SP. GRAIN-SP was positively correlated with TEM. Due to the extension of growing seasons thanks to climate warming, which increased both grain yield [41] and biomass (NPP) [42], vegetation restoration was promoted, which indirectly increased SP. Similarly, thanks to climate warming (TEM), NPP gradually increased, which in turn led to changes in the supply of SP, resulting in the increase of the synergy relationship between NPP and SP during the period investigated. The increase in AGMACH resulted in an increase in the trade-offs between SP and WR. The increase of mechanical input in crop farming directly decreased the supply of WR, and had a negative impact on vegetation restoration in the surrounding ecosystems, which in turn indirectly reduced SP provision. SP-WR was positively correlated with POP. Population growth directly led to an increase in water use (i.e., a decrease in WR), and reduced the trend of the adverse impact of vegetation restoration on SP. SP and WR both decreased at the same time, thus increasing the synergy between SP and WR.
(iii) Driving factors directly affect one ES, while having at the same time a smaller effect on another ES. Changes in the supply capacity of only one ES, will eventually lead to changes in trade-offs and synergies (see Figure 4c). SHEEP-SP showed a positive relationship with WIN. The decrease in wind speed led to an increase in SP, which led to an increase in the synergistic relationship between them. FAM_LUI had a positive impact on GRAIN-NPP. Agricultural intensification promoted an increase in grain yield, thus increasing the synergistic relationship between them. The increase in grassland utilization (GRS_PRS) resulted in a decrease in NPP. Both WR, and the trade-offs between NPP and WR, recorded an increase in the 39 years investigated. SP-WR was positively correlated with WIN, indicating that the decrease of the average wind speed led to an increase in SP. Moreover, the increase in WR fostered an increase in the synergistic relationship between them.
(iv) Driving factors indirectly affect two ES, changing their supply rate and influencing their trade-offs and synergies (see Figure 4d). SP-WR was positively correlated with GRS_PRS, suggesting that the increase of grazing pressure on grassland directly reduced vegetation cover, thereby indirectly reducing SP and WR, and increasing the synergistic relationship between them.

Trade-offs and Synergies, and ES Spatial Heterogeneity
Our results demonstrate that trade-offs and synergies among ES are spatially heterogeneous. Li et al. (2017) revealed that in grasslands, the relationships between ES tend to be in the form of synergies, while in built-up land and farmland they tend to be in the form of trade-offs [43]. In our research, SHEEP-GRAIN, SHEEP-NPP, SHEEP-SP, and NPP-SP exhibited a spatial clustering in positive and negative correlations. Land use conflicts were one of the causes of trade-offs and synergies between ES [12,44,45]. The heterogeneity of ES is usually caused by land use [46], whereby different types of land use provide different levels of ES [10]. For example, the spatially heterogeneous relationship between SHEEP and NPP, caused by their regional distribution (see Figure 2-A2 and C2), resulted in an imbalance in ES trade-offs. The area with the lowest SHEEP value is located in the northwestern part of Uxin, while the areas with a lower vegetation coverage showed a decrease trend (see Figure 2-A3). The areas with high SHEEP values are located in the northeastern and southern regions, and showed an increasing trend over the 38 years analyzed (see Figure 2-A3). NPP exhibited similar spatial patterns (see Figure 2-C2). However, it showed trends of change that are opposite to those of SHEEP (see Figure 2-C3). Xu et al. (2017) also found that, due to the presence of hot-spot areas of ES in Ningxia, which contains 37% of the cultivated land and produces 57% of grain, the provision and regulation services in the region did not appear to develop trade-off relationships [44].

Management Implications
Our research suggests that it is necessary to perform a long-term analysis of ES and their drivers, to better understand the trade-offs and synergies among them. Due to the different sensibility of ES to climate change and human disturbance, changes in ES relationships can be abrupt with high inter-annual precipitation variations [24]; hence, the ES were chosen at different time periods, and this may have influenced the results of the various ES interactions (see Figure 5). Hence, combined with the statistical test, the timely assessment and monitoring of ES and their relationships are needed to robustly identify the change trends of each ES, to detect threshold or lag effects in ES interactions [43], to help avoid surprising trade-offs, and to take advantage of emerging synergies [24].
Our research suggests that it is necessary to perform a long-term analysis of ES and their drivers, to better understand the trade-offs and synergies among them. Due to the different sensibility of ES to climate change and human disturbance, changes in ES relationships can be abrupt with high interannual precipitation variations [24]; hence, the ES were chosen at different time periods, and this may have influenced the results of the various ES interactions (see Figure 5). Hence, combined with the statistical test, the timely assessment and monitoring of ES and their relationships are needed to robustly identify the change trends of each ES, to detect threshold or lag effects in ES interactions [43], to help avoid surprising trade-offs, and to take advantage of emerging synergies [24].  [43]). ES1 and ES2 have a synergy relationship in the period t0-t1, a trade-off relationship in the period t1-t2 and t2-t3, and a synergy relationship after t3.
According our study, the magnitude and direction of ES trade-offs and synergies may also vary, responding to changes in both natural and anthropogenic drivers. Our results suggest that more adaptive approaches for ecosystem management are required [47]. For example, revegetation in arid areas creates potentially conflicting demands for water between the ecosystem and humans [48]; Figure 5. Changes in ES and in their interactions in terms of trade-offs or synergies (modified from [43]). ES 1 and ES 2 have a synergy relationship in the period t 0 -t 1 , a trade-off relationship in the period t 1 -t 2 and t 2 -t 3 , and a synergy relationship after t 3 .
According our study, the magnitude and direction of ES trade-offs and synergies may also vary, responding to changes in both natural and anthropogenic drivers. Our results suggest that more adaptive approaches for ecosystem management are required [47]. For example, revegetation in arid areas creates potentially conflicting demands for water between the ecosystem and humans [48]; these processes, although enhancing biomass storage or carbon sequestration, may decrease water availability [10,49]. In this case, the relationships between GRAIN-WR and NPP-WR will change from synergies into trade-offs. Under the warming and drying trend of regional climate [16], the continuous water demand by agriculture and humans may cause a decline in WR, and eventually result in vegetation degradation. Therefore, the implementation of a high-efficiency water saving agriculture (e.g., through drip irrigation or plastic film) is necessary to reduce water consumption rates and to improve cropland productivity, thereby improving other services as complementary or insurance measures for ecological rehabilitation [10].
According to our study, the inner relationship between driving forces should be cautious about. And relationship among drivers may display both trade-offs and synergies or interaction effect. For example, population variation and technical progress would affect the change of land use, and eventually lead to the change of ES supply in this region. Restricted by prohibition of open grazing policy since the early 2000s, herdsmen families in Uxin increased mechanical power in artificial pasture. Such changes in land management increased sheep breeding and reduced grazing pressure on nature grassland [40]. Farm families increased their investment in crop farming and cultivated land during the last three decades, which resulted in cropland expansion hot-spot areas and increase of grain production. It was founded in "Grain for Green" project region in Ningxia, China, the hot-spot areas of 37% of well managed land had contributed 57% of grain output [44]. Similarly, Pretty et al. (2006) found that farming systems adopting sustainability enhancing practices had increased productivity and improving the supply of ES [50].

Uncertainties in ES Assessment
Our study has some limitations and uncertainties. Input data accuracy is one of the sources of errors, such as those stemming from the 250-m MODIS data products and 1-km NOAA NDVI data used in the NPP calculation process. We used a regression method to bridge the resolution gap between these two data sources. However, due to the high contrast of sand land when in proximity to other land covers, uncertainties are unavoidable when using NOAA data, which have a coarse spatial resolution. Our analysis showed a higher variability of NPP in the period 1979-1999 than in the period 2000-2016. The models used in this study also have some limitations. SP was based on the RWEQ model, which was developed in the 1980s for the Great Plains region in the USA. The use of recommended parameters in our study is another source of uncertainty. Therefore, measurement data should be acquired in future research, to provide a basis for the verification of parameters and results from the local application of this model. This should allow for a more accurate evaluation of ecosystem services in the sandy area of grassland.

Conclusions
The analysis of long-term trends in ES is important to understand their responses to climate change and human activities. In this paper, we analyzed the trends of ES, and the changes in trade-offs and synergies between ES and their driving forces, in the period 1979-2016. Our study highlights the importance of considering long-term trends of ES and their drivers in understanding ES interactions. Our results demonstrate that the magnitude of ES trade-offs and synergies vary in the long period, and may even switch for some pairs of ES. Trade-offs and synergies also showed regional distribution characteristics. The analysis of the driving forces of ES relationship change found that climate factors, land use change, technological progress, and population are the main factors behind changes in ES and in the relationships between them.