*Article* **Enhancing Ecosystem Services in the Agro-Pastoral Transitional Zone Based on Local Sustainable Management: Insights from Duolun County in Northern China**

**Luwei Dai 1,2, Haiping Tang 1,2,\*, Yunlong Pan 1,2 and Dalin Liang 1,2**


**Abstract:** Ecosystem and associated ecosystem services (ESs) in the agro-pastoral transitional zone of northern China (APTZNC) are sensitive to climate change and human activities. Essential to designing targeted policy interventions toward achieving sustainability in the APTZNC is a comprehensive understanding of the spatiotemporal changes in ESs and their drivers. This study identified the spatiotemporal changes in six ESs in Duolun County from 2000 to 2017. The impacts of drivers temperature, precipitation, wind speed, vegetation cover (FVC), land use/cover (LULC), soil type, altitude, and slope—on the changes in the ESs in the county and its ecological production zones were then explored. The results indicated that the six ESs improved during the study period. The drivers influencing changes in ESs over time exhibited similarities across regions. Although FVC contributed to improvements in the food supply, grass production, carbon sequestration, and soil wind erosion (SLwind), it also reduced water yield, which may exacerbate the water shortage in arid and semi-arid areas. In regions where the ecology was in the recovery phase, especially in slope farmland, the inhibition of soil water erosion (SLwater) by FVC was easily offset by the higher SLwater potential from increased precipitation. The decrease in wind speed improved the regional ESs, whereas the increase in temperature posed a threat to SLwind. The drivers affecting the spatial patterns of ESs varied among zones. Across the three zones, the greater influential drivers of ESs were FVC and LULC. The impacts of topographic drivers and soil type on the distribution of ESs should also be noted in the agro-zone and agro-pastoral zone, respectively. Our study advocated that ES management should be adjusted to local conditions, and differentiated planning policies should be implemented in line with the ecological characteristics in the APTZNC, which will contribute to regional ecological sustainable development.

**Keywords:** ecosystem services; spatiotemporal changes; driving factors; agro-pastoral transitional zone; management

#### **1. Introduction**

The world is facing various crises, such as climate change and loss of biodiversity, which seriously threaten the survival and development of human beings [1,2]. In response to global challenges such as climate change and environmental degradation, the United Nations General Assembly declared 2021–2030 to be the Decade on Ecosystem restoration, dedicated to promoting and restoring ecosystem services (ESs) [3]. Nature-based solutions (NbS) also refer to addressing a wide range of human challenges by protecting and sustainably using the vital ESs provided by natural ecosystems [4,5]. However, rapid population growth and the overconsumption of natural resources have led to a deterioration in the capacity of ecosystems to provide ESs, and sustainable development is under serious

**Citation:** Dai, L.; Tang, H.; Pan, Y.; Liang, D. Enhancing Ecosystem Services in the Agro-Pastoral Transitional Zone Based on Local Sustainable Management: Insights from Duolun County in Northern China. *Land* **2022**, *11*, 805. https:// doi.org/10.3390/land11060805

Academic Editor: Manuel Pulido Fernádez

Received: 9 May 2022 Accepted: 26 May 2022 Published: 28 May 2022

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

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

threat [6–8]. Therefore, it is imperative to monitor and evaluate the spatial and temporal changes in ESs and explore the driving mechanisms behind the observed changes to help policymakers formulate effective policies for managing ESs [9–12].

There is abundant evidence that the generation of ESs depends on socio-biophysical factors and ecological processes, and that ESs have scale dependence in space and time [1,13,14]. It is necessary to quantify and map the ESs at broad temporal and spatial scales to help determine restoration priorities and sustainable development management of ecosystem services [1,15,16]. Multiple methods have been developed and deployed to quantify ESs [17–21]. With advances in modeling techniques, the use of ecological models to assess and spatially present ESs can better relate findings to ecosystem structures and functions, which can also provide more quantitative evidence for decision-making [8,20,21]. Despite rapid progress in quantifying and mapping ESs, there are still some fundamental issues that have not been adequately addressed. One of these issues is that existing studies have mostly focused on quantifying ESs at a single point in time or over several time intervals [22,23]. The results of these studies did not allow for general conclusions and are highly likely to mask some uncertainties due to extreme environmental disturbances or human activities during the year [5,24], especially in ecosystems where the landscape is characterized by rapid change and heterogeneity. To successfully manage natural resources and the related ESs, research aimed at long time series perspectives may provide deeper insights into ES changes and underlying ecological and anthropogenic drivers than time-point analyses [25,26].

The effectiveness of ES management is affected by many biophysical and anthropogenic drivers, such as terrain diversity, soil type, climate change, and land use/ cover [6,13,24,27–29]. Different drivers affect the supply of ESs in time and space in different ways. For example, climate change affects the supply of ESs by influencing biomes through temperature and precipitation [28]. Natural background conditions (e.g., topography, soils, and other drivers) directly influence regional ecosystem structures, resulting in the spatial heterogeneity of ESs [30]. Whether the impact of these drivers is positive or negative depends on the ES, and the magnitude of these impacts varies from place to place [2,16,31,32]. Identifying the impacts of the drivers on ESs on a time scale can help government departments optimize local response strategies in the context of climate change and intense human activity [6,33,34]. Furthermore, identifying the drivers of spatial heterogeneity of the ESs has facilitated the implementation of regional ecosystem planning and decision-making [23,35]. However, most studies have only analyzed the drivers of the temporal changes or spatial changes in ESs [28,32,33], and the results could not provide more detailed information for the formulation of ecological restoration policies. More efforts are needed to deeply explore the spatial and temporal changes and driving mechanisms of ESs in combination with time-series data to manage and improve multiple ESs by manipulating the drivers.

We focused on the agro-pastoral transitional zone of northern China (APTZNC), which divides agricultural areas from pastoral areas [36] and is situated in the transition from semiarid to arid regions [37]. APTZNC includes highly diverse ecosystems (e.g., farmland ecosystems, grassland ecosystems, and woodland ecosystems) and supports the provision of rich ESs (e.g., food supply, grass production, carbon sequestration, and soil and water conservation) [38]. To date, studies of the APTZNC have provided valuable information on ESs at large spatial scales, such as cities [39], urban agglomerations [36], and even the entire APTZNC [38,40]. However, studies that have identified drivers of change in multiple ESs at local scales (e.g., county administrative districts) are scarce. There are significant differences in precipitation, temperature, biological communities, and human activities within the APTZNC [37], and ESs are prone to change at temporal and spatial scales [40]. It is difficult to apply the results of watershed or urban-scale studies to ecological restoration management at county-level administrative units (the most basic administrative unit), which thus may not be effective in improving the supply capacity of ESs [41]. Within the county administrative districts, combining the field research situation and ES assessment

results can provide a more in-depth and accurate analysis of the ecological status of the region and detailed information for site-specific ecosystem management. On the other hand, analyzing the supply and drivers of ESs in different ecological production zones (e.g., the agricultural advantage zone, the pastoral advantage zone, the semiagricultural and semipastoral zone) is essential to balance the regional ecological protection and economic development. Therefore, it is necessary to carry out county-scale research regarding ESs and their drivers in the APTZNC to serve local policy decisions.

Combined with many years of field research, we focused on Duolun County, a typical county administrative district located in the APTZNC. The ultimate goal of our study was to identify the drivers of spatial and temporal changes in ESs at the local scale based on timeseries data. Specifically, we aimed to understand the following: (1) What are the changing trends in ESs from 2000 to 2017?; (2) How do drivers influence the supply of ESs over time?; and (3) What are the main drivers influencing the spatial distribution of ESs in county administrative districts and different ecological production zones? Our intent was that our results could contribute to an improved understanding of the key anthropogenic and biophysical processes underlying the supply of the ESs in a county administrative district of the APTZNC and provide scientific support for promoting sustainable development management in each ecological production zone.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Duolun County (41◦46 –42◦36 N, 115◦51 –116◦54 E), which covers 3 towns, 2 townships, and 65 administrative villages, is located in the middle of the APTZNC and covers a total area of 3.95 × <sup>10</sup><sup>3</sup> km<sup>2</sup> (Figure 1b). The topography of the study area is semi-annular basin-like, with an elevation range from 1149 to 1796 m (Figure 1c). The soil types can be divided into 7 soil types and 14 subtypes, and the main soil types are chestnut soil, aeolian sand soil, and meadow soil. Duolun County is in a typical agro-pastoral transitional zone. In August 2017, we visited Duolun County for field investigation, and reviewed and verified the vegetation types, land use types and soil types at the sampling sites (Figure 1c). Through investigating the ecological environment, we confirmed that the actual soil types and the land use types of the sampling sites were roughly the same as the soil data used in this study and the land use/cover data of 2015. We combined the administrative boundaries of townships and regional ecological production patterns to divide Duolun County into three regions (Figure 1c). There are large grasslands in the north and east of Duolun County, which are dominated by animal husbandry (Figure 1(cII)); the dominant plants in this grassland mainly include *Leymus chinensis*, Agropyron cristatum, Stipa krylovii, *Cleistogenes squarrosa*, *Congsheng grass* and *Artemisia frigida*. The south of the study area is an important grain-producing area; farmland mainly planting annual flax, oats, buckwheat, spring wheat and also some silage corn (Figure 1(cI)), The central part of the study area is an ecotone between agriculture and animal husbandry (Figure 1(cIII)), with a similar proportion of the development of agriculture and animal husbandry. In the 21st century, a series of ecological projects have been implemented in this area, such as the Grain for Green Project, grazing prohibition projects and Beijing–Tianjin sandstorm control engineering, to mitigate environmental pressure [41,42]. The field research found that there was still land desertification around the county town and in the northern area.

**Figure 1.** Locations in the study area and general description of geographical information: (**a**) the location of the APTZNC in China; (**b**) the location of Duolun County in the APTZNC; (**c**) a digital elevation model (DEM) of the study area; (**cI**) survey photographs of the agricultural advantage zone (agro-zone); (**cII**) survey photographs of the pastoral advantage zone (pastoral zone); and (**cIII**) survey photos of the semiagricultural and semipastoral zone (agro-pastoral zone).

Duolun County has a temperate continental arid climate with annual precipitation of 344–399 mm and a range of mean annual temperatures of 3.1 to 4.5 ◦C. From 2000 to 2017, the climate conditions and vegetation cover in Duolun County underwent obvious changes (Figure 2). The annual precipitation and annual mean temperature increased insignificantly at a rate of 4.526 mm yr−<sup>1</sup> and 0.013 C yr−1, respectively, whereas the annual mean wind speed decreased insignificantly at a rate of 0.008 m s−<sup>1</sup> yr<sup>−</sup>1. This indicated that the climate in the study area is becoming warm and humid. In addition, the annual vegetation cover increased significantly at a rate of 0.614% yr−1. The vegetation cover conditions in the study area have been improving.

