Measuring grazing pressure

The disturbance indices (DI) assume that Tasseled Cap components (greenness, brightness and wetness) respond in predictable ways after a disturbance and transform these components to maximize this response. The original DI (Eq. 1, Healey et al., 2005) was developed to capture forest disturbances, which are assumed to lead to an increase in brightness and a decrease in greenness and wetness. This index was modified for grasslands (de Beurs et al., 2016), where disturbances are assumed to lead to a decrease in all three Tasseled Cap components (Eq.2). Our study region is characterized by brighter soils in steppe areas (Beznosov and Uspanov, 1960), in contrast to the soils for which the grassland DI was developed. We therefore tested a third, steppe-DI, based on the assumption that grazing leads to an increase in brightness, because herbaceous biomass is reduced and the bare cover increases (Eq.3). We assumed greenness to be higher in grazed areas due to two mutually not exclusive reasons: the first is the resprouting of green leaves, i.e. regeneration of new shoots after a disturbance, such as grazing or fire in perennial plants (Kleinebecker et al., 2011; McNaughton, 1984). This resprouting can possibly be observed in a time series of Tasseled Cap greenness in early June (Figure IV-2) and results in a slightly higher median Tasseled Cap greenness values in our spectral-temporal metrics (see Figure IV-4). Second, grazing reduces dry biomass, which otherwise may reduce the weight of greening in the overall signal.

Equation 1 Equation 2

Equation 3 Prior to calculating the DIs, we normalized all three tasseled cap components to account for the different value ranges of each component (Healey et al., 2005, de Beurs et al., 2016). As the reference for this normalization, we used locations identified as ‘ungrazed steppe’ in the field. We then used the three resulting normalized components and calculated DIHealey and DIsteppe as outlined above. As de Beurs et al. (2016) used stratified reference datasets according to climatic zones, we followed this approach and used ten zones based on the aridity index to calculate standardization datasets for the DIde Beurs.

To quantitatively rank our metrics in their ability to capture grazing pressure, we used a random forest classifier. Specifically, we used two grazing pressure classes, based on vegetation plots visited in the field: ‘heavily grazed’ (n = 190) or ‘ungrazed’ (n = 385).

These classes represent the two contrasting ends of the grazing gradient, which allowed for a binary classification. As a result, we can interpret the class membership probability of the class ‘heavily grazed’ as a proxy for grazing pressure. We used the spectral-temporal metrics derived from years when reference data were collected (i.e., 2009, 2010, 2015, and 2016) for training and cross-validation of our generalized model, as we assumed our reference data not necessarily to represent conditions outside the year when plots were visited. We trained binary (‘heavily grazed’ vs. ‘ungrazed’) random forest models for each index individually and calculated overall and class-specific accuracies (including 95 % confidence intervals around them) using a 10-fold cross-validation. For accuracy assessment we generally followed best-practice guidelines (Olofsson et al., 2014). As a second measure, we also trained a random forest model using all 27 metrics and assessed the relative importance (mean decrease in impurity) of our metrics in this model. We used the model that performed best according to a 10-fold cross-validation, which was the model with 150 estimators (i.e., the trees in the forest) without bootstrapping (i.e., all training data were used to build each tree). Based on this, we used the resulting class membership probabilities for the class ‘heavily grazed’ as an additional grazing metric (hereafter: grazing probabilities).

For the classification, we used the scikit-learn library for python, where class membership

which in turn refers to the fraction of votes for a class by each leaf of the tree. For instance, if 100 out of 150 trees in our model would estimate a probability of a pixel belonging to

‘heavily grazed’ class as 0.7, and the remaining 50 trees would estimate this probability to be 1, the resulting ‘grazing probability’ would be 0.8. We then applied this generalized, time-aggregated model to the annual layers of spectral-temporal metrics for years between 1985 and 2017 to derive annual time series of grazing probabilities, as we assumed the model trained on the data from four years to work reliably when projected on other years.

Once all metrics were calculated, we extracted metric values for all vegetation plots for the year in which the plots were surveyed. We created boxplots of metric values for the grazing pressure classes ‘heavily grazed’, ‘moderately grazed’, and ‘ungrazed’ that were originally assigned based on field data. To assess the plausibility of our grazing probability as a proxy for grazing pressure, we visualized probabilities also for an intermediate class ‘moderately grazed’, which was identified in the field but not used in our classifier. If our grazing probability was a useful proxy of grazing pressure, we would expect grazing probabilities to be intermediate for the moderately grazed class, and probabilities to lay between the two classes ‘heavily grazed’ and ‘ungrazed’. We also compared our grazing metrics using Spearman's rank correlation ρ to herbaceous biomass yield and the number of dung piles per plot counted in the field. Finally, we plotted the best-performing metrics against a range of environmental variables (i.e., soil type, soil texture, annual mean temperature, annual precipitation, mean temperature in the driest quarter, precipitation in the driest quarter, and aridity index) in order to examine whether observed differences were more likely due to differences in grazing pressure or environmental conditions.

Using the dataset of livestock concentration points, we evaluated how our grazing pressure changed away from these points over time. We created concentric buffers of 1 km, 3 km, and 6 km around livestock concentration points and derived the mean of our grazing metrics within each buffer for each year. We then compared how these values changed over time around livestock concentration points with usage intensity < 20 % (N = 2160, close to or fully abandoned, hereafter: ‘abandoned’), 20-60 % (N = 591, largely abandoned, hereafter: ‘semi-active’), and > 60 % (N = 622, still in use, hereafter: ‘active’). Please note that we purposefully used overlapping buffers, as livestock is free-ranging. Thus, the buffers around an active settlement might contain areas adjacent to an abandoned livestock station nearby, and vice versa.

We used a grazing probability threshold of 0.8 (i.e., the threshold that included 98 % of the plots characterized as ‘heavily grazed’) as a conservative measure to estimate the proportion of heavily grazed areas in our annual maps, and we applied this threshold in our accuracy assessment. As we did not have reference data for years prior to 2009, we used the same confusion matrix as for the accuracy assessment of our model, but assuming grazing probability above 0.8 (instead of a standard 0.5 threshold) to represent ‘heavily grazed’ areas. Based on this, we calculated error-adjusted area proportions (including 95 % confidence interval) of overgrazed areas on these maps following best practices (Olofsson et al. 2014).