Grassland greening on the Mongolian Plateau despite higher grazing intensity

Changes in land management and climate alter vegetation dynamics, but the determinants of vegetation changes often remain elusive, especially in global drylands. Here we assess the determinants of grassland greenness on the Mongolian Plateau, one of the world's largest grassland biomes, which covers Mongolia and the province of Inner Mongolia in China. We use spatial panel regressions to quantify the impact of precipitation, temperature, radiation, and the intensity of livestock grazing on the normalized difference vegetation indices (NDVI) during the growing seasons from 1982 to 2015 at the county level. The results suggest that the Mongolian Plateau experienced vegetation greening from 1982 to 2015. Precipitation and animal density were the most influential factors contributing to higher NDVI on the grasslands of Inner Mongolia and Mongolia. Our results highlight the dominant effect of climate variability, and especially of the precipitation variability, on the grassland greenness in Mongolian drylands. The findings challenge the common belief that higher grazing pressure is the key driver for land degradation. The analysis exemplifies how representative wall‐to‐wall results for large areas can be attained from exploring space–time data and adds empirical insights to the puzzling relationship between grazing intensity and vegetation growth in dryland areas.

Climate change, with its effects on changes in temperature, rainfall volatility, and incoming shortwave radiation, has been deemed responsible for higher variations in vegetation growth in many global dryland areas (Gherardi & Sala, 2019;Huang et al., 2018). It has also been suggested that vegetation growth increases due to relatively higher nitrogen deposition and increasing amounts of carbon dioxide (CO 2 ; Huang et al., 2018).
However, evidence regarding the drivers for changes in grassland resources in global drylands is scarce, often contradictory, and mainly rests on results obtained from experimental sites (Dangal et al., 2016;Shinoda, Nachinshonhor, & Nemoto, 2010). The quantification of anthropogenic grassland management, such as livestock grazing, and the effects of diverging land-use policies on grassland management have rarely been addressed (Craine et al., 2013;Zhu, Chiariello, Tobeck, Fukami, & Field, 2016). This is unfortunate because an improved understanding of the determinants of changes in grassland greenness is important for informing policies that aim to balance livelihood outcomes and environmental effects. The grassland greenness is often approximated by the normalized difference vegetation index (NDVI) values obtained from remote sensing image time-series, as NDVI is a good proxy for vegetation growth, vegetation productivity, and vegetation greenness Pinzon & Tucker, 2014). The greening of grassland (i.e., an increase in NDVI) typically suggests an improvement in vegetation productivity over time, whereas grassland browning (i.e., a decrease in NDVI) implies a degradation of vegetation productivity .
The Mongolian Plateau is one of the largest contiguous global drylands and stretches approximately 2.75 million km 2 from western Mongolia, bordering Kazakhstan, to the eastern part of the Inner Mongolia Autonomous Region in China (Hilker, Natsagdorj, Waring, Lyapustin, & Wang, 2014). The region is famous for its extensive grasslands that support unique nomadic cultures (Humphrey, Sneath, & Sneath, 1999). The ecological conditions in the region are fragile due to the prevalent arid and semiarid climate, which results in a low potential for crop production and weak grassland growth (Wang, Brown, & Agrawal, 2013). In recent decades, multifaceted pressures, such as climate change and overgrazing, have arguably resulted in widespread land degradation, soil erosion, desertification, and the disappearance of rivers and lakes (Batunacun, Nendel, Hu, & Lakes, 2018;Tao et al., 2015). However, the effects of changes in grazing pressure and climate on the grassland growth in the Mongolian Plateau remain uncertain, and comprehensive analysis has been lacking to date.
Climate change is of particular concern for the sustainable development of grasslands on the Mongolian Plateau because the region has experienced more rapid climate warming than the global average, with an annual average temperature increase of 0.4 C every 10 years during the last 40 years (Pederson, Hessl, Baatarbileg, Anchukaitis, & Di Cosmo, 2014;Tao et al.,2015). Limited and volatile rainfall and, more importantly, high evapotranspiration have contributed to increasing aridity, widespread land degradation, soil erosion, and severe water scarcity (Nandintsetseg & Shinoda, 2014). These changes threaten the livelihoods of millions of herders who rely on these grasslands to feed their animals (Wang et al., 2013). Future challenges are daunting, as the growing season temperature is expected to increase by approximately 3 C between 2010 and 2100 according to the midrange scenario of the Intergovernmental Panel on Climate Change (Miao, Ye, He, Chen, & Cui, 2015). The rising evapotranspiration and plant transpiration from the strong temperature increase may negatively affect the growth of the grassland. Conversely, the projected gradual increase in precipitation may contribute to the greening of the grassland. However, the extent to which climate change has affected and will affect the grassland growth remains uncertain.
Animal husbandry remains the foundation of many rural livelihoods on the Mongolian Plateau and plays a dominant role in rural economies (Cease et al., 2015). Despite the similar biophysical conditions in Mongolia and Inner Mongolia, the density of grazing livestock in Inner Mongolia is traditionally much higher than that in Mongolia because of the higher population density, along with the higher demand for livestock products (Angerer, Han, Fujisaki, & Havstad, 2008;Sneath, 1998). In both China and Mongolia, the livestock density has increased substantially since 2000 in response to the rising demand for livestock products brought about by population growth and income increases (Hilker et al., 2014). The increasing density of grazing livestock has resulted in higher biomass extractions from pastures, which have arguably led to widespread land degradation in many areas (Middleton, 2016). However, the effect of increasing grazing pressure on grassland resources on the Mongolian Plateau remains a subject of debate (Hilker et al., 2014;Tian, Herzschuh, Mischke, & Schlütz, 2014).
The grazing patterns across the region are shaped by institutional regulations. Mongolia continues to allow the free movement of grazing livestock in line with traditional nomadic grazing practices. As a result, the majority of herders in Mongolia continue to adhere to nomadic and seminomadic lifestyles (Sneath, 1998;Wang et al., 2013). In contrast, livestock producers in China conduct sedentary grazing under the policy of the household responsibility system, under which the land-use rights of croplands and grasslands are allocated to individual households, and producers can graze animals only on their own lands. Moreover, the Chinese government has implemented various ecological protection measures, such as the fencing policy since 2000, which further regulates the grazing patterns of livestock (Wang et al., 2013). The rationale behind these policies has been to avoid common-pool resource problems (Gardner, Ostrom, & Walker, 1990), that is, to prevent land degradation and desertification caused by the excessive grazing of unmanaged grassland resources. The Chinese government considered the vicious circle of resource degradation in a situation that can best be characterized as 'open-access' as a severe threat to the sustainable development of the livestock sector. Ironically, the reduction of livestock mobility through the individualization of land-use rights is suspected to have contributed to a higher degree of grassland degradation in Inner Mongolia and thus may have jeopardized the long-term welfare of local herders (Sneath, 1998) Figure 1). Grasslands cover more than 60% of the area. The region has an arid and semiarid continental climate with extremely cold winters, dry and hot summers, and recurring droughts and dzuds (severe winter weather; . The average temperature ranges from −45 to 35 C, and the mean annual precipitation is approximately 200 mm. In both Mongolia and Inner Mongolia, grazing is the dominant land use, and products generated from livestock are the major income source for the rural population (Zhang et al., 2018). The sustainable use of grassland resources is therefore a priority for land-use policy in the region (Liu et al., 2018;Wang et al., 2013). See Figures S1 and S2 for the field conditions.

| Land cover
To focus our analysis on all areas that can potentially be used for grazing, we used the land cover classifications derived from the moderate resolution imaging spectroradiometer (MODIS; MCD12C1.V006) from 2001 to 2015 at a spatial resolution of 0.05 (ca. 5 km, downloaded from https://lpdaac.usgs.gov/products/mcd12c1v006/). We merged the categories of grassland and savanna from this MODIS land cover product, along with closed and open shrublands, into a single category that captures potential grassland. To arrive at a conservative estimate of the grassland extent and to reduce the confounding effects of changes in the extent of the grassland area, we selected all pixels from the MODIS land cover map that were grasslands throughout the whole study period, namely, between 2001 and 2015. We considered the variation in greenness only for this grassland mask in the subsequent analysis.

| Normalized difference vegetation index
We used the third-generation NDVI data (NDVI3g.v1) from NASA's GIMMS group (https://nex.nasa.gov/nex/projects/1349/). This dataset was derived from the advanced very-high-resolution radiometer (AVHRR) onboard several satellites from the National Oceanic and Atmospheric Administration (NOAA). The GIMMS NDVI 3g.v1 dataset has a spatial resolution of approximately 8 km and provides the F I G U R E 1 Location of the study area and land cover. The administrtive boundaries are from https://gadm.org/data.html. Land cover information was extracted from the IGBP classification of MODIS land cover product, available from https://modis.gsfc.nasa.gov/data/dataprod/ mod12.php [Colour figure can be viewed at wileyonlinelibrary.com] longest time-series of data monitoring vegetation growth and vegetation phenology (Fensholt & Proud, 2012). We analyzed the complete temporal coverage from 1982 and until 2015 with 15-day composites.
We excluded pixels with mean annual NDVI values of less than 0.1, because these regions are likely covered with sparse vegetation or bare soil and are thus not able to support substantive numbers of livestock. To determine the growing season of vegetation for each pixel during each year, we used polynomial curve fitting, which allowed us to calculate the start and end of the growing season and thus its length (Miao, Müller, Cui, & Ma, 2017;Piao, Mohammat, Fang, Cai, & Feng, 2006). To ensure comparability among the years, we used the average starting date and end date of the growing season over the study period and extracted NDVI data during this time-averaged growing season.
To approximate grassland greenness, we calculated the mean NDVI during the growing season for each year beginning in 1982, which was the first year for which livestock data were available, through 2015, when the time-series of livestock data that we used ends. We thus aggregated the NDVI to the county level (the unit of our analysis, see Figure S3) by averaging the mean NDVI during the growing season for all pixels within the grassland mask. We used the county level as the unit of analysis in our regressions because this is the lowest administrative level in both China and the Republic of Mongolia for which there are available statistical data. Moreover, we performed a linear regression using ordinary least squares between the NDVI and time (i.e., year) and regarded the slope as the temporal trend of the NDVI. We used least-squares method to estimate the slope and t test to determine whether it is statistically significant at p < .05.

| Grazing density
We obtained livestock statistics for Mongolia from the National Statistical Office of Mongolia (http://www.en.nso.mn/index.php) and for Inner Mongolia from statistical yearbooks. All livestock statistics were available at the county level ( Figure S3). The data include 74 counties in Inner Mongolia (ranging in size from 103 to 85,089 km 2 with an average size of 12,518 km 2 ) and 263 counties in Mongolia (ranging from 101 to 28,211 km 2 with an average size of 4,673 km 2 ). Among them, 14 counties in Inner Mongolia and 64 counties in Mongolia were excluded due to the absence of grassland coverage or livestock data ( Figure S3). We extracted the livestock species that are predominantly nourished by grazing on grassland, namely, cattle, horses, camels, donkeys, mules, sheep, and goats. To make the grazing pressures from these grazing livestock species comparable, we converted all livestock numbers into animal unit (au) equivalents, with 1 au was defined relative to one mature cow (i.e., one cow = 1.00 au). Other conversion factors were one horse = 1.80 au, one sheep = 0.15 au, one goat = 0.10 au, one camel = 1.25 au, one donkey = 1.05 au, and one mule = 0.15 au (Holechek, 1988). We multiplied all animal numbers by their respective conversion factors to obtain the animal units for each type of animal. We then summed the animal units of all types and those for each year and divided them by the area of grassland in every county to approximate the annual grazing density across the study area.
We resampled the climate data and MODIS land cover product (i.e., grassland layer) into the same spatial resolution (8 km) as GIMMS NDVI 3g.v1 . To maintain consistency among different datasets, both the climate data and NDVI data were clipped with the same grassland mask. We then aggregated all variables to the county level for the spatial panel analysis.

| Quantifying determinants for changes in grassland greenness
We used spatially explicit panel regressions to assess the drivers of changes in the annual NDVI at the county level due to the presence of spatial autocorrelation in the variables used in this analysis (the Moran's I test results can be found in Table S2). The NDVI, which was the dependent variable in our regressions, was strongly spatially clustered, that is, it exhibited positive spatial autocorrelation. A spatially autocorrelated dependent variable can lead to incorrect and biased estimation results in conventional regression analysis. To account for the spatial structure and for the repeated observations over time, we estimated spatial panel regressions (Elhorst, 2014;Kopczewska, Kudła, & Walczyk, 2017).
In our model selection procedure, we started with a general spatial panel specification (Belotti, Hughes, & Piano Mortari, 2017): Where: The variable y it denotes a vector of the dependent variable (the natural logarithm of the county-level NDVI for grassland) for the ith county in year t, the coefficient τ captures the lag effect, and w ij is a matrix of spatial weights that describes the spatial neighborhood structure. We used binary contiguity weights, wherein each entry of w ij equals one if i and j are neighbors (these are labelled first-order neighbors and denote observations that share a common boundary), and all other matrix elements are zero. We used rook continuity weights for the estimations but also tested the queen contiguity (i.e., shared borders and shared vertexes), as well as the second-order rook contiguity (i.e., also including the neighbors of the first-order neighbors) weights, but the results were qualitatively very similar and are thus not reported here. The variable w ij y jt is the spatially lagged dependent variable (NDVI), and ρ indicates the spatial autoregressive parameter, which depicts the size of the effect of the spatial lag on the outcome. The variables x itk and w ij x jtk represent the explanatory variables (temperature, precipitation, radiation, and animal density) and their spatially weighted forms, respectively. The coefficient β k denotes the statistical impact of the explanatory variables, and θ k indicates their spatially lagged version. The variable μ i symbolizes the county, and γ t indicates time-fixed effects. The variable v it is the error term, and λ is the parameter of the spatially autocorrelated errors.
Finally, ϵ it is a disturbance term that is assumed to be normally distributed.
Depending To do so, we tested whether the coefficients of the spatially lagged explanatory variables could be substituted with the negative product of the explanatory variables and the spatial autoregressive parameter (θ k = − β k ρ). If so, the SDM was favoured over the SEM. We decided between the SDM and SAC by comparing the Akaike information criterion (AIC). Finally, we used the Durbin-Wu-Hausman test to determine whether fixed effects or random effects were preferred (Hausman, 1978). These selection procedures suggested that the SDM with fixed effects was the preferred model for both regions.
The dependent variable was the natural logarithm of the mean NDVI during the growing season for the grassland in a county ( and Table S3). Across the study area, 68% of regions exhibited an increasing trend in the grassland NDVI during 1982-2015, and half of the increases were statistically significant (p < .05, Figure 2b). However, approximately 32% of the area experienced decreasing trends in the NDVI, which were significant in only 9% of the area (Figure 2b). The average temperature and total precipitation during the growing season were higher in Inner Mongolia than in Mongolia ( Figure 3a). However, the temperature increase was higher in Mongolia, at 0.56 C every 10 years (p < .01), than in Inner Mongolia, at 0.43 C every 10 years (p = .01; Table S3). In both regions, the observed temperature increase was substantially larger than the global average (0.06 C every 10 years) from 1880 to 2012 and the increase of 0.1 C every 10 years in recent decades (IPCC, 2013). Over the entire study period, the annual precipitation during the growing season decreased by 18 mm every 10 years in Mongolia (p < .01) and

Regions with increasing NDVI are observed in eastern
by 12 mm every 10 years in Inner Mongolia (p = .15; Figure 3b and Table S3). The average annual incoming shortwave radiation during the growing season was on average higher in Mongolia than in Inner Mongolia. The radiation increased in Mongolia at a rate of 0.07 MJ m −2 every 10 years (p = .12) and decreased in Inner Mongolia at 0.02 MJ m −2 every 10 years (p = .71; Figure 3c and Table S3). The spatial patterns of changes in the growing season temperature, precipitation, and radiation across the study area are presented in Figure S4.

| Determinants of changes in NDVI
The model selection procedure (see Section 2, Table S4 and Table S5) suggested that the spatial Durbin model with fixed effects was the The vegetation greening observed across northern latitudes has arguably been caused by higher photosynthetic activity due to rising temperatures during the growing season . Rising temperatures have also led to the lengthening of the growing season.
In addition, cropland expansion, agricultural intensification, and increasing levels of atmospheric CO 2 and nitrogen have increased the maximum rates of photosynthetic activity . In this research, our conservative mask of permanent grasslands derived from 15 years of MODIS-based land cover data (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) largely eliminated the contribution of changes in the NDVI related to land cover changes on the greening of the grasslands on the Mongolian Plateau, such as those exerted through cropland expansion onto former grasslands. We also fixed the start date and end date of the growing season in the analysis; thus, the results cannot reflect the effect of lengthening of the growing season (Miao et al., 2017). The positive effects of precipitation on the NDVI in the study region are the dominating climate factors contributing to vegetation greening. In both Inner Mongolia and Mongolia, surprisingly, the animal density was found to have a positive effect on the grassland greenness.
A priori expectation was that increasing numbers of animalsmostly ruminant grazers-would exert negative pressure on the vegetation growth on the Mongolian Plateau (Hilker et al., 2014;Mirzabaev, Ahmed, Werner, Pender, & Louhaichi, 2016) and that the negative effects of increasing grazing densities would be visible in the GIMMS NDVI time-series. To test this hypothesis, we used annual data for animal numbers of all grazers at the county level over three decades, which allowed the approximation of the footprint of the grazing livestock on the grasslands. Our measurement of the grazing intensity, which was the density of grazing livestock per grassland area, was imperfect but credible and validated spatial data concerning the distribution of animals within the counties. The animal breeds and the age distributions of the animals are, to the best of our knowledge, not available. Hence, we believe that our measurement of the county-based grazing densities represents the most fine-scale assessment of the impact of grazing on the NDVI for the Mongolian Plateau that has been conducted to date.
At the coarse scale of AVHHR, we did not observe the expected signal of livestock-induced reduction in grassland greenness due to increased animal density. While this result seems surprising at first, it corroborates the findings from grazing experiments on the Inner Mongolian steppe (Schönbach et al., 2011). In that research, within-year differences in the aboveground net primary production (ANPP) were highly correlated with the aboveground biomass and the grassland greenness. Their main explanation was that the annual variation in precipitation and grazing intensity shows a complex and nonlinear relationship with ANPP changes. Similar results have been found in other places in the world. In mixed-grass prairie in the United States (US), the ANPP was unaffected by the grazing intensity, while it was highly correlated with rainfall (Biondini, Patton, & Nyren, 1998). Also in northern Senegal, intensified grazing tended to enhance vegetation production, possibly caused by a high concentration of nutrients contained in livestock excrement (Rasmussen et al., 2018). In other places, such as the Öland steppes in Sweden and in North Dakota, US, and in the African savanna, moderate grazing has led to increasing ANPP in comparison to nongrazing and light grazing (Holdo, Holt, Coughenour, & Ritchie, 2007;Patton, Dong, Nyren, & Nyren, 2007;Van der Maarel & Titlyanova, 1989).
According to our data, higher animal densities also contributed, albeit slightly, to vegetation greening. We presume that the grazing effect on the NDVI was partly concealed by the spatial aggregation of the GIMMS NDVI data to the county level. In some pockets of intensive grazing within the counties, the negative impact of grazing on the NDVI may, nevertheless, be substantial. Subcounty data on the grazing intensity are required for analysis at finer spatial scales, but such data are, unfortunately, not available. Another reason for the small effect of the animal density on the NDVI in our analysis may be because climate factors that enhanced vegetation growth through increasing photosynthetic activity overcompensated the biomass extraction from grazing. Finally, our measure of the animal density may have overestimated the grazing pressure because we assumed that all animals grazed on grassland. In reality, to date, many livestock remain in sheds and are fed with fodder, especially in Inner Mongolia.
Regardless, these results suggest that the grassland vegetation of the region, at least that in locations with high carrying capacities due to favorable natural endowments, may currently be able to sustain larger amounts of grazing livestock, but an analysis with finer spatial details and management data would be necessary to corroborate this hypothesis.
The differences in the impactss of the covariates on the dynamics of grassland greenness between the two countries were subtle, except that the grazing intensity was higher and increased more in Inner Mongolia than in Mongolia. This result is surprising given that most of the animals in Inner Mongolia were confined to individually assigned and fenced pasture areas, which theoretically should have exacerbated land degradation because herders cannot roam with their animals to pastures with plentiful resources (Sneath, 1998). We observed a slightly larger increase in the grassland NDVI in Inner Mongolia than in Mongolia (cf. Figure 2a). Surprisingly, the diverging land-use policies of the two countries regarding pasture management were not visible in the GIMMS NDVI time-series because they either had a small negligible effect on vegetation growth in the region, as the county aggregation averaged out subcounty patterns, or the coarseresolution satellite observations failed to detect the overgrazing that is arguably occurring in the fenced pasture areas. Moreover, the fencing policy in China also includes directives to exclude some areas from grazing, and these areas are allocated for recovery-which may have positive effects on the vegetation growth (Li, Verburg, Lv, Wu, & Li, 2012;Zhang et al., 2015). Unfortunately, we could not control for this effect due to a lack of spatial data regarding these exclusion zones. However, the increasing uptake of CO 2 from the atmosphere by terrestrial ecosystems may soon be outweighed by the negative effects of climate change (Peñuelas et al., 2017). Climate-induced reduction of vegetation growth and the associated potential desertification trends may jeopardize the integrity of the grasslands on the Mongolian Plateau in the future (Mu, Li, Yang, Gang, & Chen, 2013;Sternberg, Tsolmon, Middleton, & Thomas, 2011). These impacts are particularly relevant in light of the region's arid and semiarid conditions, with adequate precipitation during the growing season becoming an increasingly decisive factor in vegetation growth (Tong et al., 2017;Zewdie, Csaplovics, & Inostroza, 2017;Zhou et al., 2015).
Using the best available historic livestock data and climate data, our results underscore that the vegetation greening on the Mongolian Plateau was mainly driven by higher precipitation in some years over the recent decade (cf. Figures 2a, 3, and 5). However, future climate change effects, including higher evaporation caused by rising temperatures and increasingly volatile precipitation, may endanger the observed greening trend.
While such insights are important for advancing science and informing policy, many questions currently remain unanswered. First, remotely sensed satellite data on greenness dynamics at a finer spatial resolution can better expose the subtle footprints that grazing may have on grassland vegetation in the region to a greater degree. Second, the effects of weather on grassland vegetation and thus on livestock farmers require more attention because climate change will likely lead to further increases in average temperatures, more temperature extremes, and possibly higher rainfall volatility, including more frequently recurring droughts. Third, controlled field experiments examining the grazing pressure on grasslands with and without fencing restrictions would permit a better understanding of the relationship between grassland greenness, climate variations, and anthropogenic impacts. Fourth, disentangling the climatic and human factors in influencing the variations of the greenness of grassland remains a challenge. The spatial panel model shows promising potentials. Our approach can also be supplemented with other alternative methods, such as the Residual Trend (RESTREND) method (Li, Wu, & Huang, 2012;Wessels et al., 2007). In view of the importance of livestock in the region for local livelihoods and the extensive environmental services that grasslands provide globally, such as sequestering carbon in soils and biodiversity conservation, more attention is required to strengthen understanding of the intertwined effects of climate change and land use on the integrity of the grasslands in global drylands.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.