**Figure 2.** Temporal variation in the annual (**a**) precipitation, (**b**) mean temperature, (**c**) mean wind speed and (**d**) vegetation cover in Duolun County during 2000–2017. The gray areas represent the 95% confidence intervals.

#### *2.2. Methodological Framework and Data Sources*

After the field investigation in 2017, we selected and evaluated six key ESs of high relevance to stakeholders in the region for assessment and analysis, including food supply (FS), grass production (GP), water yield (WY), carbon sequestration (CS), soil water erosion (SLwater) and soil wind erosion (SLwind). Then, combined with the time-series data, the spatial and temporal changes in the ESs were analyzed in depth. Finally, the spatiotemporal change drivers of the ESs from the perspective of the zones were analyzed. We selected four major categories of eight drivers to research the causes of the spatial and temporal changes in the ESs. The selected drivers were soil type (ST), meteorological drivers (annual precipitation (Pre), annual mean temperature (Tem), and annual mean wind speed (WS), whereas topographic drivers were altitude (Alt), slope (Slo)), and the anthropogenic drivers were vegetation cover (FVC), land use/cover (LULC). Table 1 shows the datasets used in this study. All raster maps were converted to the UTM coordinate system at a spatial resolution of 100 m.


**Table 1.** Descriptions and sources of the study data.

#### *2.3. Mapping Ecosystem Services*

#### 2.3.1. Food Supply

According to the positive NDVI, the yields of crops such as grains, oilseeds, and vegetables, as well as the yields of meat and milk in the statistical yearbook of the corresponding year, were allocated to the cultivated land and grassland grids, respectively. The formula can be expressed as follows [29]:

$$P\_i = \frac{NDVI\_i}{NDVI\_{sum}} \times P\_{sum} \tag{1}$$

where *Pi* represents the crop yield or meat and milk yield (kg/hm2) of grid *i*, *Psum* is the crop yield or meat and milk yield (kg/hm2) in the county, *NDV Ii* is the NDVI value of farmland grid *i* or grassland grid *i*, and *NDV Isum* is the sum NDVI of the farmland or grassland in the county.

#### 2.3.2. Grass Production

In this study, an inversion model of regional grass yield was used to estimate the grass yield per unit area. The calculation formula is as follows [43]:

$$AGB\_{\rm i} = 1741.089NDVI\_{\rm i}^2 + 2130.557NDVI\_{\rm i} - 424.757$$

where *AGBi* represents the grass yield (kg/hm2) of grid *i* and *NDV Ii* is the NDVI value of grassland grid *i*.

#### 2.3.3. Water Yield

In this study, water yield was simulated using the Integrated Valuation of ESs and Trade-offs model (InVEST, version 3.3.0), which was developed by Stanford University, The Nature Conservancy (TNC) and the World Wildlife Fund (WWF). The model is expressed as follows [17]:

$$Y\_x = \left(1 - \frac{AET\_x}{P\_x}\right) \times P\_x \tag{3}$$

$$\frac{AET\_x}{P\_x} = 1 + \frac{PET\_x}{P\_x} - \left[1 + \left(\frac{PET\_{xj}}{P\_x}\right)^{w\_x}\right]^{1/w\_x} \tag{4}$$

$$w\_{\mathbf{x}} = Z \frac{AW\mathbf{C}\_{\mathbf{x}}}{P\_{\mathbf{x}}} + 1.25\tag{5}$$

where *Yx* is the annual water yield (mm) for grid point *x*, *AETxj* is the annual actual evapotranspiration (mm) for pixel *x*, *Px* is the annual precipitation (mm) in pixel *x*, *Wx* is a nonphysical parameter that characterizes the natural climatic-soil properties, *PETx* is the potential evapotranspiration (mm), *AWCx* is the volumetric plant-available water content (mm), and *Z* is an empirical constant which is used to characterize the seasonal distribution of precipitation [41].

#### 2.3.4. Carbon Sequestration

In this study, NPP was used as a proxy for carbon sequestration and estimated using the Carnegie–Ames–Stanford Approach (CASA) model [44]. The model is expressed as follows:

$$NPP(\mathbf{x},t) = APAR(\mathbf{x},t) \times \varepsilon(\mathbf{x},t) \tag{6}$$

$$APAR(\mathbf{x}, t) = SOL(\mathbf{x}, t) \times FPAR(\mathbf{x}, t) \times 0.5 \tag{7}$$

$$
\varepsilon(\mathbf{x},t) = \varepsilon\_{\max} \times T\_{\varepsilon1}(\mathbf{x},t) \times T\_{\varepsilon2}(\mathbf{x},t) \times \mathcal{W}\_{\varepsilon}(\mathbf{x},t) \tag{8}
$$

where *NPP*(*x,t*) is the net primary productivity of pixel *x* in month *t* (gC/m2), *APAR*(*x,t*) is the photosynthetically active radiation absorbed by pixel *x* in month *t* (MJ/m2), *ε*(*x,t*) is the light utility efficiency of pixel *x* in month *t* (gC/MJ), *SOL*(*x,t*) is the total solar radiation on pixel *x* in month *t* (MJ/m2), and *FPAR*(*x,t*) is the ratio of the absorption of the incoming photosynthetically active radiation by the vegetation layer (dimensionless) [41]. The constant 0.5 reflects the proportion at which the effective solar radiation accounted for the total solar radiation, *εmax* is the maximum light use efficiency of the vegetation (gC/MJ), whereas *Tε*1(*x,t*), *Tε*2(*x,t*), and *Wε*(*x,t*) refer to parameters describing the stress coefficients at the highest temperature, the stress coefficients at the lowest temperature and the water stress coefficient in cell x in month t, respectively [45].

#### 2.3.5. Soil Water Erosion

The soil water erosion was calculated using the revised universal soil loss equation (RUSLE) [46]. The formula is expressed as follows:

$$Lwater = R \times K \times LS \times \mathbb{C} \times P \tag{9}$$

where SLwater represents soil erosion water (t/(hm2)), *R* is the rainfall erosion factor (MJ·mm/(hm2·h)), *<sup>K</sup>* is the soil erosion index (t·h/(MJ·mm)), *<sup>C</sup>* is the vegetation cover index (dimensionless), *P* is the soil erosion control practice factor (dimensionless), and *LS* is the slope length and slope gradient factor (dimensionless). Detailed information on the parameter localization of RUSLE is presented in Supplementary Table S1.

#### 2.3.6. Soil Wind Erosion

The soil wind erosion was calculated using the revised wind erosion equation (RWEQ) model [47]. The basic equations are as follows:

$$SLwind = \frac{2z}{S^2} \times Q\_{\text{max}} \times e^{-\left(z/S\right)^2} \tag{10}$$

