Spatiotemporal Patterns of Desertification Dynamics and Desertification Effects on Ecosystem Services in the Mu Us Desert in China

Abstract

Degradation of semi-arid and arid ecosystems due to desertification is arguably one of the main obstacles for sustainability in those regions. In recent decades, the Mu Us Desert in China has experienced such ecological degradation making quantification of spatial patterns of desertification in this area an important research topic. We analyzed desertification dynamics for seven periods from 1986 to 2015 and focused on five ecosystem services including soil conservation, water retention, net primary productivity (NPP), crop productivity, and livestock productivity, all assessed for 2015. Furthermore, we examined how ecosystem services relate to each other and are impacted by desertification. Three major conclusions are drawn from the study. First, the eastern part of the study area experienced overall improvement while desertification in the west first increased and then reversed its trend during those periods between 1986 and 2015. Second, significant synergistic relationships are observed for three regulating services (soil conservation, water retention, NPP) and two provisioning services (crop productivity and livestock productivity). Strong relationships across different types of ecosystem services were found only between crop productivity and NPP. Third, in response to increasing desertification, the three regulating services exhibit a monotonically decreasing trend, while the two provisioning services follow a hump-shaped response.

1. Introduction

People benefit from ecosystems in the form of provisioning (food, water, wood, and fiber), regulating (climate, floods, diseases, waste, and water quality), and supporting (soil formation, photosynthesis, and nutrient cycling) ecosystem services, on one hand, and entertainment, aesthetic enjoyment and spiritual benefits and other cultural needs, on the other [1]. All those services are important for maintenance and improvement of human well-being and, eventually, for the sustainability of regions [2]. The Millennium Ecosystem Assessment [1] reports that in the last 50 years 15 out of 24 ecosystem services have been degrading worldwide. Such degradation will directly threaten the safety and health of human populations, as well as regional and global ecological security [1,3]. Induced by a combination of climate change and human activities, desertification has become the main reason for ecosystem degradation in many dry regions. Rightfully so, it was named one of the main impediments in sustainable trajectories of development of those regions by the Ordos Declaration issued by the United Nations Convention to Combat Desertification’s 13th conference [4]. If desertification continued to occur at the current rate, it could seriously endanger and jeopardize the existence of approximately 350 million people [5]. It is, therefore, important to study the impact of desertification on ecosystem services, which will have far-reaching implications for regional sustainability.

Due to the diversity of ecosystem services and the spatial-temporal heterogeneity of land uses, their interactions are not only characterized by a combination of natural environmental dynamics and human decisions and experiences but also by reciprocal tradeoffs and synergies [6]. For example, when availability of ecosystem services is highly constrained by natural conditions, the increase in provisioning services often results in the reduction of regulating and cultural services thus exhibiting trade-offs [1,7]. With respect to ecosystem services, synergy refers to the simultaneous increase or reduction of two kinds of ecosystem services. For example, biodiversity conservation can cause an increase in provisioning services [8,9]. There are two main types of maintenance mechanisms for tradeoffs or synergies in ecosystem services. The first one is based on a single factor driving multiple ecosystem services. Changes in this driving factor will cause feedbacks from ecosystem services, resulting in trade-off or synergy between them. The second mechanism stems from interactions between ecosystem services. Changes in one ecosystem service will be immediately reflected in the other [10]. Therefore, to achieve maximum benefits from ecosystem services, we need to understand the relationships between multiple services and use the knowledge to anticipate reciprocal effects rather than just pursuing the revenue from one service [11].

