Quantitative Assessment for the Dynamics of the Main Ecosystem Services and their Interactions in the Northwestern Arid Area, China

: Evaluating changes in the spatial–temporal patterns of ecosystem services (ESs) and their interrelationships provide an opportunity to understand the links among them, as well as to inform ecosystem-service-based decision making and governance. To research the development trajectory of ecosystem services over time and space in the northwestern arid area of China, three main ecosystem services (carbon sequestration, soil retention, and sand ﬁxation) and their relationships were analyzed using Pearson’s coe ﬃ cient and a time-lagged cross-correlation analysis based on the mountain–oasis–desert zonal scale. The results of this study were as follows: (1) The carbon sequestration of most subregions improved, and that of the oasis continuously increased from 1990 to 2015. Sand ﬁxation decreased in the whole region except for in the Alxa Plateau–Hexi Corridor area, which experienced an “increase–decrease” trend before 2000 and after 2000; (2) the synergies and trade-o ﬀ relationships of the ES pairs varied with time and space, but carbon sequestration and soil retention in the most arid area had a synergistic relationship, and most oases retained their synergies between carbon sequestration and sand ﬁxation over time. Most regional relationships of sand ﬁxation and soil retention with a weak trade-o ﬀ during the before-2000 stage turned into weak synergies at the after-2000 stage; (3) the ES pairs of carbon sequestration and sand ﬁxation had signiﬁcant time-lagged trade-o ﬀ relationships in most arid areas, and these time trade-o ﬀ s were generally maintained for 2–6 years, where a higher absolute value of the maximum time-lagged correlation coe ﬃ cient corresponded to a longer time-lag. This study could assist us in understanding the trade-o ﬀ and synergy of ESs and will provide quantitative zonal information for the ecosystem-services-based ecomanagement of arid areas.


Introduction
Ecosystem services (ESs) are defined as the benefits of the natural environment that sustain the survival of people and the products and services directly or indirectly through ecosystem functions. These services contribute greatly to human welfare and goods [1]. Generally, ESs fall into four main categories: provisioning services, supporting services, regulating services, and cultural services [2]. Lower-level classifications (widely recommended by previous research) include the system of Costanza et al. [3] and the scheme of Xie et al. [4].

Study Area
The northwestern arid area of China (34°20′10″ N ~ 49°10′30″ N, 73°33′27″E ~ 106°50′34″ E) lies to the north of the Kunlun Mountains-Qilian Mountains terrain belt and to the west of the Helan Mountains, with an area of 235 × 10 4 km 2 , around one quarter the size of China. This arid area includes the whole area of the Xinjiang Uygur Autonomous Region, the Alxa plateau of the Inner Mongolia Autonomous Region, and the Hexi corridor (i.e., towards the west of the Yellow River) of Gansu province, with a total of 134 counties. The study area has a temperate continental climate, except for high-altitude localities with highland mountain climates. Here, the strong evapotranspiration capability is 8~10 times as great as its extremely scarce precipitation (230 mm). Mountain ranges (M) surround the two basins, and extensive deserts (D) are interwoven with scattered oases (O), a unique type of topography denoted a mountain-oasis-desert system (MODS). Referring to the research of Wang et al. [22], our study divides the northwestern arid land into five large parts: northern Xinjiang (NX), southern Xinjiang (SX), the Yili basin (YL), eastern Xinjiang (EX), and the Alxa Plateau-Hexi Corridor area (APHC). Furthermore, mountains, oases, and deserts are extracted for each area based on the elevation range and land use type. Here, Yili is only zoned into mountains and plains. After these partitions, the study area is divided again into 14 sub-regions ( Figure 1).