$$Q\_{\text{max}} = 109.8 \times \left( \text{WF} \times EF \times \text{SCF} \times K' \times \text{COG} \right) \tag{11}$$

$$S = 150.71 \times \left( \text{WF} \times EF \times \text{SCF} \times K' \times \text{COG} \right)^{-0.3711} \tag{12}$$

where SLwind is soil wind erosion under the conditions of the ground cover vegetation (kg/m2), *z* is the distance from the upwind edge of a field (m), *Qmax* is the maximum transport (kg/m), *S* is the critical field length (m), *K* is the surface roughness, *WF* is the climatic factor (kg/m), *EF* is the soil erodibility factor (%), *SCF* is the soil crust factor, and *COG* is the combined vegetation factor. Detailed information on the parameter localization of the RWEQ is presented in Supplementary Table S2.

#### *2.4. Trend Analysis of the Time Series Data*

The Theil–Sen median (Sen) estimation and Mann–Kendall nonparametric test [48,49] were used to detect the statistically significant changing trends and trend slopes in the long time series ESs at a significance level of 0.05. Sen's slope indicates the direction and amplitude of the variables' changes with time, whereas the positive and negative slopes indicate increasing and decreasing trends, respectively. The Mann–Kendall test does not require the samples to have normal distributions and is less sensitive to outliers, which is extensively employed to analyze long time series data [33]. The calculation was performed with MATLAB R2020b programming.

#### *2.5. Estimating the Impact of the Drivers on the Ecosystem Services*

#### 2.5.1. Drivers of the Temporal Changes in ESs

Correlation analysis and partial correlation analysis were used to detect the relationship between ESs and drivers (Pre, Tem, WS and FVC). When multiple drivers are related to ESs at the same time, the use of partial correlation analysis allows removal of the influence of the remaining drivers, and the relationship between a single driver and ESs is analyzed separately. The calculation was performed with R statistical software and OriginPro 2022 software [50]. In addition, the annual average values of each ES in individual land cover categories (only unchanged land cover) in Duolun County and different ecological production zones were calculated to compare the ES supply capacity of the major land cover (farmland, woodland, and grassland). The calculations were performed with ArcGIS 10.5 software and OriginPro 2022 software [51].

#### 2.5.2. Drivers of the Temporal Changes in ESs

A geographical detector is a statistical method used to detect the spatially stratified heterogeneity of geographic phenomena and reveal nonlinear associations between potential drivers and geographic phenomena [52]. A geographical detector assumes that if there is spatial consistency between independent variable X and dependent variable Y, then a statistical association is present between them. The advantages of this method are that it does not need a linear hypothesis, and its physical meaning is clear. It contains four formulas: a factor detector, an interaction detector, a risk detector and an ecological detector. In this study, the factor detector module was selected to evaluate the explanatory power of the independent variable X to the dependent variable Y. The independent variable X was assigned to the eight drivers, and the dependent variable Y was assigned to the seven types of ES supply changes. The formula to measure the *q* value is as follows:

$$q = 1 - \frac{\sum\_{h=1}^{L} N\_h \sigma\_h^2}{N\sigma^2} \tag{13}$$

where *q* is the explanatory power of variable X to the spatial variation in variable Y. The *q* value ranges from 0 to 1, and the value means that X explains *q* × 100% of Y. The larger (or smaller) the q-statistic is, the stronger (or weaker) the explanatory power of the independent variable to explain the dependent variable. L is the number of classifications or partitions of X; *N* and *Nh* are the numbers of units in the entire study area and subregion h, respectively. *σ*<sup>2</sup> and *σ*<sup>h</sup> <sup>2</sup> are the variances of variable Y over the entire study area and subregion h, respectively.

The input data of geographic detectors must be in the form of categorical layers (such as soil type and land use type); therefore, continuous datasets (such as precipitation and temperature) must be categorized. In this study, the precipitation, temperature, wind speed, vegetation coverage, elevation, and slope were divided into six strata by the natural break method.

#### **3. Results**

#### *3.1. Changing Trends in Ecosystem Services during 2000–2017*

#### 3.1.1. Temporal Changes in Ecosystem Services

During 2000–2015, the six ESs showed similar trends in annual mean values in the county administrative district, agro-zone, pastoral zone and agro-pastoral zone (Figure 3). The amount of FS exhibited an insignificant upward trend (*p* > 0.05), with the largest slope in the agro-zone (slope = 7.1 kg m−<sup>2</sup> yr<sup>−</sup>1) and the smallest slope in the agro-pastoral zone (slope = 5.025 kg m−<sup>2</sup> yr<sup>−</sup>1). There was a trend of significant increase (*p* < 0.05) in GP both in the county administrative district and different ecological production zones. Notably, there are strong interannual fluctuations in CS and WY with a non-significant increasing trend (*p* > 0.05), with the largest increase rate in WY in the agro-zone (slope = 2.007 mm yr<sup>−</sup>1) and CS in the pastoral zone (slope = 1.623 gC m−<sup>2</sup> yr−1). Generally, SLwater and SLwind showed a downward trend, especially in the agro-pastoral zone, and the rate of decrease reached −0.3 t hm−<sup>2</sup> yr−<sup>1</sup> (*<sup>p</sup>* = 0.068) and −0.016 kg m−<sup>2</sup> yr−<sup>1</sup> (*<sup>p</sup>* < 0.05), respectively.

Combined with the results of the pixel-by-pixel trend analysis (Figure 4), we found that the trends in the six ESs were spatially different. Except for the two negative ESs (SLwater and SLwind), all four ESs showed a good trend in that the gain areas were larger than the loss areas. Interestingly, the decrease in WY (6.2%) was concentrated in the central part of this study area, where GP and CS increased significantly. The increase in CS was mostly located in the pastoral zone, whereas the losses were mainly located in the agro-zone and the agro-pastoral zone. SLwater and SLwind showed good trends in that the areas of improvement were much larger than the areas of their gain (Figure 4), with the areas of increases in SLwater scattered throughout the county.

**Figure 3.** *Cont.*

**Figure 3.** The temporal dynamics of the six ESs during 2000–2017. We determined the average annual values of each ES and plotted their overall trends. The first column (**a**), second column (**b**), third column (**c**) and fourth column (**d**) show the trends in the six ESs in Duolun County, the agricultural advantage zone (agro-zone), the pastoral advantage zone (pastoral zone) and the semiagricultural and semipastoral zone (agro-pastoral zone), respectively. The first to sixth rows indicate the trends of food supply (FS), grass production (GP), water yield (WY), carbon sequestration (CS), soil water erosion (SLwater), and soil wind erosion (SLwind), respectively. The gray areas represent the 95% confidence intervals.

**Figure 4.** Spatial patterns of the six key ESs change trends and significance levels (*p* < 0.05) from 2000 to 2017. Zone I: the agricultural advantage zone (agro-zone), Zone II: the pastoral advantage zone (pastoral zone), Zone III: the semiagricultural and semipastoral zone (agro-pastoral zone). FS: food supply, GP: grass production, WY: water yield, CS: carbon sequestration, SLwater: soil water erosion, and SLwind: soil wind erosion.

#### 3.1.2. Spatial Changes in Ecosystem Services

The spatial patterns of the ESs in Duolun County exhibited regional heterogeneity (Figure 5), and the spatial distribution patterns of each ES remained stable (Figures S1–S6). Taking 2017 as an example (Figure 5), FS shows a decreasing trend from the agro-zone to the pastoral zone (Figure 5a). In contrast, GP and CS show an opposite pattern compared with FS, i.e., the high supply was concentrated in the northern and eastern areas of the pastoral zone (Figure 5b,d). The spatial distribution of WY patterns changed slightly in different years, but most of the high-value zones were distributed in the agro-pastoral zone, with a few in the northern and eastern parts of the pastoral zone (Figure 5c). SLwater showed a decreasing characteristic from southwest to northeast, with most of the higher areas concentrated in the agro-zone and the eastern part of the pastoral zone, whereas SLwater was less concentrated in the agro-pastoral zone (Figure 5e). Although the SLwind has decreased significantly over the past 18 years (Figure S6), the relatively higher areas can still be found mainly in the northern pastoral zone of Duolun County, and the lower areas are more stably distributed in the southern and eastern regions (Figure 5f).

**Figure 5.** Spatial patterns of the six ESs in 2017. Zone I: the agricultural advantage zone (agro-zone). Zone II: the pastoral advantage zone (pastoral zone). Zone III: the semiagricultural and semipastoral zone (agro-pastoral zone). (**a**) FS: food supply (in kg/hm2), (**b**) GP: grass production (in kg/hm2), (**c**) WY: water yield (in mm), (**d**) CS: carbon sequestration (in gC/m2), (**e**) SLwater: soil water erosion (in t/hm2), and (**f**) SLwind: soil wind erosion (in kg/m2).

#### *3.2. Drivers of Temporal Changes in Ecosystem Services*

Figure 6 represents the correlation coefficients based on the time series data between ESs and drivers (including Pre, Tem, WS, and precipitation) for Duolun County and different regions. The results of the study showed that the important factors affecting changes in ESs over time were essentially the same in Duolun County and different zones. The results showed that (1) Pre, Tem and FVC exhibited a positive correlation with FS and GP, whereas WS exhibited a negative correlation with FS and GP; (2) FVC and Pre were the two drivers that positively affected the changes in WY, whereas Tem and WS showed insignificant inhibitory effects on WY; (3) Pre, Tem and FVC had a positive effect on CS, whereas WS had an inhibitory effect on CS; (4) Pre, WS and FVC mainly showed an insignificant promotion effect on SLwater, whereas Tem showed an inhibitory effect on SLwater; (5) FVC and Pre had suppressive effects on SLwind, whereas Tem and WS had facilitating effects on SLwind.

**Figure 6.** Correlation analysis between the ESs and the four drivers (vegetation cover (FVC), precipitation (Pre), temperature (Tem), and wind speed (WS)) in different regions from 2000 to 2017. The blue and red circles above the diagonal indicate positive and negative correlations, respectively. The asterisks in the circles show the significance degree (∗ for *p* < 0.05). The numbers below the diagonal indicate Spearman's correlation coefficients, with their color matching those of the corresponding circles. Darker colors demonstrate stronger correlations. FS: food supply, GP: grass production, WY: water yield, CS: carbon sequestration, SLwater: soil water erosion, SLwind: soil wind erosion. agro-zone: the agricultural advantage zone, pastoral zone: the pastoral advantage zone, agro-pastoral zone: the semiagricultural and semipastoral zone.

Combined with the results of partial correlation analysis (Figure 7): (1) FVC was significantly and positively correlated with FS and GP, whereas WS showed a negative correlation with FS and GP; (2) the correlation between Pre and WY was still significant and positive, but the correlation between FVC and WY became significantly negative from significantly positive; (3) the relationship between CS and drivers did not change, but the driver with the greatest correlation between CS in agro-zone, pastoral zone, and agro-pastoral zone was different, namely, Pre, WS, and WS; (4) FVC and SLwater changed from a positive correlation to a negative correlation, whereas Pre and SLwater were still positively correlated; (5) Tem and WS were significantly and positively correlated with SLwind, whereas Pre and FVC were still negatively correlated with soil wind erosion.

**Figure 7.** Partial correlation coefficients (∗ for *p* < 0.05) between the ESs and the four drivers (vegetation cover (FVC), precipitation (Pre), temperature (Tem), and wind speed (WS)) in different regions from 2000 to 2017. (**a**) Partial correlation coefficients in Duolun County, (**b**) partial correlation coefficients in the agricultural advantage zone (agro-zone), (**c**) partial correlation coefficients in the pastoral advantage zone (pastoral zone), (**d**) partial correlation coefficients in the semiagricultural and semipastoral zone (agro-pastoral zone). FS: food supply, GP: grass production, WY: water yield, CS: carbon sequestration, SLwater: soil water erosion, SLwind: soil wind erosion.

From 2000 to 2017, the sum of farmland, grassland and woodland in Duolun County exceeded 80% of the total area of the study area; thus, this paper analyzed the changes in the annual average ES supplies of the three key land use/cover types (only unchanged land cover) over the past 18 years (Figure 8). The results indicate that ES supplies differ by land-use type. Specifically: (1) The FS of farmland fluctuated widely from year to year, with a multiyear average of approximately 491.35 kg/hm<sup>2</sup> in Duolun County. In contrast, the FS of grassland was more stable, with a multiyear average of approximately 69.1 kg/hm2; (2) All the GP in this study was provided by grassland. The multiyear average value of GP in Duolun County was 1633.4 kg/hm2 and the highest GP (1672 kg/hm2) was in the pastoral zone; (3) There were significant differences in the effects of land-use types on the WY, specifically farmland > grassland > woodland. The differences in the supply of WY by similar land-use types in different zones were not significant; (4) The interannual fluctuations of the CS in the three land-use types were smaller and show the characteristics of woodland > grassland > farmland, with multiyear average values of 502.6 gC/m2, 392.6 gC/m2, and 347.5 gC/m2 in Duolun County, respectively. The supply of CS on woodland in the agro-zone (415.62 gC/m2) was significantly lower than in other zones; (5) 'The characteristics of SLwater in the three land-use types are the same as those of CS, i.e., the interannual fluctuations were smaller and woodland > grassland > farmland. The highest SLwater was found in woodland in the agro-zone. (6) In terms of the multiyear average values in the county, the highest SLwind erosion (0.1 kg/m2) was found in grassland, followed by farmland and woodland. Among the three zones, the highest SLwind was found in the agro-pastoral zone, followed by the pastoral zone.

**Figure 8.** *Cont.*

**Figure 8.** The supply of the six ESs in a different land-use cover (farmland, grassland, and woodland) from 2000 to 2017. The first column (**a**), second column (**b**), third column (**c**) and fourth column (**d**) indicate the distribution of annual mean values of ESs on different land use cover in Duolun County, the agricultural advantage zone (agro-zone), the pastoral advantage zone (pastoral zone) and the semiagricultural and semipastoral zone (agro-pastoral zone), respectively. The first to sixth rows indicate the distribution of annual mean values of food supply (FS), grass production (GP), water yield (WY), carbon sequestration (CS), soil water erosion (SLwater), and soil wind erosion (SLwind)

#### *3.3. Drivers of Spatial Changes in Ecosystem Services*

across different land use cover, respectively.

The dominant drivers affecting the spatial heterogeneity of the six ESs and their qvalues differed significantly within the county administrative district (Figure 9a). The top three drivers affecting the distributions of FS and GP were the same (i.e., FVC > ST > Tem), with FVC having the largest explanatory power for FS and GP, indicating that the distributions of the two ESs were more influenced by FVC. The largest explanatory power for WY was LULC (q = 0.73), and the following drivers were FVC, WS, ST, Alt, Pre, Tem, and Slo, in order of q. The important drivers for the distribution of CS were FVC and LULC, which explained more than 49.77% and more than 23.2% of the distribution, respectively, whereas the other drivers had weaker explanatory power. The slope was the dominant driver determining the distribution of SLwater (q = 0.36), followed by ST (q = 0.136) and Alt (q = 0.129), indicating that soil type and topographic drivers are important environmental drivers affecting SLwater. The FVC had the greatest effect on SLwind (q = 0.15), followed by ST (e.g., sandy soil distribution areas with high SLwind). The explanatory power of meteorological drivers, LULC, and topographic drivers on SLwind did not exceed 5%. In a comprehensive view, the explanatory powers of the four major categories of drivers on the distribution of SLwater are ranked as follows: topography drivers > soil type > anthropogenic drivers > meteorological drivers. The explanatory powers of the four major drivers on the spatial distribution of the other five ESs are as follows: anthropogenic drivers > soil type > meteorological drivers > topographic drivers.

The dominant drivers for the distributions of ESs in the three zones remained consistent with those of the county, and the importance of the other drivers changed with the ecological production zone. FVC was the dominant driver in the distributions of FS and GP in each ecological production zone. Furthermore, in the agro-zone, Slo (q = 0.136) and Pre (q = 0.104) contributed more to the distributions of FS, and Pre (q = 0.221) and Alt (q = 0.175) contributed more to the distributions of GP. In the pastoral zone, Tem (q = 0.097) and Alt (q = 0.094) had greater explanatory powers for FS, and Tem (q = 0.091) and WS (q = 0.071) had greater explanatory power for GP. The Slo had a greater explanatory power for the distributions of FS and GP in the agro-pastoral zone. LULC was still the dominant driver of WY in different ecological production zones, and the other important influencing drivers were meteorological drivers. FVC and LULC were the two most important drivers affecting the distributions of CS in all three ecological production zones, indicating that anthropogenic drivers had an important role in the distribution of CS. Comparing the three ecological production zones, the explanatory power of the slope on the distribution of SLwater exceeded 30% in all three zones. The explanatory power of Alt for SLwater was higher in the agro-zone and the agro-pastoral zone, 26.7%, and 34.1%, respectively, but only 1.5% for SLwater in the pastoral zone. Different from SLwater, SLwind in the different ecological production zones was influenced less by topographic drivers. The dominant drivers of SLwind in the agro-zone and pastoral zone were FVC (q = 0.388) and FVC (q = 0.21), respectively, whereas the dominant driver of SLwind in the agro-pastoral zone was ST (q = 0.246), followed by FVC (q = 0.175). In addition, ST and WS had greater explanatory power for SLwind in the pastoral zone and agro-pastoral zone, whereas Pre and Tem had greater explanatory power in the agro-zone

**Figure 9.** The q values of drivers affect the spatial distributions of ESs in Duolun County and different zones. The drivers include precipitation (Pre), temperature (Tem), wind speed (WS), soil type (ST), altitude (Alt), slope (Slo), vegetation cover (FVC), and land use/cover (LULC). (**a**) The q values of drivers affecting the spatial distributions of the ESs (FS, GP, WY, CS, SLwater, and SLwind) in Duolun County. (**b**) The q values of drivers affecting the spatial distributions of the ESs (FS, GP, WY, CS, SLwater, and SLwind) in the agricultural advantage zone (agro-zone). (**c**) The q values of drivers affecting the spatial distributions of the ESs (FS, GP, WY, CS, SLwater, and SLwind) in the pastoral advantage zone (pastoral zone). (**d**) The q values of drivers affecting the spatial distributions of the ESs (FS, GP, WY, CS, SLwater, and SLwind) in the semiagricultural and semipastoral zone (agropastoral zone). FS: food supply, GP: grass production, WY: water yield, CS: carbon sequestration, SLwater: soil water erosion, SLwind: soil wind erosion.

#### **4. Discussion**

*4.1. The Impacts of Anthropogenic and Meteorological Drivers on the Temporal Variations in the ESs*

As a result of the Grain for Green Project and the Beijing–Tianjin Sand Source Control Project (initiated in 2000), vegetation coverage has increased in Duolun County [42], effectively increasing grass production. Although there was a large amount of farmland conversion in the implementation of the ecological project [53], the food supply in Duolun County gradually increased through the gradual development of unused land and the construction of agricultural mechanization demonstration parks. In addition, our results showed that the increase in FVC reduced the probability of SLwater and SLwind (Figure 10). On the one hand, precipitation intercepted by the vegetation canopy can reduce the direct erosion of soil by rainfall, which, in turn, reduces the probability of SLwind [33]. On the other hand, the vegetation canopy can reduce wind speed, and the vegetation root system has the function of consolidating soil, which effectively enhances the soil resistance to erosion [32,54]. However, it is worth noting that WY is negatively influenced by FVC. Studies on the Loess Plateau [55] also pointed out that extensive vegetation restoration threatens the water supply required for human survival and vegetation growth. The rapid growth of vegetation in Duolun County since 2000 has consumed a large amount of soil water and increased evapotranspiration, thus significantly reducing the WY [32,34]. In particular, the study area is located in an arid and semi-arid region, and the actual precipitation and available water resources in recent years need to be taken into consideration when carrying out vegetation restoration projects to minimize the pressure on water resources caused by vegetation restoration.

**Figure 10.** Visualization of the correlation between temporal changes in ecosystem services and drivers (precipitation, temperature, wind speed, vegetation cover) from 2000 to 2017.

Although vegetation cover is considered to be the main factor in conserving soil from water erosion [54], this protection may be less than the higher SLwater potential due to increased precipitation in the ecologically fragile APTZNC. We noted that increased precipitation can significantly contribute to the occurrence of SLwater. Ecosystem structure is fragile within the APTZNC, and despite revegetation efforts on slope farmland, the damaged soil structure has not been fully restored due to the long history of farming [8,34], and soil erodibility is still high. Despite SLwater in the study area having improved, it is necessary to continue to strengthen the water and fertility retention capacity of the soil through ecological protection in the future, especially on sloping farmland. In addition to precipitation, changes in temperature and wind speed lead to changes in the ecological

environment and ES supplies. A warmer climate can increase evaporation from the surface, with drier soil surfaces and greater susceptibility to wind erosion [56]. As our results show, there is a significant positive correlation between temperature and SLwind. However, the mean annual temperature in Duolun County is low, and warmer temperatures can offset some of the wind erosion by promoting plant growth and increasing surface roughness [6,57]. In line with other studies [58], our results show that wind speed is an important driver of the occurrence of SLwind. The increase in wind speed not only tends to disperse the soil, but also accelerates the evaporation of soil water, which inhibits vegetation growth [33] and increases the risk of SLwater. The decrease in wind speed in the study area directly reduces the potential for SLwater and SLwind.

Land use/cover changes are also the greatest pressures affecting the provision of ESs [34,59]. Over the last 20 years, the growth of plantation forests in Duolun County has gradually entered semimature and mature stages, and the strong transpiration of the canopy consumes a large amount of water [34,54], resulting in a lower WY of woodland than that of farmland and grassland. In contrast, farmland has a higher capacity of WY due to less evaporation from vegetation, although its CS is lower than that of woodland and grassland [29,60]. Grassland has a lower water demand than woodland, and its CS capacity is greater than that of farmland. Increasing the area of grassland in arid and semiarid areas may be a compromise in terms of increasing CS and WY at the same time. SLwater and SLwind erosion were higher in the woodland and grassland in the study area than in the farmland, which is similar to our hypothesis. This is related to the spatial geography of this study area, and our numerous field studies have revealed that most plantation forests and grasses in the agro-zone have been planted in erodible landscapes characterized by relatively steep slopes and poor soils. Simultaneously, previous studies have shown that if landscape patches are disturbed, a patch may have difficulty blocking any erosive action [8]. Tree planting in the study area (especially in the agro-zone) is dispersed, leading to the fragmentation of woodland and grassland landscape, which, to some extent, weakens the soil and water conservation capacity of the two land use types. Although SLwater and SLwind were higher in woodland and grassland, the decreasing trends in SLwater and SLwind were higher in both land types than in farmland, which indicates that woodland and grassland are more capable of erosion control, especially woodland. Appropriately increasing woodland area and aggregation degree can improve the benefit of soil and water conservation.

#### *4.2. The Impacts of the Eight Drivers on the Spatial Changes in the ESs*

Ecological control measures and approaches can be explored in a targeted manner by clarifying the driving characteristics of ES spatial changes [61]. Our study shows that anthropogenic drivers (FVC and LULC) have a stronger influence on the spatial distributions of multiple ESs than other environmental factors in county administrative districts and different ecological production zones (Figure 11), indicating that the supplies of the ESs in the APTZNC depend to a large extent on the degree of restoration of ecological engineering and the rational allocation of land resources [62]. The regional heterogeneity of ecological production zones means that other drivers (meteorological drivers, topographic drivers, ST) affecting the spatial distributions of the six ESs show significant differences, which are closely related to the human and physical geography of the different regions. For example, topographic drivers (Alt and Slo) have a greater influence on the spatial distribution of CS in the agro-zone. However, in the pastoral zone and agro-pastoral zone, the influence of topographic drivers decreased, and the influence of the meteorological drivers increased. This pattern is formed because the spatial heterogeneity of Slo and Alt in the pastoral zone and agro-pastoral zone is small, and the water and heat conditions required for vegetation growth are almost entirely dependent on nature. In contrast, in the agro-zone, humans intervene and regulate moisture and temperature according to vegetation growth needs, and meteorological drivers such as natural precipitation have relatively little influence on the ecological environment. At the same time, the topographies

of the agro-zone are complex (high altitudes and steep slopes) and crops are grown with significant spatial heterogeneity, which can have a significant impact on the FVC and ESs [62]. Therefore, meteorological drivers have stronger explanatory power on the spatial distribution of CS in the pastoral zone and agro-pastoral zone, and topographic drivers have more influence on the spatial distributions of the ESs in the agro-zone. In addition, ST has a greater explanatory power than other drivers for SLwind in the agro-pastoral zone. Although soil types are abundant in Duolun County, they are mostly reflected in the agro-pastoral zone. As the foundation of terrestrial ecosystems, soil provides many important benefits and ESs to society, and changes in the physical properties of soils can affect ecosystem processes, and thus, the ESs [63,64]. Therefore, the spatial heterogeneity of soil type should be considered in regional ecosystem management to ensure a precise fit between management measures and soil type, thus contributing to the improvement of the various ESs.

**Figure 11.** Important drivers of the spatial distributions of ESs in different zones. The drivers include precipitation (Pre), temperature (Tem), wind speed (WS), soil type (ST), altitude, slope, vegetation cover (FVC), and land use/cover (LULC). FS: food supply, GP: grass production, WY: water yield, CS: carbon sequestration, SLwater: soil water erosion, SLwind: soil wind erosion.

Although the drivers of ES spatial heterogeneity change with the region, the effects of some environmental factors on ESs still show consistent characteristics, and this effectively reduces information redundancy needed for managing multiple ESs simultaneously to some extent. For example, LULC and meteorological drivers are always the dominant drivers affecting the distributions of WY, FVC is the dominant driver for FS, GP, CS and SLwind. More importantly, Slo has greater influences on SLwater than FVC. Soil erosion is more likely to occur in areas with steep slopes, and regional topographic features should be considered when developing management and soil erosion prevention measures.

#### *4.3. Implications for Integration of Ecosystem Services in APTZNC Management*

As a bridge between natural ecosystems and socio-economic systems, ESs are critical to conserving biodiversity, maintaining ecological security, and meeting human livelihood needs [37,38,40]. This study focused on the APTZNC, a transition zone where agricultural and pastoral areas are connected, which plays a critical role in sustaining the stability of natural ecosystems in northern China and provides an important guarantee for the livelihood of local and surrounding residents [37,41]. The sustainable management of ESs in the APTZNC is conducive to restoring the ecological environment and enhancing human well-being in the region [38,40]. Based on widely used biophysical models and remote sensing data, this study analyzed the drivers of spatiotemporal changes in ESs in the most basic administrative unit (i.e., county), which can provide some clues for ES management studies in other administrative units in the APTZNC. Our research suggests that as compared with the use of time node variations, the continuous change in ESs over a long period provides more meaningful results. Attempts to use discontinuous years as a study period to assess ESs for management purposes must be taken cautiously, and researchers should consider modeling ESs over continuous time periods and finer spatial scales [1,32]. In addition, ES management should be tailored to local conditions and zoning.

In the face of ecological engineering and climate change, ecosystems in the agro-zone still exhibit vulnerability. Although FS was highest in the agro-zone, CS in farmland, woodland, and grassland was the lowest of the three zones, and SLwater was significantly higher than in the other zones. Our field survey found that although the implementation of ecological projects has enhanced the restoration of vegetation on slopes, the damage caused by long-term tillage to the ecosystem has not been fully remediated. The study area should continue to strengthen the ecological protection of sloping land and reclassify areas to avoid the fragmentation of ecological woodlands and grasslands [8,65] to give full play to the water and fertility retention capacity of both land types. To promote the synergistic development of FS and the other ESs in the agro-zone, it is recommended to ensure the irrigated area of arable land through water conservancy engineering measures, to improve the water and fertility retention capacity of the soil through deep plowing and deep loosening techniques and to increase the application of organic fertilizers to maintain and enhance the ES supply capacities.

On the other hand, the supply of the six ESs in the agro-pastoral zone was at a low level in the whole county. There are bare soils and abandoned farmlands in the dryland regions of the agro-pastoral zone, which is the county seat, and it is necessary to actively promote vegetation restoration in unused land to improve the various ESs. It is important to pay attention to the characteristics of diversified soil types and restore vegetation in this area. At this time, attention needs to be paid to the recovery of sandy vegetation, which may result in a shortage of water resources in the local and surrounding areas [34]. To address this phenomenon, in addition to relying on natural precipitation, there is a need to reasonably exploit groundwater and take advantage of the county's large impervious layer to properly collect surface runoff as reserve water.

Finally, the distributions of the ESs in the non-intensively managed pastoral zone are susceptible to the influence of soil type and meteorological drivers, in addition to human activities. The northern part of the pastoral zone is distributed with large areas of grassland sandy soil, which still needs attention to prevent SLwind. We believe that the northern region should be well protected from the wind and that the expansion of farmland should be prohibited to avoid negative impacts on the other ESs due to deteriorating soil conditions [64]. In addition, we recommend that the northern part should increase the degree of grassland aggregation by planting some drought-tolerant or less waterconsuming grass species [33], and attention should be given to the negative impact of revegetation on terrestrial water storage (especially in arid and dry areas). The eastern part of the pastoral zone is rich in water and has larger areas of natural forests and grasslands, which are tourist destinations. The region should protect the existing vegetation, improve the quality and stability of the forest ecosystems, and avoid overexploitation.

#### *4.4. Limitations and Future Research Directions*

Some constraints in our analysis should be considered. First, due to the lack of annual high-resolution land use data, we had to use land use/cover data from one point in time to assume a gradual change in land use over a five-year period, which may have had some modest effect on the WY and some modest effect of the soil P-factor on SLwater. In the future, the use of consecutive years of land use/cover data for interannual WY and SLwater assessments deserves further study. Second, the correlation analysis only characterized the numerical relationships between the ESs and the drivers, and cannot characterize the threshold values of these relationships [54]. Combining analytical methods such as multiple regression and constraint lines to investigate the nonlinear relationships and interaction thresholds between driving factors and ecosystem services will be further explored in future research.

#### **5. Conclusions**

This study simulated spatiotemporal changes in six ESs over the last 18 years (2000–2017) and explored the drivers of changes in ESs of the county administrative district. This provides more detailed information for local ecosystem service studies and has the potential to assist decision-making processes. FS, GP, WY, and CS increased over time, and the capacity of the landscape to control water erosion and wind erosion was enhanced. The drivers of temporal changes in a single ES over time showed similarities across zones. The significant increase in FVC since 2000 has improved the FS, GP, CS, SLwater, and SLwind, but at the same time, has put pressure on water resources. Precipitation contributed to the improvement of WY and CS, but increased SLwater. The reduction in wind speed improved the six ESs, and the temperature had a significant promoting effect on SLwind. Exploring the drivers influencing the spatial distribution of the six ESs underscored the importance of anthropogenic drivers for the spatial distribution of ES over other environmental factors. Our findings also suggest the importance of integrating ES management with the ecosystem characteristics of each zone, because the influence of meteorological drivers, topographic drivers, and ST on the spatial distribution of ESs varied in different zones. For example, it is important in the agro-zone and agro-pastoral zone to account for Alt having a greater influence on SLwater, and in the agro-pastoral zone to account for ST being the main driver of SLwind. The study area should continue to strengthen the soil's ability to retain water and fertilizer in sloping fallow areas. The agro-pastoral zone needs to actively promote the revegetation of unused land and strive to improve the six ESs. The overall level of ESs in the pastoral zone is high, but attention is still needed to prevent and control SLwind.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/land11060805/s1, Table S1: Computational methods of parameters for RUSLE; Table S2: Computational methods of parameters for RWEQ; Figure S1: Spatial distribution of food supply from 2000 to 2017; Figure S2: Spatial distribution of grass production from 2000 to 2017; Figure S3: Spatial distribution of water yield from 2000 to 2017; Figure S4: Spatial distribution of carbon sequestration from 2000 to 2017; Figure S5: Spatial distribution water of soil erosion from 2000 to 2017; Figure S6: Spatial distribution of wind soil erosion from 2000 to 2017; Figure S7: Ecological background of the study area [56,66–72].

**Author Contributions:** Conceptualization, L.D. and H.T.; methodology, L.D., D.L. and H.T.; validation, L.D.; formal analysis, L.D. and D.L.; resources, H.T.; data curation, L.D.; writing—original draft preparation, L.D.; writing—review and editing, H.T.; visualization, L.D. and Y.P.; supervision, H.T.; project administration, H.T.; funding acquisition, H.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the National Natural Science Foundation of China (Grant No. 31972945).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

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

#### **References**


### *Article* **Spatial Distribution of Agricultural Eco-Efficiency and Agriculture High-Quality Development in China**

**Mingjia Chi 1, Qinyang Guo 2, Lincheng Mi 2, Guofeng Wang <sup>2</sup> and Weiming Song 1,\***


**Abstract:** Agricultural ecological efficiency is not only the key link between green development and high-quality development of agriculture, but also an important regulatory indicator for China's rural revitalization. Based on provincial panel data of China from 2000 to 2019, using land, mechanical, labor, fertilizer, pesticide, and agricultural film as input variables and economic output and agricultural carbon emissions as output variables, the inter-provincial agricultural ecological efficiency is calculated by a super-efficient SBM model, and the traditional spatial Markov probability transfer matrices are constructed based on time series and spatial correlation analyses. By exploring the spatial and temporal dynamic evolution characteristics of agricultural ecological efficiency, it is found that the agricultural ecological efficiency of China increased steadily with fluctuations. In addition, the provincial gap has been narrowing, but the overall level is still at a low level; thus, there is still a large space for improvement in agricultural ecological efficiency. The overall trend of agricultural ecological efficiency shifting to a high level in China is significant, but its evolution has the stability to maintain the original state, and achieving leapfrog transfer is relatively hard. The geospatial pattern plays an important role in the spatial-temporal evolution of agricultural ecological efficiency, with significant spatial agglomeration characteristics. Provinces with high agricultural ecological efficiency enjoy positive spillover effects, while provinces with low agricultural ecological efficiency have negative spillover effects; thus, gradually forming a "club convergence" phenomenon of "high agglomeration, low agglomeration, high radiation, and low suppression" in the spatial pattern. In addition, support for the improvement of agricultural ecological efficiency will be provided in this study.

**Keywords:** agricultural ecological efficiency; spatial-temporal evolution; SBM model; spatial Markov chain

#### **1. Introduction**

"If the nation is to be revitalized, the countryside must be revitalized", the explanation of General Secretary Xi Jinping about China's agriculture and rural areas shows the importance of agricultural and rural work to China's development [1]. Judging from the global development experience, agriculture is an important cornerstone of national security and stability [2]. Agricultural ecological efficiency (AEE) refers to the exchange of less natural resource consumption and environmental costs for more quantities and higher quality of agricultural products or agricultural services within the affordability of the ecosystem [3,4]. Since the reform and opening-up in 1978, based on the production model of petroleum agriculture, the agricultural output level in China has witnessed a continuous growth [5], but it has also paid a high price for resources and the environment [6]. Since 1978, China has made great achievements in agricultural development [7]. In 2021, the total grain output of China reached 683 million tons, increasing by 2% over the previous year. The annual grain output reached a new high level, maintaining more than 0.665 trillion kg for seven consecutive years. In 2015, the Ministry of Agriculture and Rural Affairs organized a zero-growth

**Citation:** Chi, M.; Guo, Q.; Mi, L.; Wang, G.; Song, W. Spatial Distribution of Agricultural Eco-Efficiency and Agriculture High-Quality Development in China. *Land* **2022**, *11*, 722. https://doi.org/ 10.3390/land11050722

Academic Editor: Hossein Azadi

Received: 27 April 2022 Accepted: 8 May 2022 Published: 11 May 2022

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

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

action for the use of chemical fertilizers and pesticides, and by the end of 2020, the reduction and efficiency improvement of chemical fertilizers and pesticides in China had successfully achieved the expected goals [8]. The use of chemical fertilizers and pesticides has been significantly reduced [9], while their utilization rate has been significantly improved, and the effect of promoting high-quality development of the planting industry was obvious [10]. In 2020, formula fertilizer accounted for more than 60% of the total application of the three major grain crops, the area of mechanical fertilization exceeded 46.66 million hectares, and the integration of water and fertilizer exceeded 9.33 million hectares. In addition, the green prevention and control area covered nearly 66.66 million hectares, and the green prevention and control rate of major crop diseases and insect pests was 41.5% in 2020, 18.5% higher than that in 2015. However, the green development of agriculture is still the main theme of the current agricultural development [11]. In the context of tightening resource and environmental constraints, the balance between economic and ecological benefits of agricultural development has become a major challenge facing China [12].

The green development of agriculture needs to clarify the relationship between agricultural input and output [13] and resource consumption and environmental protection [14], so it is necessary to consider the ecological externalities generated by agriculture in the measurement of input-output processes [15] and use key indicators of ecological efficiency to measure the production efficiency of agriculture [16]. The concept of eco-efficiency was proposed by German scholars [17], promoted by the World Business Council for Sustainable Development and the OECD, and recognized by scholars and governments [18]. The specific embodiment of the concept of agricultural ecological efficiency is to obtain as much agricultural output as possible with the smallest possible resource consumption and environmental pollution [19] and to ensure the quality of agricultural products [20]. Based on the concept of agricultural production efficiency, it not only needs to pay attention to maximizing agricultural economic benefits (desirable output) [21], but also minimizing resource consumption and environmental damage (undesirable output) to conform to the concept of low-carbon green agricultural development [22].

Scholars have conducted in-depth research on agricultural production efficiency and agricultural ecological efficiency from the micro [23], meso, and macro levels [24], involving index construction [25], influencing factor methods [26], and other aspects [27]. The concept of embedded input-output has been constructed. The input side of the agricultural eco-efficiency measurement index system includes labor, land, capital, fertilizers, pesticides, machinery, etc. [14–17]. The output side is different due to research. Specifically, the agricultural eco-efficiency with undesirable output includes the desirable output of total agricultural output value, total agricultural output, agricultural ecosystem service functions, etc. The undesirable output includes agricultural carbon emissions, agricultural pollution residues, and agricultural greenhouse gas emissions [4,11]. The measurement methods mainly include Stochastic Frontier Analysis (SFA) [28], Data Envelopment Analysis (DEA) [29], Super-efficiency DEA [30], Three-stage DEA [31], Undesirable SBM [32], etc. [33], of which DEA is the basic evaluation method. In addition, the undesirable SBM incorporates negative externalities into the model, which has effectively solved the problem of input-output relaxation and has gradually become the main model for measuring eco-efficiency [34]. Considering that the efficiency measurement method has gradually matured, the undesirable SBM is used in this paper to measure agricultural ecological efficiency. In addition, for the study of the spatial-temporal evolution characteristics of ecological efficiency, most of the literature adopts the method of combining DEA and exploratory spatial data analysis (ESDA) [35], but the ESDA focuses more on the comparative analysis of cross-sectional data in different years.

Therefore, the study of existing literature can be expanded in the following aspects. On the one hand, the SBM model that incorporates undesirable outputs such as agricultural environmental pollution has gradually matured, but the literature using the undesirable SBM model has less consideration for the further comparison of the decision-making unit with an efficiency of 1 [12–14], and the super-efficient SBM model can effectively solve this

problem. On the other hand, considering the limitations of the existing ESDA method for spatial panel data research, there is little literature focusing on predicting the evolution of agricultural ecological efficiency, which is the problem that the spatial Markov Chain (MC) can solve. Based on this, the provincial panel data of China from 2000 to 2019 is chosen as the research unit, the agriculture in the narrow sense (planting industry) is regarded as the research object, the agricultural carbon emissions are considered the undesirable output, the super-efficient SBM is used to analyze the agricultural ecological efficiency, and the traditional and spatial Markov probability transfer matrix based on time series analysis and spatial correlation analysis are constructed. Through the comparison of the transfer matrix, the spatial-temporal dynamic evolution characteristics of agricultural ecological efficiency can be analyzed, and the long-term evolution trend can be predicted, which aims to provide support for narrowing the agricultural ecological efficiency between various provinces.

#### **2. Methods and Data**

#### *2.1. Super-Efficient SBM Model Based on the Undesirable Output*

In the process of agricultural production, it is generally expected that the environmental pollution caused by the excessive use of chemical products such as fertilizers, pesticides, and agricultural films should be as small as possible [36]. The SBM model based on undesirable output is the model firstly proposed by Tone in 2001 to measure eco-efficiency. The SBM model can effectively solve the "crowding" or "relaxation" phenomenon of input elements caused by the radial and angled traditional data envelopment model (DEA) model [37], but the SBM has the same problem as the traditional DEA, that is, for DMU with an efficiency of 1, it is difficult to further distinguish the difference between these efficient DMUs. Based on the SBM model, Tone further defined the super-efficient SBM. It is a model combining the super-efficiency of the DEA and SBM models, which integrates the advantages of the two models. The super-efficiency SBM model can further compare and distinguish the efficient DMU in the forefront compared to the general SBM model, and the model is constructed as follows.

$$Min\rho = \frac{\frac{1}{m}\sum\_{i=1}^{m} \left(\frac{\overline{x}}{x\_{ik}}\right)}{\frac{1}{r\_1 + r\_2} \left(\sum\_{s=1}^{r\_1} \frac{\overline{y}^d}{y\_{sk}^d} + \sum\_{q=1}^{r\_2} \frac{\overline{y}^u}{y\_{qk}^u}\right)}\tag{1}$$

$$\left\{ \begin{array}{c} \overline{\boldsymbol{\pi}} \ge \sum\_{j=1, \neq k}^{n} \boldsymbol{\lambda}\_{j} \boldsymbol{\lambda}\_{j}; \overline{\boldsymbol{\mathcal{Y}}}^{d} \le \sum\_{j=1, \neq k}^{n} \boldsymbol{y}\_{sj}^{d} \boldsymbol{\lambda}\_{j}; \overline{\boldsymbol{\mathcal{Y}}}^{d} \ge \sum\_{j=1, \neq k}^{n} \boldsymbol{y}\_{qj}^{d} \boldsymbol{\lambda}\_{j}; \overline{\boldsymbol{\mathcal{X}}} \ge \boldsymbol{x}\_{k}; \overline{\boldsymbol{\mathcal{Y}}}^{d} \le \boldsymbol{y}\_{k}^{d}; \overline{\boldsymbol{\mathcal{Y}}}^{\boldsymbol{\mu}} \ge \boldsymbol{y}\_{k}^{\boldsymbol{\mu}} \right. \tag{2}$$
  $\boldsymbol{\lambda} \ge \boldsymbol{0}, \boldsymbol{i} = 1, 2, \dots, m; \boldsymbol{j} = 1, 2, \dots, n$ ;  $\boldsymbol{\eta} \ne 0$ ;  $\boldsymbol{s} = 1, 2, \dots, r; \boldsymbol{q} = 1, 2, \dots, r\_{2}$ 

where it is assumed that there are *n* decision units, each decision unit has input *m*, desirable output *r*1, and undesirable output *r*2. *x*, *yd*, *yμ*, respectively represent the corresponding elements of the input matrix, desirable output matrix, and undesirable output matrix and *ρ* is the ecological efficiency values.

Agricultural ecological efficiency refers to obtaining the greatest possible agricultural economic output and ecological protection with as little agricultural resource input and low environmental cost as possible, which comprehensively reflects the coordinated and winwin relationship between the agricultural economy, resource utilization, and environmental protection. An indicator system of agricultural ecological efficiency in China is constructed in this paper; land, labor, mechanical power, irrigation, fertilizer, pesticide, etc. are regarded as the input index of regional agricultural resources, while the total agricultural output value is taken as the desirable output index, and the agricultural carbon emission is chosen as the undesirable output (Table 1).


**Table 1.** Agricultural ecological efficiency index system.

Undesirable output: carbon emissions. The following carbon emission sources and their emission coefficients were selected: fertilizer 0.8956 (kg/kg), pesticide 4.9341 (kg/kg), agricultural film 5.18 (kg/kg), diesel 0.5927 (kg/kg), agricultural sowing 312.6 (kg/km2), agricultural irrigation 20.476 (kg/km2).

#### *2.2. Nonparametric Kernel Density Estimation*

Kernel density estimation is a kind of density mapping. In essence, it is a process of surface interpolation through discrete sampling points; that is, through the smoothing method, the histogram is replaced by continuous density curves to better describe the distribution pattern of the variables, which has excellent statistical characteristics and is more accurate and smoother than the histogram estimation. As a nonparametric estimation method, Kernel density estimation can use a continuous density curve to describe the distribution pattern of random variables. The density function of random variables is set to *f*(*x*), for the random variable *y*, there are *n* independent observations of the same distribution, with *y*1, *y*2,..., *yn*, respectively. The Kernel density function is as follows.

$$f(x) = \frac{1}{nh} \sum\_{i=1}^{n} K(\frac{y\_i - y}{h}) \tag{3}$$

where *n* is the number of study areas, *h* is the window width, and *K*(•) is a random Kernel function, which is a weighted function or smoothing function, including Gaussian (Normal) Kernel, Epanechnikov Kernel, Triangular Kernel, Quartic Kernel, etc. The window width determines the degree of smoothness of the estimated density function. To be specific, the larger the window width, the smaller the variance of the kernel estimation, the smoother the density function curve, but the greater the estimated deviation. Therefore, the choice of the best window width must be weighed between the variance and deviation of the kernel estimation so that the mean squared error can be minimized. The kernel density function of the Gaussian Kernel distribution is used in this paper, with the window width set to *h* = 0.9*SeN*<sup>−</sup>0.2.

#### *2.3. Spatial Correlation Analysis*

According to the first law of geography, each thing or phenomenon in the spatial unit is not isolated, but related. The degree of connection between neighboring things or phenomena is closer, and the agricultural production activities in adjacent areas may affect each other more. Spatial autocorrelation can indicate the impact of neighboring areas, while differences in regional spatial distribution may have spatial autocorrelations; that is, the geographical location of a region affects not only its own agricultural ecological efficiency, but also affects the efficiency of its neighbors. In this case, it is necessary to measure the spatial autocorrelation of regional agricultural ecological efficiency. Spatial autocorrelation analysis includes global spatial autocorrelation and local spatial autocorrelation. Global spatial autocorrelation is used in this study to clarify the spatial correlation and spatial difference of regional agricultural ecological efficiency. In spatial statistics, the commonly used statistical indicator to measure the degree of spatial autocorrelation is Moran's I index, and its calculation formula is as follows.

$$I = \frac{n\sum\_{i=1}^{n}\sum\_{j=1}^{n}w\_{ij}(\mathbf{x}\_{i}-\overline{\mathbf{x}})\left(\mathbf{x}\_{j}-\overline{\mathbf{x}}\right)}{\sum\_{i=1}^{n}\sum\_{j=1}^{n}w\_{ij}(\mathbf{x}\_{i}-\overline{\mathbf{x}})^{2}}\tag{4}$$

where *n* is the sample size, *xi* and *xj* are the observations of the spatial positions *i* and *j*. *wij* represents the proximity of spatial positions *i* and *j* when *i* and *j* are adjacent, *wij* = 1; otherwise, it is 0. The global Moran's I index has a range of values with [−1,1], greater than 0 means spatial positive correlation, less than 0 means negative correlation, and equal to 0 means uncorrelation.

#### *2.4. Space Markov Chains*

The traditional Markov Chain is derived from the Russian mathematician Markov's theory of stochastic processes, which measures the state of occurrence of events and their development trend by constructing a state transition probability matrix. In which, given current knowledge or information, the past (i.e., the historical state before the current period) is irrelevant to predicting the future (i.e., the future state after the current period); that is, "no after-effect", also known as "Markovity". The evolution of many economic phenomena has the characteristic of "no after-effect", and the evolution process of agricultural ecological efficiency also has no exception. Assuming that *Pij* is the transfer probability of agricultural ecological efficiency in a certain area from state *i* in year *t* to state *j* in year *t* + 1, the transfer probability can be estimated by the frequency of transfer *Pij* = *nij*/*ni*. Where *nij* represents the number of provinces that have shifted from state *i* in year *t* to state *j* in year *t* + 1 during the sample investigation period and satisfies the formula (∑*<sup>j</sup> Pij* = ∑*<sup>j</sup> P*{ *Xn*+<sup>1</sup> = *j*|*Xn* = *i*} = 1). If the agricultural ecological efficiency is divided into *N* types, that is, *N* states, the state transition probability matrix *N* × *N* can be constructed. In addition, the transfer direction can be defined according to the upward (increasing), constant, and downward (decreasing) changes of agricultural ecological efficiency types (Table 2).


**Table 2.** Markov chain transfer probability matrix (*N* = 4).

Spatial Markov Chain analysis introduces the concept of spatial lag into the transfer probability matrix because the evolution of regional economic growth and other economic phenomena is not geographically isolated and random, but closely related to and interacts with neighboring regions. Spatial MC compensates for the neglect of the spatial correlation effects of traditional MC analysis on the study area and is used to reveal the intrinsic relationship between the spatial-temporal evolution of an economic phenomenon and the spatial background of the region. Taking the spatial lag type of a region in the initial year as the condition, the traditional *N* × *N* state of the transition probability matrix is decomposed into an *N* transfer condition probability matrix (*N* × *N*). It enables the analysis of the possibility of improving or decreasing the agricultural ecological efficiency of a certain region under different geographical background conditions. Taking the *N* condition matrix as an example, *Pij*(*N*) indicates the spatial lag type *N* of a certain region in year *t*, and the one-step spatial transfer probability of the year is transferred from state *i* in year *t* to state *j* in year *t* + 1. The spatial lag type of a region is classified by the spatial lag value of its ecological efficiency value in the initial year, and the spatial lag value is the spatially weighted average of the agricultural ecological efficiency value in the adjacent area of the region, which is calculated by the product of the regional agricultural ecological efficiency value and the spatial weight matrix, that is ∑*<sup>j</sup> WijYj*. *Yj* represents the ecological efficiency value of a certain region and *Wij* represents the element of the spatial weight matrix *W*. The principle of common boundaries is used in this paper to determine the spatial weights matrix. If the regions are adjacent, *Wij* = 1; otherwise, *Wij* = 0. Due to the special geographical location of Hainan Province, the calculation of the weight matrix assumes that Hainan is adjacent to Guangdong.

By comparing the corresponding elements of the traditional Markov transfer matrix and the spatial Markov transfer matrix (Table 3), the relationship between the magnitude of the probability of upward or downward transfer of agricultural ecological efficiency in a certain area and the surrounding neighborhood can be understood, so as to explore the differences and the spatial spillover effects between spatial contexts and the transfer of agricultural ecological efficiency. If *P*<sup>12</sup> > *P*12/1, it means that the probability of a province's ecological efficiency shifting from state 1 to state 2 without considering that the neighborhood is greater than the probability of transferring from state 1 to state 2 when considering the neighborhood and the province is adjacent to the province in state 1. If the spatial background has no effect on the state transfer, there is *P*<sup>12</sup> = *P*12/1.


**Table 3.** Spatial Markov transfer condition probability matrix (*N* = 4).

The Markov process occurs after a long transition. If the system has an equilibrium state, i.e., the probability that the system is in the same state when it is balanced, it does not depend on the initial state and no longer changes over time. The probability distribution at this time is a stationary distribution. Based on the Markov probability transfer matrix, a smooth distribution of the stochastic process can be obtained, which can predict the dynamic evolution trend of an economic phenomenon (agricultural ecological efficiency). It is assumed that the traditional Markov chain is {*Xn*, *n* = 0, 1, 2, . . .} , *pij* is a one-step transfer probability, {*π*, *i* ∈ *S*} is a smooth distribution of the traditional Markov chain. Extended to spatial Markov chains, according to the similarity principle, the spatial stationary distribution under each spatial lag type is obtained, and the maximum transition probability is taken as the possible evolution trend of the corresponding state.

#### *2.5. Research Data*

Agriculture in the broad sense includes agriculture, animal husbandry, and fishery, and in the narrow sense refers to the planting industry. This paper measured the agricultural ecological efficiency in the narrow sense, and the research objects were 31 provinces, autonomous regions, and municipalities directly under the central government (excluding Hong Kong, Macao, and Taiwan). Starting from 2000, a total of 20 years was used as the study interval to 2019. The agriculture-related data used in this paper were derived from the "China Rural Statistical Yearbook", the statistical yearbook of various provinces, and the carbon emission data was used in the calculation.

Through the spatialization analysis of the total power of agricultural machinery, fertilizer input, land input, and labor input, there was spatial heterogeneity between the inputs in different provinces. Specifically, in 2020 (Figure 1), the provinces with a total mechanical power greater than 6000 KWH were mainly clustered in the main grainproducing areas, of which Heilongjiang, Henan, Shandong, and other provinces were representatives. The provinces with fertilizer inputs greater than 300 tons were Henan and Shandong; the provinces with larger agricultural land inputs were Henan, Shandong, Sichuan, and Heilongjiang; the provinces with large population inputs were Shandong, Henan, Guangxi, and Sichuan.

**Figure 1.** Input factors for provincial agricultural production in China in 2020. (**a**) Machine input; (**b**) fertilizer input; (**c**) land input; (**d**) population input.

#### **3. Results**

#### *3.1. Measurement and Temporal and Spatial Pattern of Agricultural Ecological Efficiency in China*

DEA SOLVER PRO 5.0 (Beijing, China) software, non-Oriented, and super-efficient SBM model with variable returns to scale were used to measure the agricultural ecological efficiency of 31 provinces, autonomous regions, and municipalities directly under the central government in China from 2000 to 2019 (Figure 2). The agricultural ecological efficiency in each year was basically below 0.7, and the overall agricultural ecological efficiency in China was at a low level, which means that there is a large space for resource conservation and environmental protection in the two-type development of China's agriculture. Between 2001 and 2019, the agricultural ecological efficiency of China showed a steady upward trend in fluctuations, especially from 2006 to 2010. Since the 21st century, the No. 1 central document has focused on agriculture for many years, paying attention to the issues of "Agriculture, rural areas and farmers", and clearly proposing to "encourage the development of circular agriculture and ecological agriculture", showing the government's emphasis on the sustainable development of agriculture and improved agricultural ecological efficiency. By comparing the agricultural ecological efficiency of the three regions in the east, central, and west, the agricultural ecological efficiency ranked in order of the eastern, the western, and the central, but the gap between the central and western regions was small, showing that with the continuous development of the agricultural economy, the agricultural technology level in the eastern region made significant progress. More attention has been paid to agricultural modernization, the coordination between agricultural production, resource conservation, and environmental protection. However, the level of agricultural technology in the central and western regions has developed slowly, the degree of agricultural mechanization is low, and the mode of agricultural economic development is relatively extensive.

In order to explore the agglomeration differences in the evolution of agricultural ecological efficiency over time in various provinces, the nonparametric Kernel density function of the Gaussian normal distribution was used. The estimated Kernel density at five time points of observation in 2001, 2005, 2010, 2015, and 2019 was selected to obtain the distribution of different time points. The height of the crest reflects the degree of agglomeration of agricultural ecological efficiency in various provinces, and it can be seen in (Figure 3) that the agricultural ecological efficiency of China as a whole presents the evolution and distribution characteristics of a "double peak" from left to right and crests from high to low, showing the increasing trend of China's agricultural ecological efficiency with time. Most provinces have gradually changed from low-level agglomeration to the trend of narrowing the difference between "high-low" quantities. In 2001, the agricultural ecological efficiency of most provinces was at low levels. After 2010, with the enhancement of agricultural environmental awareness, the transfer of rural labor, and the acceleration of the progress of agricultural mechanization, the agricultural ecological efficiency of various provinces has improved to varying degrees, but there are still differences in resource endowments and economic strength among provinces. The gap in agricultural ecological efficiency in various provinces and cities has begun to increase, forming a number of peaks of different amplitudes. However, the peaks of low-level agglomeration have gradually declined until 2019. The peak height gap of the bimodal distribution has narrowed, indicating that the gap between low-level agricultural ecological efficiency and high-level agricultural ecological efficiency has been further narrowed, and a "bimodal" evolution pattern similar to "club convergence" of "low-level agglomeration and high-level agglomeration" has been gradually formed.

**Figure 2.** Provincial agricultural ecological efficiency from 2001 to 2019 (x-axis: agricultural gross domestic production).

**Figure 3.** Kernel density function plot from 2001 to 2019.

#### *3.2. Temporal Evolution Characteristics of Agricultural Ecological Efficiency in China*

The temporal analysis of agricultural ecological efficiency and the density Kernel analysis can only characterize the temporal change trend and evolution difference of ecological efficiency, but it is hard to deeply reflect its inherent spatial-temporal evolution law. To solve the probability transfer matrix of the traditional Markov chain, considering the provincial view measurements of each type are roughly the same, according to the quantile division method, the agroecological efficiency values are divided into four complete intervals that are adjacent but not intersecting according to the quantile division method, with the quantiles of 1/4, 1/2, and 3/4 as the boundaries. The complete intervals of these four states are defined as *k* = 1, 2, 3, 4, respectively; the larger the *k,* the higher the agricultural ecological efficiency of the region. In addition, according to the trend evolution diagram, the whole research period is roughly divided into two stages: 2001–2010 and 2011–2019. According to the division of state types, the traditional Markov probability transition matrix is obtained (Table 4).

**Table 4.** Markov Chain state transition probability matrix of agricultural ecological efficiency in China from 2001 to 2019.


The elements on the diagonal line indicate the probability of the non-transfer of agricultural ecological efficiency between different state types, reflecting the stability characteristics of agricultural ecological efficiency in the region. The elements on the non-diagonal line indicate the probability of the transfer of agricultural ecological efficiency between different state types, showing the evolution characteristics of agricultural ecological efficiency

without considering the geospatial pattern. The agricultural ecological efficiency of each province has the stability of maintaining its original state. All elements on the diagonal are significantly larger than those on the non-diagonal. The minimum value of the elements on the diagonal is 0.7643, and the maximum value is 0.9872. It means, regardless of the period, the agricultural ecological efficiency of a certain province belongs to a certain type in the current year and the probability of still belonging to the type in subsequent years. In addition, the types at both ends of the diagonal have the greatest probability of maintaining stability (type 1 and 4), and there is a trend of agricultural ecological efficiency converging to low and high levels, that is, "club convergence". On the whole, agricultural ecological efficiency has a significant trend of transferring to a high level. In view of the transfer frequency of the two stages, the number of type 1 and type 2 in 2011–2019 was significantly less than that in 2001–2010, while the number of type 3 and type 4 was more than that in 2001–2010. However, there were also differences in transfers between different phases. From 2001 to 2010, the probability gap of the upward or downward transfer between type 2 and type 3 was small, and the gap between *P*<sup>22</sup> and *P*<sup>33</sup> was also small, indicating that the change of agricultural ecological efficiency of provinces in this stage was relatively stable. While during 2011–2019, *P*<sup>21</sup> =0< *P*<sup>23</sup> = 0.1638, *P*<sup>32</sup> = 0.027 < *P*<sup>34</sup> = 0.0495, the provincial agriculture ecological efficiency was more likely to increase during this period. Besides, it is difficult to achieve a leapfrog transfer of agricultural ecological efficiency. The transfer of agricultural ecological efficiency occurs almost on both sides of the diagonal line, and from the perspective of non-diagonal elements, the probability on the non-diagonal line is significantly smaller than the probability on the diagonal line, and its maximum value is 0.1638, which means that between two consecutive years, the probability of achieving a leapfrog transfer of agricultural ecological efficiency is extremely small. It proves that the improvement of agricultural ecology in various provinces is a relatively stable and sustained process, and achieving leapfrog development and evolution in the short term is relatively hard.

#### *3.3. Spatial Evolution Characteristics of Agricultural Ecological Efficiency in China*

With the improvement of the agricultural market economy and the expansion of regional openness in China, the spatial mobility of agricultural production factors has become more frequent, the spatial connection between agricultural production has become closer, and the location effect of agricultural ecological efficiency between adjacent provinces will become more significant. In order to deeply research the temporal and spatial evolution differences of agricultural ecological efficiency in different provinces, the spatial Markov chain probability transfer matrix is constructed. On the premise of whether there is a spatial background effect, by comparing the corresponding elements in the two matrices, a relationship was found between the transfer probability of agricultural ecological efficiency types of a province and its neighboring provinces. The spatial correlation test results of the global Moran's I index of agricultural ecological efficiency show that the Moran's I index of agricultural ecological efficiency in different years was positive (0.1245~0.2780). Except for a few years, all of them passed the test at the significance level of 5% (or 10%), which showed that there was a significant positive correlation in the spatial distribution of agricultural ecological efficiency in China. The geospatial spatial pattern is an important factor affecting the agricultural ecological efficiency of China, and the impact of agricultural ecological efficiency between neighboring provinces was spatially dependent. From 2001 to 2010, 28 provinces witnessed an increase in agricultural ecological efficiency. From 2011 to 2019, 7 provinces saw an increase in it (Figure 4). Although the number of growing provinces decreased, the gap between provinces became smaller.

**Figure 4.** Changes in agricultural ecological efficiency in China from 2001 to 2019: (**a**) 2001–2010; (**b**) 2011–2019.

The traditional Markov probabilistic transfer matrix does not take into account the effects of the transfer of the surrounding neighborhood type, but the upward or downward transfer of agricultural ecological efficiency is not spatially isolated; however, it is paved with the surrounding neighborhood and practically related. On the basis of the traditional Markov probability transfer matrix, the geospatial background factor is introduced. The spatial Markov probability transfer matrix is constructed based on the spatial lag type of each province in the initial year. The spatial Markov probability transfer matrix is constructed on the condition that the spatial lag type in different regions in the initial year is obtained. The division of the two periods in the above study shows that the type transfer of agricultural ecological efficiency is unstable in time, and the 2001–2010 and 2011 are established, respectively (Table 5).

**Table 5.** Spatial Markov probability transfer matrix of agroecological efficiency in China from 2001 to 2019.


#### **4. Conclusions**

This paper takes the provincial panel data of China from 2001 to 2019 as the research unit, regards the narrow agriculture (planting industry) as the research object, uses the super efficiency SBM model to measure the inter-provincial agricultural ecological efficiency, and constructs the traditional and spatial Markov probability transfer matrix based on the time analysis of Kernel density estimation and the spatial correlation analysis of global Moran's I index. Through the comparative analysis of the matrix, the temporal and spatial dynamic evolution characteristics of agricultural ecological efficiency are analyzed, and the influence of geospatial patterns on the temporal and spatial evolution of agricultural ecological efficiency is concluded. The main conclusions are as follows.

From the perspective of the difference in temporal evolution, the agricultural ecological efficiency of China based on line charts and Kernel density estimation shows a steady upward trend in fluctuations, and the volatility was mainly clustered around 2001–2019. From 2011 to 2019, there was still much space to improve agricultural ecological efficiency. The Kernel density estimation chart shows that the crest of the evolution of agricultural ecological efficiency in China has the distribution characteristics of a "double peak" from high to low, and the gap in wave crest height has narrowed. A "club convergence" pattern of "low agglomeration and high agglomeration" is gradually formed. This result is the same as Wang (2022) [4].

From the perspective of spatial evolution, the results of the global Moran's I index show that there was a significant positive correlation between the agricultural ecological efficiency of China and spatial distribution, and there was spatial dependence between interprovincial agricultural ecological efficiency. The traditional Markov probability transfer matrix presents that the overall trend of agricultural ecological efficiency shifting to a high level is significant, the evolution of agricultural ecological efficiency has the stability of maintaining the original state, and it is difficult to achieve leapfrog transfer. The probability of the type at both ends of the matrix remaining unchanged is the largest, and there is a possibility of "low agglomeration and high agglomeration". The comparison between the spatial Markov probability transfer matrix and the traditional transfer matrix shows that, in addition to common characteristics, the geospatial pattern plays an important role in the spatial-temporal evolution of agricultural ecological efficiency and the spatial spillover effect is obvious. In different spatial contexts, the probability of the transfer of agricultural ecological efficiency types in different provinces has differences. If a province is adjacent to a province with a high agricultural ecological efficiency, the probability of upward transfer increases, while if it is adjacent to a province with low agricultural ecological efficiency, the probability of downward transfer increases. It means that provinces with high agricultural ecological efficiency enjoy positive spillover effects, while provinces with low agricultural ecological efficiency have negative spillover effects, so as to gradually form a "club convergence" of "high agglomeration, low agglomeration, high radiation and low inhibition" in the spatial pattern, which can echo the characteristics of time evolution. This result is the same as Guo (2020) [14] and Liu (2020) [15].

#### **5. Discussion**

Based on the super-efficiency SBM model and the spatial Markov probability transfer matrix, the spatial and temporal evolution characteristics of agricultural ecological efficiency in China are summarized as follows.

The super-efficient SBM model based on undesirable output can fully consider the impact of the resource environment on agricultural ecological efficiency and can further compare DMU with an efficiency of 1, but this paper does not discuss the various socioeconomic factors that lead to the loss of agricultural ecological efficiency, because the loss of agricultural ecological efficiency is the result of the joint influence of various factors [38]. Future research can consider structural changes, policy changes, technological changes, urbanization [39], foreign direct investment (FDI), and other factors, as well as profoundly discussing the influencing factors of agricultural ecological efficiency loss.

The long-term evolution trend of China's agricultural ecological efficiency shows that there is a huge potential for improvement in provinces with low agricultural ecological efficiency. The agricultural modernization of China is still facing the arduous task of resource conservation and environmental protection [40]. In accordance with the requirements of green and low-carbon development, it is necessary to continue to deeply change the mode of agricultural economic development, pay attention to resource conservation and protection and the control of agricultural pollution emissions, learn from the experience of ecological agriculture management in neighboring provinces with high agricultural ecological efficiency, and reduce the provincial gap of agricultural ecological efficiency.

**Author Contributions:** Conceptualization, M.C. and W.S.; methodology, M.C.; software, M.C. and L.M.; validation, W.S; formal analysis, M.C.; investigation, M.C.; resources, W.S.; data curation, M.C.; writing—original draft preparation, M.C.; writing—review and editing, W.S. and Q.G.; visualization, M.C.; supervision, G.W. and W.S.; project administration, G.W. and W.S.; funding acquisition, W.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 72003111.

**Institutional Review Board Statement:** Not applicable.

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

#### **References**