Desertification is referred to as land degradation triggered by climate change and human activities in arid and semi-arid regions. Such degradation takes form of physical, chemical, and biological degradation of soil characteristics and the long-term loss of natural vegetation [12,13]. The latter is often used as the indicator of desertification in addition to documenting the reduction in soil depth due to erosion and in the decline in productive capacity and water-holding properties of the soil. Desertification is accompanied by the decrease of soil nutrients, such as nitrogen, phosphorus, and potassium resulting in the decline of soil fertility [14,15]. Another negative process is the decline in soil water holding capacity, which directly impacts many ecosystem services. On broader scales, desertification can lead to the change in productivity of the whole region by affecting regional photosynthesis by dust particles suspended in the air. As a result, soil conservation, water retention and net primary productivity become a matter of concern [16]. Clearly, desertification has mainly negative impacts on provisioning of ecosystem services [17].

To address desertification trends and their impact on ecosystem services we chose to study the Mu Us Desert transitioning from the Ordos Plateau to the Loess Plateau in China. It is characterized as an ecotone of pasture, forest, and agriculture, mainly pasture grazing. The area is deemed to be quite fragile and sensitive to environmental changes [18]. Being exposed to effects of recent climate changes the region has experienced the development of agriculture and animal husbandry and the increase in energy production, among other human activities, which contribute to the accelerated desertification in the area. While those threats and their effects have been extensively studied in the Mu Us Desert, only a few studies focused on the effects of desertification on ecosystem services. Our study was designed to simultaneously assess the process of desertification and dynamics of five ecosystem services—soil conservation, water retention, net primary productivity (NPP), crop productivity, and livestock productivity—in the area. This allowed us to assess the relationships between desertification and these ecosystem services, which we consider a contribution to the theory and analysis of ecosystem services in dry regions. The work is important in that it provides scientific knowledge for controlling desertification and achieving regional sustainability in the Mu Us Desert.

2. Materials and Methods

2.1. Study Area

The Mu Us Desert (106°11′–110°54′ E, 36°49′–40°12′ N) cuts across Shaanxi Province, Ningxia Hui and Inner Mongolia Autonomous Regions and includes eleven counties (banners) and districts (Figure 1). The climate of the area is temperate arid, semi-arid continental monsoon climate with the annual average temperature of 6–9 °C. The average annual precipitation increases from about 250 mm in the northwest to 400 mm in the southeast. Elevation generally changes from 1000 to 1600 m following the northwest–southeast trend. The variable topography includes such landforms as active dunes, semi-fixed dunes, and fixed dunes. The main soil contains Kastanozems, Arenosols, Histosols and Solonchaks. The zonal vegetation is dominated by Stipa bungeana and Thymus serpyllum plant communities, and sand vegetation is dominated by Artemisia ordosica, Salix psammophyla and Sabina vulgaris.

2.2. Data Sources

Based on field investigation, the vegetation coverage in the study area was obtained at the end of July 2016; a total of 71 sites were set up (Figure 1), three 10 × 10-m sample plots were established within each site, all the shrubbery clump width was measured, and the vertical projection area of each plant was calculated by the ellipse calculation method. Five 10 × 10-m quadrats were established within each sample plot, using digital cameras to shoot photos vertically above 1 m, ENVI 5.3 software (Exelis Visual Information Solutions; Boulder, CO, USA; 2015) is used to extract the herbaceous coverage. The sum coverage of shrubs and herbaceous as the measured coverage of one sample plot, the average value of coverage of three sample plots as the measured coverage.