Data and Preprocessing
Meteorological observations of the 86 stations in and around the study area were collected from the National Meteorological Information Center (http://data.cma.cn/). The meteorological elements include daily wind speed, maximum air temperature, minimum air temperature, sunshine hours,

Data and Preprocessing
Meteorological observations of the 86 stations in and around the study area were collected from the National Meteorological Information Center (http://data.cma.cn/). The meteorological elements Sustainability 2020, 12, 803 4 of 17 include daily wind speed, maximum air temperature, minimum air temperature, sunshine hours, relative humidity, and precipitation. After removing outliers (about 0.1%), the daily meteorological data over a half month and a whole month were averaged from records at day 15 and day 30 of the corresponding month. By applying AUNSPLIN 4.3, half-monthly and monthly grid meteorological data (5 km × 5 km) were generated by spatial interpolation. The Shuttle Radar Topographical Mission (SRTM) digital elevation model (DEM) data, with a 3 arc-second spatial resolution, was download from the Resource and Environment Data Cloud Platform (http://www.resdc.cn/). The gridded data of the soil particle diameter fraction was generated by the National Second Soil Survey of China, whose units were transformed from International Units (clay, silt, and sand: 0.002-0.02-2.000 mm) into American Units (clay, silt, and sand: 0.002-0.050-2.000 mm) using the Logarithmic normal distribution function. Besides soil texture data, the land cover map from the Climate Change Initiative (https://www.esa-landcover-cci.org/) (CCI-LC, v2.0), snow depth data (http://westdc.westgis.ac.cn/), and the LTDR (Land Long Term Data Record) v5 NDVI (https://ltdr.modaps.eosdis.nasa.gov/) were resampled into a 5 km × 5 km spatial solution.

Carbon Sequestration
In this study, NPP (gC·m −2 ·a −1 ) was selected as an indicator of carbon sequestration ESs. The terrestrial Carnegie Ames-Stanford Approach (CASA) model [23] was used to estimate NPP based on a light-use efficiency model. Compared with other NPP-models based on biophysical processes (e.g., VPM (Vegetation Photosynthesis Model) [24], BESS (Breathing Earth System Simulator Model) [25], and GLO-PEM (Global Production Efficiency Model) [26]), CASA demands fewer inputs and consumes less time for computing due to its simple and efficient algorithm. Its formulations are as follows: where APAR(x,t) is the photosynthetically absorbed active radiation (MJ·m −2 ) on the grid cell x in month t; ε is the light-use efficiency (gC·MJ −1 ) of the vegetation at location x at time t; SOL is the total solar radiation (MJ·m −2 ); and FPAR is the fraction of the photosynthetically active radiation absorbed by the canopy. Light-use efficiency refers to the canopy's efficiency in transforming absorbed radiation into organic carbon, as a process controlled by low temperature, high temperature, and wetness. Therefore, ε is allowed to be represented by the stress of low temperature and high temperature (T ε1 ), the stress effects of the gap between the varying ambient air temperature, and the optimal temperature closely related to the accumulation of NPP (T ε2 ), as well the wetness stress (W ε ). These three variables are dimensionless. ε max is the maximum possible efficiency. Here, the studies of Jiao Wei et al. [27] and Zhu Wenquan et al. [28] serve as the main references to calculate the key parameters (i.e., NDVImax, NDVImin, SRmax, SRmin, and ε max ) based on the 16 reclassified land-cover types and 8 land-use types from the CCI-LC map.

Soil Retention
The Revised Universal Soil Loss Equation (RUSLE) was employed to calculate soil retention. It can be formulated as where SC, SE P , and SE V represent the soil's retention capacity (t·km −2 ·a −1 ), potential soil loss (t·km −2 ·a −1 ), and actual soil loss (t·km −2 ·a −1 ), respectively; R is the rainfall erosivity factor (MJ·mm·km −2 ·h −1 ·a −1 ); K is the soil erodibility factor (t·km 2 ·h·km −2 ·MJ −1 ·mm −1 ); L, S, C, and P are the slope length factor, slope factor, vegetation factor, and management factor, and all factors are dimensionless. The value of P can be retrieved from different land use/covers and slopes. In this study, a reformulated scheme of the rainfall erosion factor from Yang Fengbo and Lu Changhe's study [29] was used to compute the R of the study area, which is described as where R is the rainfall erosivity in a half month (MJ mm ha −1 h −1 ), P j is the erosive rainfall (mm) for day j in a given half-month, and α and β are the empirical coefficients that are practical for the northern area of China. P d12 and P y12 are the average daily and annual rainfall for the days with rainfall ≥12 mm. The annual erosivity is the summation of all half-month erosivity during the year. (8) where sand, silt, clay, and oc are the content fraction of sand, silt, clay, and organic carbon in the top soil layer (%); sni = 1 − sand/100.
where λ, θ, and m are the slope length, slope, and slope index, respectively, while fvc is the vegetation coverage.

Sand Fixation
The Revised Soil Erosion Equation (RWEQ) model was used to simulate the regional wind erosion modulus (defined as the soil loss per unit area) [30] of the northwestern arid area of China. This model combines empirical and process modeling and has been extensively tested under broad field conditions [11,13,14].
Sustainability 2020, 12, 803 6 of 17 where ST max 0 is the maximum transport capacity (kg/m), WF is the climatic factor (kg/m), EF is the soil erodibility factor (%), SCF is the soil crust factor (dimensionless), K is the surface roughness (dimensionless), S 0 is the field length scale (m), and Z is the distance from the upwind edge of the field (m).
where VEG is the vegetation coverage factor; S is the critical field length, at which 63% of maximum transport capacity occurs. SF (sand fixation (t·hm −2 ·a −1 )) is the difference between potential soil erosion (SL 0 , t·hm −2 ·a −1 ) and actual soil erosion (SL, t·hm −2 ·a −1 ) per unit area by wind.

Trend Analysis
The linear slopes of the analytical factors (carbon sequestration, soil retention, and sand fixation) over time can reflect changing trends. The slope of the linear model was calculated for each ES at every time period (1990-2000, 2000-2015, and 1990-2015) by using the regress function of MATLAB R2016a. The formulas are as follows: where xi is the i th value of the ES time series for a certain time period, and n is the total number of years for a certain time period. A negative value of slope (α) represents the occurrence of a decrease, while a positive value of the slope means an increasing trend. Meanwhile, a t-test was used to identify the linear significance.

Correlation Analysis
To uncover the relationships of the ESs, a Pearson's correlation analysis (using the corr function in MATLAB R12a) was conducted among the pairs of ESs for each time period (1990-2000, 2000-2015, and 1990-2015). A positive correlation coefficient indicates that paired ESs have a synergistic relationship, whereas a negative correlation coefficient indicates a certain trade-off between the paired ESs.

Time-Lagged Cross-Correlation Analysis
A time-lagged cross correlation analysis was used to research the interaction time of each ES pair. With the assumption that every ESs pair is related to all other pairs at any time-lagged point, their time-lagged correlation coefficient can be represented as Sustainability 2020, 12, 803 where r k is the correlation coefficient at time-lagged point k, c k is the covariance of time series x t and time series y t+k , δ x is the mean square deviation of time series x t , and δ y+k is the mean square deviation of time series y t+k . These statistical parameters are calculated as follows: where the time-lagged k is equal to 0, ±1, ±2, . . . , ±n/4; y t+k is the mean of the time-lagged series y t+k , and x t is the mean of time series x t .

Spatial Pattern and Temporal Trends of ESs
The ecosystems services of NW from 1990 to 2015 had an average of 159.563 gC m −2 a −1 , 70.077 t km −2 a −1 , and 80.256 t hm −2 a −1 in annual carbon sequestration, soil retention, and sand fixation, respectively. The spatial patterns of the ecosystem services during the two sub-periods (i.e., before 2000 and after 2000) and over the whole study period are shown in Figure 2 However, the southern Xinjiang oasis and its desert, as well as the Alxa Plateau-Hexi Corridor desert, had the lowest soil retention. Besides being the highest and the lowest areas, the northern Xinjiang mountains, Alxa Plateau-Hexi Corridor mountains, and Yili plains experienced higher soil retention. In a word, most mountains had higher soil retention than deserts and oases. The boundary between the soil retention of the desert and that of the oasis, therefore, seems to disappear. As shown in Figure 2c Least-squares linear regression was used to predict the slope of each ES in every sub-region for different periods (Table 1, Figure 3). In terms of carbon sequestration, the whole arid area, mountains, and desert had similar temporal trends: The before-2000 series significantly increased (p < 0.01), and the after-2000 series significantly decreased (p < 0.05), with a slight variation over the whole period. Carbon sequestration of the oasis significantly increased during all periods. Carbon sequestration of the MODS for each large region resembled that of the MODS for the whole NW (mountains: increasedecrease; oasis: increase-increase; desert: increase-decrease) in the two sub-periods (Table 1, Figure  3(a1-a3)). In terms of soil retention, although the regional average of the oases and deserts decreased, the average increased in the mountains and for the whole NW during 1990-2015 ( Figure 3(b1-b3)). None of the trends of these regional averages passed a significance level of a = 0.05. As Table 1 details, during 1990-2015, soil retention of the oasis and deserts in SX, EX, APHC, and YL increased, but that of the MODS sub-regions in the NX decreased (Table 1). Before 2000, the soil retention of mountains and oases in all five large regions slightly declined, except for the MODS of NX, which showed increasing trends. After 2000, the soil retention of all MODS sub-regions insignificantly decreased, except for SX-O and APHC-D, which showed increasing trends. In terms of the sand fixation, mountains, oases, deserts, and the whole arid area all showed a monotonous downward trend (p < 0.01) during 1990-2015, with a slope of −2.080 t hm −2 a −2 , −2.591 t hm −2 a −2 , −4.535 t hm −2 a −2 , and −3.666 t hm −2 a −2 , respectively. Although most MODS sub-regions experienced significant upward trends during the whole study period, when viewed by the two sub-periods, half of the MODS sub-regions Least-squares linear regression was used to predict the slope of each ES in every sub-region for different periods (Table 1, Figure 3). In terms of carbon sequestration, the whole arid area, mountains, and desert had similar temporal trends: The before-2000 series significantly increased (p < 0.01), and the after-2000 series significantly decreased (p < 0.05), with a slight variation over the whole period. Carbon sequestration of the oasis significantly increased during all periods. Carbon sequestration of the MODS for each large region resembled that of the MODS for the whole NW (mountains: increase-decrease; oasis: increase-increase; desert: increase-decrease) in the two sub-periods (Table 1, Figure 3(a1-a3)). In terms of soil retention, although the regional average of the oases and deserts decreased, the average increased in the mountains and for the whole NW during 1990-2015 ( Figure 3(b1-b3)). None of the trends of these regional averages passed a significance level of a = 0.05. As Table 1 details, during 1990-2015, soil retention of the oasis and deserts in SX, EX, APHC, and YL increased, but that of the MODS sub-regions in the NX decreased (Table 1). Before 2000, the soil retention of mountains and oases in all five large regions slightly declined, except for the MODS of NX, which showed increasing trends. After 2000, the soil retention of all MODS sub-regions insignificantly decreased, except for SX-O and APHC-D, which showed increasing trends. In terms of the sand fixation, mountains, oases, deserts, and the whole arid area all showed a monotonous downward trend (p < 0.01) during 1990-2015, with a slope of −2.080 t hm −2 a −2 , −2.591 t hm −2 a −2 , −4.535 t hm −2 a −2 , and −3.666 t hm −2 a −2 , respectively. Although most MODS sub-regions experienced significant upward trends during the whole study period, when viewed by the two sub-periods, half of the MODS sub-regions in the entire NW area (i.e., the MODS of NX, EX and YL, and SX-M) continuously decreased, while the other half (i.e., the MODS of APHC, SX-O, and SX-D) experienced an increase-decrease process (Table 1, Figure 3(c1-c3)). Note: the * symbol on the top-right of the p value means that the tendency is significant (p < 0.05), and the ** symbol means that the tendency is very significant (p < 0.01).
in the entire NW area (i.e., the MODS of NX, EX and YL, and SX-M) continuously decreased, while the other half (i.e., the MODS of APHC, SX-O, and SX-D) experienced an increase-decrease process (Table 1, Figure 3(c1-c3)).

Spatial Patterns and Temporal Trends of the Relationships between ES Pairs
The interactions determined by regional average correlation coefficients are shown in Figure 4, alongside those that were checked according to significance (p < 0.05). The relationships between carbon sequestration and soil retention, carbon sequestration and sand fixation, and sand fixation and soil retention were classified into six categories: trade-off** (r < 0, p ≤ 0.01), synergy** (r > 0, p ≤ 0.01), trade-off*(r < 0, 0.01 < p ≤ 0.05), synergy* (r > 0, 0.01 < p ≤ 0.05), trade-off (r < 0, p > 0.05), and synergy (r > 0, p > 0.05). These trade-offs and synergies ( Figure 4) have great spatial and temporal heterogeneity throughout the whole NW. From the aspect of the relationship between carbon sequestration and soil retention, during 1990-2015, synergies were present in all sub-regions of NW, except APHC-D, with a strong trade-off (r < 0, p ≤ 0.05) (Figure 4a). Before 2000, weak synergies (r > 0, p > 0.05) were mostly spread across the whole arid area, with the exception of weak trade-offs (r < 0, p > 0.05) in the NX-D, EX-O, and APHC-D. After 2000, only SX-D and APHC-D had a weak trade-off relationship, while the other sub-regions established synergistic relationships between carbon sequestration and soil retention, which presented all the strong synergies for the oases and the whole YL (r > 0, p ≤ 0.05).

Time-Lagged Correlation of ESs Pairs
The corresponding correlation coefficients of carbon sequestration, soil retention, and sand fixation at each time-lag were calculated based on the principle of the time-lagged cross-correlation method. Then, the maximum time-lagged correlation coefficients were selected based on their significant linear relationships (p < 0.05), and their corresponding time-lags were regarded as the interaction time range among them. These results are illustrated by the radar charts ( Figure 5). According to the study of Li Boyan et al. [17], the time trade-off is the maximum time-lagged correlation coefficient that is significant and less than 0, and its corresponding time is the trade-off time; the time synergy is the maximum time-lagged correlation coefficient that is significant and greater than 0, and its corresponding time is synergistic time.
The maximum correlation coefficients of carbon sequestration and sand fixation had obvious regional differences (Figure 5a), and most regions (i.e., 18 out of 23 regions) passed a significance test among various spatial scales. Here, all oasis sub-regions, including the plains of YL and NW-O, had significant time trade-offs. Half of the mountainous sub-regions had a significant time trade-off like the mountains across the whole NW. Two out of the four desert sub-regions with significant relationships had time trade-offs, while the other two had time synergies. These differences among the oases themselves, among the mountains themselves, and among deserts themselves seem not to influence the time-lagged correlation between the carbon sequestration and sand fixation of MODS at the spatial scale of the whole NW, and they all presented time trade-offs with the maximum timelagged correlation coefficients of −0.704, −0.422, and −0.434, respectively. Correspondingly, the timelag of the oases ranged from 1-4 years (Figure 5b), and half of them were delayed by four years. All oasis sub-regions had longer time-lags than the oases for the whole NW (i.e., one year), which means that the time-lag of the whole NW does not equal or approximate the average of the values of all its components. For the mountainous sub-regions and desert sub-regions, the time-lag was mainly

Time-Lagged Correlation of ESs Pairs
The corresponding correlation coefficients of carbon sequestration, soil retention, and sand fixation at each time-lag were calculated based on the principle of the time-lagged cross-correlation method. Then, the maximum time-lagged correlation coefficients were selected based on their significant linear relationships (p < 0.05), and their corresponding time-lags were regarded as the interaction time range among them. These results are illustrated by the radar charts ( Figure 5). According to the study of Li Boyan et al. [17], the time trade-off is the maximum time-lagged correlation coefficient that is significant and less than 0, and its corresponding time is the trade-off time; the time synergy is the maximum time-lagged correlation coefficient that is significant and greater than 0, and its corresponding time is synergistic time.
The maximum correlation coefficients of carbon sequestration and sand fixation had obvious regional differences (Figure 5a), and most regions (i.e., 18 out of 23 regions) passed a significance test among various spatial scales. Here, all oasis sub-regions, including the plains of YL and NW-O, had significant time trade-offs. Half of the mountainous sub-regions had a significant time trade-off like the mountains across the whole NW. Two out of the four desert sub-regions with significant relationships had time trade-offs, while the other two had time synergies. These differences among the oases themselves, among the mountains themselves, and among deserts themselves seem not to influence the time-lagged correlation between the carbon sequestration and sand fixation of MODS at the spatial scale of the whole NW, and they all presented time trade-offs with the maximum time-lagged correlation coefficients of −0.704, −0.422, and −0.434, respectively. Correspondingly, the time-lag of the oases ranged from 1-4 years (Figure 5b), and half of them were delayed by four years. All oasis sub-regions had longer time-lags than the oases for the whole NW (i.e., one year), which means that the time-lag of the whole NW does not equal or approximate the average of the values of all its components. For the mountainous sub-regions and desert sub-regions, the time-lag was mainly concentrated at the one year, four year, and six year time ranges, and at the one year, two year, and six year time ranges, respectively.
A few of the regional maximum correlation coefficients between carbon sequestration and soil retention passed the given significance level, including NX-M (0.004), APHC-M (0.005), APHC-O (0.007), NW-O (0.010), and NW (0.017) (Figure 5a). All these coefficients showed weak time synergies due to their values being much less than 1. The relationships among mountains and among oases displayed great spatial differences, and there were few significant interactions in the mountains (i.e., 2 out of 6 mountainous sub-regions), few significances in the oases (i.e., 1 out of 5 oasis sub-regions), and no significant relationship in the deserts. However, based on the scale of the whole NW, both NW-O and the whole NW had significant time synergies. Moreover, only APHC-O had a 1-year time lag among these significant time synergy regions (Figure 5b). There was no significant difference among all regional maximum correlation coefficients between soil retention and sand fixation, regardless of the sub-regions, parts, or MODS of the whole arid area (Figure 5a). Sustainability 2020, 12, x 12 of 17 concentrated at the one year, four year, and six year time ranges, and at the one year, two year, and six year time ranges, respectively. A few of the regional maximum correlation coefficients between carbon sequestration and soil retention passed the given significance level, including NX-M (0.004), APHC-M (0.005), APHC-O (0.007), NW-O (0.010), and NW (0.017) (Figure 5a). All these coefficients showed weak time synergies due to their values being much less than 1. The relationships among mountains and among oases displayed great spatial differences, and there were few significant interactions in the mountains (i.e., 2 out of 6 mountainous sub-regions), few significances in the oases (i.e., 1 out of 5 oasis sub-regions), and no significant relationship in the deserts. However, based on the scale of the whole NW, both NW-O and the whole NW had significant time synergies. Moreover, only APHC-O had a 1-year time lag among these significant time synergy regions (Figure 5b). There was no significant difference among all regional maximum correlation coefficients between soil retention and sand fixation, regardless of the sub-regions, parts, or MODS of the whole arid area (Figure 5a).

Explanation of the Interactions between ESs
The spatial scale was not only reflected by the trend of each ES but also by relationships among ESs. At the sub-regional scale, most relationships among ESs in the MODS of the whole NW were not consistent with those in the MODS of each large region; this difference varied over time and space. For example, during 1990-2015, trade-offs between carbon sequestration and soil retention were found for the oases and deserts of the APHC and the mountains of the APHC with synergy, but synergies occurred at the MODS for the whole NW. Therefore, interactions between ESs at microscopic scales may be hidden at the macroscopic scale. Although synergies between carbon sequestration and soil retention occurred in most areas of the NW and during most periods, continuous trade-offs happened in APHC-O and EX-O. These conflicting relationships between carbon sequestration and soil retention could be explained by the constraint threshold theory, proposed by Hao Ruifang et al. [31], where this relationship is characterized by a hump-shaped line. Higher carbon sequestration indicates good vegetation coverage that could protect soil from erosion

Explanation of the Interactions between ESs
The spatial scale was not only reflected by the trend of each ES but also by relationships among ESs. At the sub-regional scale, most relationships among ESs in the MODS of the whole NW were not consistent with those in the MODS of each large region; this difference varied over time and space. For example, during 1990-2015, trade-offs between carbon sequestration and soil retention were found for the oases and deserts of the APHC and the mountains of the APHC with synergy, but synergies occurred at the MODS for the whole NW. Therefore, interactions between ESs at microscopic scales may be hidden at the macroscopic scale. Although synergies between carbon sequestration and soil retention occurred in most areas of the NW and during most periods, continuous trade-offs happened in APHC-O and EX-O. These conflicting relationships between carbon sequestration and soil retention could be explained by the constraint threshold theory, proposed by Hao Ruifang et al. [31], where this relationship is characterized by a hump-shaped line. Higher carbon sequestration indicates good vegetation coverage that could protect soil from erosion by water or indicates more precipitation, which could increase the possibility of soil erosion. These multiple factors may jointly affect the relationship of carbon sequestration and soil retention. During 1990-2015, the carbon sequestration of APHC-O and EX-O exceeded the constraint threshold, which restricted the improvement of soil retention [20]. Meanwhile, the relationship among multiple ecosystem services reveals some degree of evolution, and this development is multistage and distinct [32,33]. For example, the relationship between carbon sequestration and soil retention in NX-D converted from a trade-off to synergy relationship. Carbon sequestration with the other two ESs both have time-lags; these delayed effects may result from the changes of vegetation coverage in the study area, based on the historical influence of land use/cover on the current provisions of ESs [13,34]. Further, a greater absolute value of the maximum time-lagged correlation coefficient corresponded to a longer time-lag.

Implications of Ecosystem Services Assessment for Ecological Management
The carbon sequestration of the most arid area has improved over the past 26 years and has shifted in the mountains and deserts around the year 2000. The Tianshan Mountains-Altai Mountains Natural Forest Protection Project launched in 2000, effectively protecting state-owned natural forests by expanding ecological forest land based on rational zoning, changing the business direction of forest industry enterprises to silvicultural treatments, and reducing the timber output [35]. Carbon sequestration of the whole arid area changed around 2000, possibly because of the increase in biological sand barriers to combat sand movement in the deserts of APHC. Although the protection of desert Populus diversifolia was strengthened in 1980, most areas of the protected region were relatively small, such as gardens or belts alongshore rivers or lakes. However, the declining carbon sequestration of desert eco-regions needs to be paid more attention to, particularly EX-O and SX-O. Carbon sequestration reflects vegetation density and growth conditions, and vegetation coverage has degraded severely across these two regions [20]. Therefore, the key to increasing the ability of carbon sequestration is an improvement of vegetation coverage. Increasing the amount of water availability (e.g., precipitation, evapotranspiration, soil wetness, and humidity) is a widely accepted and recommended approach [20,36]. Due to a deficiency in natural water, eco-environment water (3.6%) is used by agricultural irrigation (about 92%-95%) [37]. Breeding native plant species (e.g., Eurotia arborescens ALos, Caragana Korshinsii Kom, Psammochloa mongolica, and Calligunum mongolicum Turcz) and applying drip irrigation technology may improve this degeneration of future ecosystem management. Due to the declined provision of soil retention and the warming climate in the NW [38], more runoff has been yielded by melting high-altitude glaciers (e.g., the Tian Shan in the Xinjiang) and snow water flows to the plains or basins [39]; this process can remove soil and its fertilizer. In addition, increasing the area of farmland results in greater damage to the soil's surface after tillage. Wind-break and sand-fixation efforts in the APHC in the 1990s, may have caused an increase of sand fixation in the APHC, which was mainly supported by the implementation of the Gansu Integrated Desert Control and Sustainable Agriculture Project [40]. This project combined physical sand-fixation technologies (e.g., straw checkerboard barriers and an enclosure zone for natural vegetation) and biological technologies (e.g., a farmland shelterbelt net and a vegetation arrangement for windbreak). To encourage people to positively participate in ecological work, the plants selected as vegetation to prevent land desertification could include drought-enduring medicinal herbs with excellent self-survival or grass for forage. The considerable economic benefits of these plants could stimulate participants to continue the work. As one of the sources of sandstorms, the arid area of northwestern China carries sand and greatly impacts the sandstorm frequency downwind of Beijing and Tianjin, but the Badain Jaran Desert (i.e., the desert to the west of the Alxa desert, Inner Mongolia) was not covered by the second round of the improvement project of the Beijing-Tianjin sandstorm source region [41]. In the future, plans for this arid area may improve the effectiveness of controlling sandstorm projects. For the whole northwest arid area, availability water and existing vegetation are both important for the ecological environment. Therefore, the punishments for damaging natural forests and polluting water could be increased in severity. Meanwhile, raising the price of agricultural water could restrict water-wasting behaviors. Ecological recovery and reconstruction need to be optimized stepwise, along with problem diagnoses, planning, implementation, monitoring, and evaluation [42]. Considering the dynamic nature of the ecosystem's components, structure, progress, and functions, the evolution of the ecosystem and landscape present various uncertainties in ecological recovery applications, even those under artificial control. The costs of each plan, and their participants, also influence the implementation of ecological projects. Our results only offer some preliminary suggestions from the perspective of ecosystem services. The dynamic and continuous multi-aspect (e.g., investment and costs, willingness of participants, and ecosystem quality) activities of supervision, analysis, and assessment should be further enhanced.

Limitations and Improvements Methods
We recognize some limitations of this study. The first lies in the linear correlation analysis of the relationship among ESs, which assumes that the relationships between ecosystem services are monotonic. However, considering the complexity and diversity of the ecosystem process, there may be more types of interactions among them, such as constraint effects with diverse shapes [31], as well as trade-offs and synergistic relationships depicted by a bagplot [43]. Therefore, it is necessary to adopt some nonlinear statistical analyses (e.g., a segmented linear regression, an accumulative probability density curve, and logarithmic function fitness) to explore these relationships. The second limitation is that the effects of spatial scale have been shown in the current research, but a finer spatial scale would allow land-use policy-makers and ecological planners to acquire more explicitly spatial information. In future research, ESs and their relationships with the others, as well as their associations with potential contributors (like climate change and anthropic activities), could be assessed at a grid-cell scale. Moreover, more detailed ecological levels (like ecosystem class and landscape level) could be evaluated to enrich our understanding of ESs. The third limitation is that errors in the regionalization of MODS may influence determinations of the relationship type between pair-wised ESs. On-the-spot investigations could be conducive to reducing error and may help further reveal the regular patterns in different eco-regions. Last but not least, there are problems in the validations of the ESs estimated by the biophysical progress model; these difficulties result from the lack of observations in the study regions and the inconsistencies between observations and scaled extrapolation results. This is an important and challenging matter for ES assessment.

Conclusions
In summary, our study provided a quantitative assessment of the dynamics of ecosystem services (i.e., carbon sequestration, soil retention, and sand fixation) and their interrelationships over space and time in the arid area of northwestern China. We found that the ecosystem services of most of the arid area declined. Carbon sequestration in the oases continued to increase, and sand fixation of the Alxa Plateau-Hexi Corridor area, before 2000, improved markedly due to the great efforts of wind-break and sand-fixation. The trade-offs between carbon sequestration and soil retention mostly affected all periods and all sub-regions. The time-lagged trade-offs between carbon sequestration and sand fixation occurred widely in the arid area, and a greater maximum correlation coefficient resulted in a longer time-lag. These results highlight that interactions between ecosystem services vary over temporal and spatial scales. The enhancement of carbon sequestration in deserts and soil retention and sand fixation in most arid areas could be future targets for ecological restoration.
Author Contributions: T.K. designed the research and carried out the results analysis with T.K. and S.Y.; J.B. collected meteorological data and filled the gap; J.C. pretreated the remote sensing data; Y.G. was the academic adviser. All authors took part in the writing and revision of the manuscript. All authors read and approved the final manuscript.