63 Landsat images (TM and OLI, 30-m spatial resolution) from 1986 to 2015 (1986, 1990, 1995, 2000, 2005, 2010, 2015; path/row: 127/32-34, 128/32-34, 129/32-34; http://landsat.usgs.gov) were used in the study. All images were selected at the end of the July, due to the limits of cloud and time resolution, some images have to be replaced by near time: in 1995, images of 128/32-34 were replaced by the same period in 1996; in 2000, image of 128/34 was replaced by the same period in 1999; in 2010, images of 128/32 and 128/33 were replaced by the same period in 2009. Radiometric calibration, atmospheric, topographical, and geometric corrections were performed using ENVI 5.3 software with Digital Elevation Model (DEM). We obtained land use from the global land cover maps at 300-m spatial resolution published by the European Space Agency (ESA). DEM data were obtained from ASTER satellite images at 30-m spatial resolution (http://www.jspacesystems.or.jp/en_/). Data on precipitation and temperature were provided by the China Meteorological Data Service Center. Maps of soil type, soil texture, and soil depth were prepared by the Inner Mongolia Grassland Ecosystem Research Station, Chinese Academy of Sciences. Crop productivity and livestock productivity are from published statistical yearbook.

2.3. Data Analysis

2.3.1. Climate Change in Mu Us Desert from 1985 to 2013

Due to the lack of data from some meteorological stations, finally we obtained annual precipitation and mean annual temperature data from five meteorological stations in 1985–2013, two meteorological stations in the Inner Mongolia Autonomous Region (in Otog and Wushen banner), two meteorological stations in the Shaanxi Province (in Yuyang district, Hengshan county) and one in the Ningxia Autonomous Region (in Yanchi county). The distance between five meteorological stations is far enough and is evenly distributed in the study area. Pearson correlation analysis was conducted to determine the differences of trends in climate change between meteorological stations in Mu Us Desert for 1985–2013.

2.3.2. Desertification Difference Index (DDI)

The DDI was developed to quantify the degree of desertification based on remotely sensed data. The index combines the amount of vegetation with albedo, an important characteristic relevant to the surface radiation exchange [19,20]:

DDI = 2.234·NDVI − Albedo

where NDVI is the Normalized Difference Vegetation Index, which is often used as a proxy for vegetation cover or biomass. It is calculated as follows:

$$\mathrm{NDVI}= \frac{NIR-RED}{NIR+RED}$$

where NIR and RED correspond to spectral reflectance in the near-infrared band and the red band of a multispectral satellite image, respectively.

Albedo used in DDI calculation is the land surface albedo derived from satellite data. We used for following empirical formulas developed for Landsat imagery [21]:

Albedo (TM) = 0.356·ρ1 + 0.130·ρ3 + 0.373·ρ4 + 0.085·ρ5 + 0.072·ρ7 − 0.0018
Albedo (OLI) = 0.356·ρ2 + 0.130·ρ4 + 0.373·ρ5 + 0.085·ρ6 + 0.072·ρ7 − 0.0018

where ρ1, ρ2, ρ3, ρ4, ρ5, ρ6, ρ7 are spectral reflectance of corresponding Landsat TM and OLI spectral bands.

The range of NDVI and Albedo was normalized to 0–100:

N = [(NDVI − NDVI_min)/(NDVI_max − NDVI_min)] × 100

A = [(Albedo − Albedo_min)/(Albedo_max − Albedo_min)] × 100

To investigate the transformation of desertification in two-dimensional feature space constituted by Albedo and NDVI, the scatter diagram of Albedo-NDVI feature space is constructed (Figure 2a), and Zeng et al. [20] summarized the all situation of desertification to Figure 2b. In Albedo-NDVI feature space, Albedo is not only a function of vegetation coverage but also a function of soil moisture content. In Figure 2b, point A represents dry bare soil, point B represents wet bare soil, point C represents a high vegetation covered area with low soil moisture content, point D represents high vegetation covered area, and the soil moisture content is high [20].

To obtain the quantitative relationship between Albedo and NDVI under different desertification, the statistical regression analysis was carried out by using the Albedo and NDVI values of the 71 field samples, and various degrees of desertification can be found in these samples. The results show that there is a significant negative linear relationship between Albedo and NDVI (Figure 3):

Albedo = −0.44762 NDVI + 48.344

According to Verstraete and Pinty’s study [22], if the Albedo-NDVI feature space is divided on the vertical direction of desertification change trend, different desertification land can be separated effectively, and the desertification can be expressed by a binary linear equation:

DDI = a·NDVI − Albedo

where DDI is Desertification difference index, a is the reciprocal of the slope in the equation, so in this study the a is 2.234.

The accuracy of obtained DDI was compared with the other 9 desertification monitoring indices, it was verified by field measured vegetation coverage data through correlation analysis (

Table 1. Correlation analysis of measured vegetation coverage and 10 desertification monitoring indices.

	DDI	Albedo	SBI	SFOI	FVC	NDVI	WI	GI	LST	RUE
MVC	0.721 **	−0.714 **	−0.709 **	−0.656 **	0.583 **	0.580 **	0.566 **	0.393 **	−0.232	0.012

** represents highly significant correlation. MVC is measured vegetation coverage; Albedo is surface albedo; SBI is soil brightness index; SFOI is soil ferric oxide index; FVC is fractional vegetation coverage; NDVI is normalized difference vegetation index; WI is wetness index; GI is greenness index; LST is land surface temperature; RUE is rain use efficiency.

).

2.3.3. Quantification of Ecosystem Services in the Mu Us Desert

Soil Conservation (SC) was assessed as the difference between potential and actual soil losses calculated using the Revised Universal Soil Loss Equation (RUSLE) as follows [23]:

A_p = R·K·LS

where A_p is the potential soil loss (t/ha·y); R is the rainfall erosivity factor (MJ·mm/ha·h·yr):

$$\mathrm{R}={\displaystyle \sum }_{i=1}^{12}\left[1.735\times {10}^{\left(1.5lg\frac{P_i^2}{P}-0.8188\right)}\right]$$

where P is annual precipitation (mm); P_i is i month preipitation (mm).

K is the soil erodibility factor (t·ha·h/ha·MJ·mm):

$$\begin{array}{c}\mathrm{K}=\left\{0.2+0.3exp\left[-0.0256SAN\left(1-\frac{SIL}{100}\right)\right]\right\}{\left(\frac{SIL}{CLA+SIL}\right)}^{0.3}\\ \left(1-\frac{0.25C}{C+exp\left(3.72-2.95C\right)}\right)\left(1-\frac{0.7SN}{SN-exp\left(-5.51+22.9SN\right)}\right)\end{array}$$

where SAN, SIL, CLA and C is percentage content (%) of sand (0.05~2 mm), silt (0.002~0.05 mm), clay (<0.002 mm) and organic carbon in soil; SN = 1 − SAN/100.

LS is the topographic factor:

$$\mathrm{L}={\left(\frac{\lambda }{22.1}\right)}^{\alpha }$$

$$\mathrm{S}=\{\begin{array}{c}10.8\sin \theta +0.03\hspace{1em}\theta <9°\\ 16.8\sin \theta -0.50\hspace{1em}9°\le \theta <14°\\ 21.91\sin \theta -0.96\hspace{1em}\theta \ge 14°\end{array}$$

α = β/(1 + β)

$$\beta =\left(\frac{\sin \theta }{0.0896}\right)/\left(3.0\sin {\theta }^{0.8}+0.56\right)$$

where L is slope length factor; λ is horizontal slope length; S is slope factor; θ is slope extract from DEM; α is slope length index; β is the ratio of rill erosion to surface erosion.

A_r is the actual soil loss (t/ha·y):

A_r = R·K·LS·C·P

where C is the vegetation cover and management factor:

$$\begin{gathered} \mathrm{C}=\{\begin{array}{c}\hspace{1em}1\hspace{1em}0\le fc<0.1\%\\ 0.6508-0.3436\lg \left(fc\right)\hspace{1em}0.1\%\le fc<78.3\%\\ 21.91\sin \theta -0.96\hspace{1em}fc\ge 78.3\%\end{array} \\ fc=\frac{NDVI-NDV{I}_{min}}{NDV{I}_{max}+NDV{I}_{min}} \end{gathered}$$

where fc is vegetation coverage (%); P is the conservation practices factor, refers to the ratio of soil loss of sloping land under specific soil conservation measures to soil loss with no soil conservation measures, which reflects the inhibition of soil erosion.

Soil conservation (A_c) can be estimated as follows:

A_c = A_p − A_r = R·K·LS·(1 − C·P)

Water Retention (WR)

Water retention was estimated using the water yield sub model of the InVEST model [24], the equation is:

$$\mathrm{WR}=\mathrm{Min}(1,\frac{249}{Velocity})·\mathrm{Min}(1,0.3\mathrm{TI})·\mathrm{Min}(1,\frac{K_{sat}}{300})·\mathrm{WY}$$

where Velocity is the flow coefficient, TI is the topographic index extract from DEM, K_sat is the soil saturated hydraulic conductivity, and WY is the water yield, which is defined as the amount of water calculated using the InVEST model. The formula for the water yield model is based on the Budyko curve and the annual average precipitation [25]:

$${\mathrm{WY}}_{xj}=(1-\frac{AE{T}_{xj}}{P_x})·P_x$$

where WY_xj is the annual water yield at pixel x of land cover type j, AET_xj is the annual actual evapotranspiration for pixel x with land cover j, P_x is the annual precipitation for pixel x, AET_xj/P_x is based on the expression of the Budyko curve given as:

$$\frac{AE{T}_{xj}}{P_x}=1+\frac{PE{T}_{xj}}{P_x}-{\left[1+{\left(\frac{PE{T}_{xj}}{P_x}\right)}^{\omega }\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\omega $}\right.}$$

where PET_xj is the potential evapotranspiration for pixel x with land cover j, w is a non-physical parameter that characterizes the natural climatic-soil properties, which is defined as:

$${\omega }_x=\mathrm{Z}·\frac{AW{C}_x}{P_x}+1.25$$

where AWC_x is the volumetric (mm) plant available water content at pixel x, which is estimated using soil texture, soil depth, and root depth of vegetation, Z is the seasonal rainfall factor known as Zhang coefficient, which represents annual precipitation distribution and can be any real number in the range from 1 to 30. In our study, area, we set Z to 7, according to the InVEST User’s Guide.

Net Primary Productivity (NPP)

We used the CASA model to estimate NPP [26] based on the intercepted photosynthetically active radiation (IPAR) and light utilization efficiency (ε) as:

$$\begin{array}{c}\mathrm{NPP}\left(\mathrm{x},\mathrm{t}\right)=\mathrm{IPAR}\left(\mathrm{x},\mathrm{t}\right)·\epsilon \left(\mathrm{x},\mathrm{t}\right)\\ \mathrm{IPAR}\left(\mathrm{x},\mathrm{t}\right)=\mathrm{SOL}\left(\mathrm{x},\mathrm{t}\right)·\mathrm{FPAR}\left(\mathrm{x},\mathrm{t}\right)·0.5\\ \mathrm{FPAR}\left(\mathrm{x},\mathrm{t}\right)=\frac{NDVI\left(\mathrm{x},\mathrm{t}\right)-NDV{I}_{i, min}}{NDV{I}_{i, max}-NDV{I}_{i, min}}\left({\mathrm{FPAR}}_{max}-{\mathrm{FPAR}}_{min}\right) +{\mathrm{FPAR}}_{min}\end{array}$$

where t is time, and x is the grid cell; SOL(x, t) is the total solar radiation (gC·m⁻²·month⁻¹) in t month at pixel x; FPAR(x, t) is the absorption ratio of the vegetation layer to the incident photosynthetic effective radiation; the values of FPAR_max and FPAR_min were independent of the vegetation types, which were 0.001 and 0.95, respectively.

ε(x, t) = T_ε1(x, t)·T_ε2(x, t)·W_ε(x, t)·ε_max

where T_ε1(x, t) and T_ε2(x, t) indicate the stress effects of low temperature and high temperature on the utilization of light energy; W_ε(x, t) is the coefficient of water stress, reflecting the effect of water conditions; ε_max is the maximum utilization rate of the light energy under the optimum conditions.

Crop Productivity (CP) and Livestock Productivity (LP)

CP and LP were calculated from statistical data. We obtained the crop production and livestock production per district and per county (banner) from the 2015 statistical yearbook and determined the spatial distribution of CP and LP in the study area. For statistical purposes we converted all large livestocks into sheep units.

From the 2015 statistical yearbook, only a specific value of CP and LP is obtained in each county, unable to get the pixel size of the spatial distribution in the entire area, so these two services are limited to the county scale. For the five ecosystem services to be analyzed together, we extend raster layers of DDI, SC, WR, and NPP to county scales, that is, calculate the mean value of all pixels in each county as the value of ecosystem service. So, we show the values of ecosystem services on a county scale on the map. Spatial patterns of DDI and ecosystem services were analyzed for the period of 1986–2015 in each of the 11 districts and counties (banners) of the Mu Us Desert using ArcGIS 10.0. Pearson correlation analysis was conducted using SPSS 21 to determine the relationship between ecosystem services and the impact of desertification (2015) on ecosystem services [27], and the differences of trends in climate change between meteorological stations in Mu Us Desert for 1985–2013.

3. Results

3.1. Climate Change in Mu Us Desert from 1985 to 2013

The annual precipitation and mean annual temperature in Mu Us Desert monitored by each meteorological station are fluctuating in 1985–2013 (Figure 4). Whether annual precipitation or mean annual temperature, the trends of change in each meteorological station are consistent (

Table 2. Correlations between annual precipitation.

	Otog	Wushen	Yuyang	Hengshan	Yanchi
Otog	1
Wushen	0.744 **	1
Yuyang	0.403 *	0.744 **	1
Hengshan	0.447 *	0.769 **	0.869 **	1
Yanchi	0.377 *	0.646 **	0.585 **	0.673 **	1

* represents significant correlation; ** represents highly significant correlation.

and

Table 3. Correlations between mean annual temperature.

	Otog	Wushen	Yuyang	Hengshan	Yanchi
Otog	1
Wushen	0.975 **	1
Yuyang	0.949 **	0.963 **	1
Hengshan	0.966 **	0.961 **	0.935 **	1
Yanchi	0.820 **	0.824 **	0.729 **	0.835 **	1

** represents highly significant correlation.

).

3.2. Dynamic of Spatial Pattern of Desertification from 1986 to 2015

As a whole, desertification reversed during the past 30 years in the center and the east of the Mu Us Desert, while it was exhibiting a growth and intensification in the western part of the region (Figure 5). Desertification decreased in two of the four banners of the Inner Mongolia Autonomous Region area, including the Ejin Horo banner and Wushen banner. The decreasing trend was relatively slow in the Wushen banner, and no change was observed in the Otog banner and Otog Front banner, where small-scale intensification and reversal can be still detected. The five counties (districts) of Shaanxi in the southeastern corner of the Mu Us Desert have exhibited gradual decrease in desertification since 1986. In Shenmu county, Yuyang district, Hengshan county and Jingbian county, the recovery is particularly remarkable. Since 2005, no significant desertification has been observed except for a minor increase in Dingbian county. The part of Mu Us Desert, which is in Ningxia Autonomous Region, first shows the increase in desertification followed by the reversed pattern. Desertification in Lingwu City and Yanchi County was the strongest in 2005 after which it exhibited the decline during the 2005–2015.

3.3. Spatial Pattern of Ecosystem Services

Except for LP (Figure 6f), spatial patterns of the other four ecosystem services are generally consistent with each other. SC, WR, NPP, and CP decrease gradually from the southeast to the northwest of the Mu Us Desert (Figure 6b–e). Among them, the highest value of SC is in Shenmu county and the lowest value is in Otog banner (Figure 7b). The highest value of WR is also in Shenmu county, but the lowest value is in Lingwu city (Figure 7c). Hengshan county has the highest NPP and the Otog banner is characterized by the high productivity (Figure 7d). The highest values of CP and LP are both in Jingbian county and lowest are in the Otog banner and Ejin Horo banner (Figure 7e,f).

3.4. Relationships between Ecosystem Services

Synergies among ecosystem services are clearly observed for regulating services SC, WR and NPP, provisioning services, CP, and LP, but only NPP and CP are correlated when the link between regulating and provisioning services is considered (

Table 4. Relationships between ecosystem services.

	Soil Conservation	Water Retention	NPP	Crop Productivity	Livestock Productivity
Soil Conservation	1
Water Retention	0.961 **	1
NPP	0.897 **	0.928 **	1
Crop Productivity	0.426	0.585	0.725 *	1
Livestock Productivity	0.334	0.389	0.539	0.657 *	1

* represents significant correlation; ** represents highly significant correlation.

). Relationships between other services are not significant.

3.5. Effects of Desertification on Ecosystem Services

With the increase of DDI, the SC, WR, and NPP also increase significantly (Figure 8a–c). This suggest that all these ecosystem services improve when desertification decreases. CP and LP represented by non-significant unimodal curve with increasing DDI. CP and LP reached a maximum when DDI was approximately 131.32 (Figure 8d,e).

4. Discussion and Conclusions

4.1. Desertification Dynamics in 1986–2015

Numerous extensive studies on desertification pattern and dynamics in Mu Us Desert have been carried out, the results of these studies show that the main desertification in Mu Us Desert basically remained unchanged with the alternation of intensification and reversal, the eastern part of the desert is characterized by largely declining desertification. However, some individual localities in the western part of the region experienced significant desertification [28,29,30]. Our findings generally support those from previous studies of desertification patterns and dynamics in the Mu Us Desert. Desertification was not largely occurring in the region as a whole, despite some spatiotemporal heterogeneity and reversal trends, especially in the western part. Reversing trends were most obvious in Shaanxi Province. Desertification had been gradually decreasing since 1986, whereas no significant trends in desertification were observed from 2005 onward. The reversal and intensification of desertification in the Inner Mongolia Autonomous Region are both present since 1986. In Ningxia Autonomous Region, the first half of the study period is characterized by increasing desertification, which reached its peak in 2005. Desertification was then replaced by the opposite trend here (Figure 5).

Dynamics of desertification are driven by both geomorphic structure and climatic conditions [31], with contribution from human activities [32,33] and policies [18,34]. According to our results, climate change among the regional meteorological sites is highly consistent (Figure 4, Table 2 and Table 3), and if the climate is the dominant factor in the desertification change in the region, the trend of desertification among the three provinces/autonomous region should be the same, but the trend of desertification in the three provinces/autonomous region is quite different in this study, so climate factors in the Mu Us Desert are not likely to play a dominant role in the change of desertification patterns. We note a strong dependence of patterns of desertification dynamics on administrative divisions of this area. For example, in Shaanxi Province, where the Grain-for-Green Program (GFGP), which is the important component of development strategy of western frontiers of China, was introduced the earliest. It promptly yielded the desired outcome in the form of slowing desertification and ecosystem recovery here [35]. A series of ecological restoration projects have started, including the Double Rights and One System and Aerial Seeding Afforestation Project in Inner Mongolia Autonomous Region. Yet, the recovery of ecosystems did occur mostly because of different natural conditions and policy implementations. Thus, the coupling effect of anthropogenic factors and climate should explain most of the dynamics of desertification in this region during this study period at least. The driving factor of desertification trend may change from one period to another [36], but in this study, the driving factors have not been divided into a finer time scale, we therefore propose further in-depth analysis of drivers of desertification in the Mu Us Desert, which should be done in a follow-up study.

4.2. Relationships between Ecosystem Services

Relationships between different ecosystem services are often characterized as tradeoffs or synergies [37]. Previous studies have shown that synergy is often attributable to regulating, as well as provisioning services. The relationship between provisioning and regulating services has often the form of trade-offs as well [38,39,40]. Regulating services in this study include soil conservation, water retention, and net primary productivity (NPP). Provisioning services we studied include crop productivity and livestock productivity. The synergies between regulating services, such as soil conservation (SC), water retention (WR) and NPP, and between provisioning services, such as crop productivity (CP) and livestock productivity (LP) are strong, which is partly supported by our results (Table 4). Trade-offs between regulating services and provisioning services were not found by our study. So was the synergy trend.

SC, WR, and NPP are synergistically related to each other through vegetation cover. Higher NPP corresponds to dense vegetation cover [41], which in turn can reduce the wind speed while dense root systems of vegetation can enhance the ability of soil to resist erosion [42]. WR primarily depends on soil moisture content [8], and high NPP corresponds to better precipitation conditions. Plants in this region have high water use efficiency indicating that the regional soil has higher water holding capacity high [43,44]. Therefore, there is also a synergistic relationship between WR and NPP. In the region, the development of agriculture and livestock husbandry often depend on the quality of natural resources [45]. Higher NPP generally reflects better resource conditions. Consequently, the two provisioning services (CP and LP) are synergistic with three regulating services (SC, WR and NPP).

4.3. Impacts of Desertification on Ecosystem Services

Negative correlation between ecosystem services and desertification has been previously identified. Regulating services tend to decrease as desertification increases, while the relationship between desertification and other ecosystem services could exhibit non-linear patterns [17]. Our results provide additional support to those findings. Three regulating services examined in this study were positively correlated with DDI (Figure 8a–c) confirming the decline in those services with the intensification of desertification. The two provisioning services and DDI showed a hump-shaped relationship. Despite the lack of statistical significance of this relationship, it is interesting to note the potential of a decreasing trend in provisioning services after desertification passes some threshold and does not affect those services. Perhaps some other, unaccounted for in the study, factors could explain the decline of services in lieu of desertification (Figure 8d,e).

The explanation of the decline of three regulating services with the intensification of desertification is straightforward. With respect to desertification, the most basic functions of the ecosystem are soil and water conservation functions [17]. In the Mu Us Desert, the loss of vegetation cover triggers the loss of soil organic matter. This reduction has direct effects on soil structure by weakening soil aggregates. Such a condition accelerates the disaggregation of those particles by raindrops and wind erosion. Land cultivation and heavy grazing can further impact soil conditions. The combined effects of these soil erosion processes can be summarized as follows. First, in sandy substrates soil nutrients, such as nitrogen, phosphorus, and potassium, are eroded fast and soil fertility decreases [14,15,46]. Second, eroded soil becomes porous and permeable and cause the water regulation function to degrade and create a negative effect for biomass production and re-vegetation [17]. Third, sandy soils, when eroded, lead to increased contamination of atmosphere by dust thus interfering with vegetation photosynthesis, respiration, and transpiration, which all results in declining primary production [47,48,49]. In summary, desertification directly reduces three regulating services examined in this study. With the intensification of desertification, the two provisioning services respond non-linearly by first increasing and then decreasing. This may be an indication that human activities play a prominent role here leading to complex relationships between provisioning services and desertification [50]. In the initial stages of desertification or recovery from it, increases in farming and grazing take advantage of increased provisioning services and exacerbate soil erosion due to human disturbance. This is consistent with the first half of the hump-shaped curve obtained in this study. As human activities continue to intensify alongside with increasing desertification, the carrying capacity threshold is passed, which disrupts water and soil conservation functions. Sustained crop growth and livestock feed cannot be maintained, and provisioning services decline as desertification increases. This process corresponds to the latter half of the hump-shaped curve. Therefore, when the two parts are combined, the overall relationship turns into the hump-shaped curve found in this study.