**Land Degradation Assessment with Earth Observation**

Editor

**Elias Symeonakis**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editor* Elias Symeonakis Manchester Metropolitan University UK

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Remote Sensing* (ISSN 2072-4292) (available at: https://www.mdpi.com/journal/remotesensing/ special issues/land degradation assessment earth observation).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-4227-0 (Hbk) ISBN 978-3-0365-4228-7 (PDF)**

Cover image courtesy of Elias Symeonakis

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**



## **About the Editor**

#### **Elias Symeonakis**

Elias Symeonakis is an Associate Professor of Earth Observation and GIS and the leader of the Remote Sensing group at Manchester Metropolitan University (MMU; https://www. remote-sensing-mmu.org/). His research is related to the assessment of land degradation in drylands, mainly in sub-Saharan Africa and the Mediterranean regions. Dr Symeonakis is one of the contributing authors of the First Mediterranean Assessment Report (https://www.medecc.org/ first-mediterranean-assessment-report-mar1/) prepared by the Mediterranean Experts on Climate and environmental Change (MedECC; https://www.medecc.org/). For this work, MedECC was awarded the 2021 North–South Prize of the Council of Europe.

### *Editorial* **Land Degradation Assessment with Earth Observation**

**Elias Symeonakis**

Department of Natural Sciences, Manchester Metropolitan University, Manchester M1 5GD, UK; e.symeonakis@mmu.ac.uk

For decades now, land degradation has been identified as one of the most pressing problems facing the planet. Alarming estimates are often published by the academic community and intergovernmental organisations claiming that a third of the Earth's land surface is undergoing various degradation processes and almost half of the world's population is already residing in degraded lands. Moreover, as land degradation directly affects vegetation biophysical processes and leads to changes in ecosystem functioning, it has a knock-on effect on habitats and, therefore, on numerous species of flora and fauna that become endangered or/and extinct.

By far the most widely used approach in assessing land degradation has been to employ Earth observation (EO) data. Especially during the last decade, with technological advancements and the computational capacity of computers on the one hand, together with the availability of open access, remotely sensed data archives on the other, numerous studies dedicated in the study of the various aspects of land degradation have been undertaken. The spectral, spatial and temporal resolution of these studies varies considerably, and multiscale, multitemporal and multisensor approaches have also evolved.

This Special Issue (SI) on "Land Degradation Assessment with Earth Observation" provides 17 original research papers with a focus on land degradation in arid, semiarid and dry-subhumid areas (i.e., desertification) but also temperate rangelands, grasslands, woodlands and the humid tropics. The studies cover different spatial, spectral and temporal scales and employ a wealth of different optical, as well as radar sensors: from the finest spatial scale of an Unoccupied Aerial Vehicle (UAV), to PlanetScope, Sentinel-1 and -2, Gaofen, Landsat, MODIS, PROBA-V, SPOT VGT and AVHRR. Some of the ancillary datasets included in the methodological framework of a number of the papers are also derived from remotely sensed imagery, e.g., the SRTM digital elevation model or the Climate Hazards group Infrared Precipitation with Stations (CHIRPS) precipitation estimates. Many studies incorporate time-series analysis techniques that assess the general trend of vegetation or the timing and duration of the reduction in biological productivity brought about by land degradation (e.g., Mann—Kenndall test, Theil–Sen's slope, BFAST, TSS-RESTREND, LandTrendR). A number of papers employ statistical approaches in their analyses (e.g., principal components analysis, ordinary least squares or geographically weighted regression) or machine learning classification/regression techniques (e.g., Random Forests, Support Vector Machines). As anticipated from the latest trend in EO literature, some studies utilise the cloud computing infrastructure of Google Earth Engine to deal with the unprecedented volume of data that current methodological approaches entail.

Geographically, the papers of this SI are mostly related with areas within Africa (9 papers), which is unsurprising, as the African continent is the most severely affected by land degradation. The Asian region is also well represented with seven papers, and one paper is focused on an area in North America. In terms of the processes addressed, both abrupt and more salient changes and degradation processes are covered, with the most studied theme being the different aspects of vegetation degradation:

• Tomaszewska and Henebry [1] investigate pasture degradation in the Kyrgyz Republic, using spatiotemporal phenometrics with MODIS land surface temperature (LST) and Landsat Normalized Difference Vegetation Index (NDVI) data;

**Citation:** Symeonakis, E. Land Degradation Assessment with Earth Observation. *Remote Sens.* **2022**, *14*, 1776. https://doi.org/10.3390/ rs14081776

Received: 3 March 2022 Accepted: 2 April 2022 Published: 7 April 2022

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

**Copyright:** © 2022 by the author. 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/).


Another common theme of the Special Issue is land degradation related to drought:


Soil erosion also appears as one of the land degradation processes of interest in the Special Issue:


Land degradation related with the salinisation of the soil is also addressed in this Special Issue:

• Yu et al. [14] use Landsat data and integrate the salinization index, albedo, NDVI and the land surface soil moisture index to establish the salinized land degradation index (SDI) and apply their approach in an area that runs through Turkmenistan and Uzbekistan;

• Moussa et al. [15] compare a salinity index to an approach that employs Sentinel-2 derived NDVI time-series for detecting salt-affected soils in irrigated systems in an area of Niger.

One of the papers of the Special Issue deals with the humid tropics. In their study, Liu et al. [16] propose a framework for the improved accounting of reference levels (RLs) for the United Nations' Reducing Emissions from Deforestation and Forest Degradation (REDD+) programme. They combine the Intergovernmental Panel on Climate Change's (IPCC) Good Practice Guidance on Land Use, Land Use Change and Forestry with a land use change modelling approach and apply this to an area in southern China. Finally, the paper by Reith et al. (2021) focuses on the issue of land degradation monitoring and the methodology suggested by the United Nations Charter to Combat Desertification (UNCCD) to inform the sustainable development goal (SDG) 15.3.1 (i.e., "Proportion of degraded land over total land area"). Aiming to optimise the land degradation neutrality (LDN) efforts of the UN, Reith et al. [17] compare the land degradation assessments for an area in Central Tanzania derived using the coarser spatial-resolution MODIS and the finer Landsat data.

A clear message stemming from this SI is the ever-increasing relevance of Earth Observation technologies when it comes to assessing and monitoring land degradation. With the recently published IPCC AR6 Working Group II Report (https://www.ipcc.ch/report/ ar6/wg2/downloads/report/IPCC\_AR6\_WGII\_FinalDraft\_FullReport.pdf, accessed on 3 March 2022), informing us of the severe impacts and risks to terrestrial and freshwater ecosystems and the ecosystem services they provide, the EO scientific community has a clear obligation to step up its efforts to address any remaining gaps—some of which have been identified in this SI—in order to produce highly accurate and relevant land degradation assessment and monitoring tools.

**Funding:** Elias Symeonakis is partly funded by a LEVERHULME TRUST INTERNATIONAL ACA-DEMIC FELLOWSHIP (Contract Number IF-2021–040).

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

#### **References**


## *Article* **Remote Sensing of Pasture Degradation in the Highlands of the Kyrgyz Republic: Finer-Scale Analysis Reveals Complicating Factors**

**Monika A. Tomaszewska <sup>1</sup> and Geoffrey M. Henebry 1,2,\***


**Abstract:** Degradation in the highland pastures of the Kyrgyz Republic, a small country in Central Asia, has been reported in several studies relying on coarse spatial resolution imagery, primarily MODIS. We used the results of land surface phenology modeling at higher spatial resolution to characterize spatial and temporal patterns of phenometrics indicative of the seasonal peak in herbaceous vegetation. In particular, we explored whether proximity to villages was associated with substantial decreases in the seasonal peak values. We found that terrain features—elevation and aspect—modulated the strength of the influence of village proximity on the phenometrics. Moreover, using contrasting hotter/drier and cooler/wetter years, we discovered that the growing season weather can interact with aspect to attenuate the negative influences of dry conditions on seasonal peak values. As these multiple contingent and interactive factors that shape the land surface phenology of the highland pastures may be blurred and obscured in coarser spatial resolution imagery, we discuss some limitations with prior and recent studies of pasture degradation.

**Keywords:** Kyrgyzstan; pastures; Landsat; MODIS; land surface phenology

#### **1. Introduction**

The Kyrgyz Republic (Kyrgyzstan) is a small, highly mountainous nation of ~6.5 million in 2019, where more than 60% live in rural areas and where agropastoralism is the predominant land use [1]. Degradation of pasture resources in highly mountainous Kyrgyz Republic has been both widely reported [2–8] and disputed [9,10]. Many studies have attempted to detect land degradation in Central Asia by using remote sensing products based on finer temporal but coarser spatial resolution imagery [4,7,8,11–16]. While the finer temporal resolution products are better able to obtain clear views of the land surface, coarser spatial resolution products blend together heterogeneous surfaces. This blending or mixing is of particular concern in mountainous terrain, where differential insolation regimes arising from the interaction of aspect and slope generates microclimates that can accommodate distinct vegetation communities exhibiting different phenologies [3,17–19].

Untangling the differential influences of climate and human activity is complicated by coarser spatial resolution [7,14,20]. Pastoralism in Kyrgyzstan is characterized by transhumance, the seasonal movement of livestock to fresh pastures, and vertical transhumance in particular, where the herds are moved from low-elevation winter pastures near villages through transitional pastures in spring (and again in fall) to higher elevation summer pastures distant from settlements. This combination of spatial heterogeneity, terrain effects, and seasonality of pasture use presents challenges for detecting pasture degradation from coarser resolution remote sensing imagery. A further complication is presented by the high interannual variation in weather arising from the land-locked location of the country, its relatively high elevation, and regional climate change [12,21–24].

**Citation:** Tomaszewska, M.A.; Henebry, G.M. Remote Sensing of Pasture Degradation in the Highlands of the Kyrgyz Republic: Finer-Scale Analysis Reveals Complicating Factors. *Remote Sens.* **2021**, *13*, 3449. https://doi.org/10.3390/rs13173449

Academic Editor: Elias Symeonakis

Received: 1 July 2021 Accepted: 26 August 2021 Published: 31 August 2021

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

**Copyright:** © 2021 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/).

Here, we explore the pasture degradation question across rural Kyrgyz Republic using phenological metrics (phenometrics) derived from finer spatial resolution, but less frequent imagery (Landsat 30 m) linked to more frequent, but coarser spatial resolution products (MODIS 1 km) across 17 years (2001–2017). The phenometrics that we use here indicate the first seasonal peak in NDVI and the amount of thermal time (measured as accumulated growing degree-days calculated from land surface temperature time series) required to reach that peak NDVI. These phenometrics can be calculated for each pixel in each year for which there are sufficient, high-quality data and for which the joint time series exhibit sufficient seasonality to enable the modeling of the land surface phenology (LSP). Pasture degradation can be evaluated with phenometrics in a variety of ways, including trends, variation, extremes, and abrupt spatial and/or temporal shifts [12,20,25,26]. However, we are interested here in three questions about patterns that can complicate interpretation of patterns and trends: (1) What are the spatial patterns of the phenometrics as a function of distance from village center? (2) How do these patterns change as a function of elevation? and (3) How do these patterns change as a function of aspect?

Although we focus on the patterns of the temporal mean phenometrics during the study period, we also examine the patterns during two years with contrasting weather: hotter, drier 2007 versus cooler, wetter 2009. This study leverages the findings from prior studies [25,26] using the same approach to modeling land surface phenology and advances our understanding of the spatio-temporal dynamics of vegetation in the socio-ecological landscapes of montane Central Asia, specifically in the Kyrgyz Republic.

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

#### *2.1. Study Area*

The study area focuses on pasture lands within the territory of the Kyrgyz Republic that neighbors Uzbekistan (west), Kazakhstan (north), China (east and southeast), and Tajikistan (southwest) (Figure 1). The total area of the country is shy of 200,000 km2 (96% in land), and the 2019 population was about 6.5 million, according to the World Bank (https://data.worldbank.org/country/kyrgyz-republic; accessed on 20 June 2021). It is a highly mountainous country, where more than 56% of the territory lies above 2500 m and where the mountain ranges of the Tien Shan, Pamir, and Alatau cover more than 90% of the total land area [27]. Pastoral rangelands constitute 87% [1] of the agricultural lands in the Kyrgyz Republic. Less than 10% of the land is used for crops, while forests cover only about 5%. Our study period extends from 2001 through 2017. The Kyrgyz Republic is divided into seven provinces (oblasts)—Talas, Chuy (including the capital city, Bishkek), and Issyk-Kul in the northern part as well as Jalal-Abad, Naryn, Osh, and Batken in the southern part—and 40 districts (rayons).

#### *2.2. Geospatial Data*

The land surface phenology metrics (or phenometrics) and terrain information used for this study were supplied from our previous work [25], where a detailed description of data, its processing, and phenometrics calculation methodology can be found. We provide an overview here.

To calculate LSP metrics, we used two products: MODIS land surface temperature (LST) and Landsat surface reflectance NDVI.

**Figure 1.** Pasture land use area (light tan) and selected settlement points (blue-green) in the Kyrgyz Republic (from [28,29]) draped over the SRTM 30 m DEM [30] (Projected coordinate system: Albers Conic Equal Area). Province (oblast) names appear in yellow.

We downloaded two tiles (h23v04 and h23v05) of 8-day MODIS Terra and MODIS Aqua Land Surface Temperature (MOD11A2/MYD11A2 V006) products at 1 km spatial resolution [31] from 2001 (MODIS/Terra) and from 2002 (MODIS/Aqua) up to the end of 2017. We merged tiles, removed poor quality pixels, converted units from Kelvin to ◦C, and reprojected data to Albers Conic Equal Area at 30 m spatial resolution using bilinear resampling.

The surface reflectance NDVI dataset from Landsat Collection 1 Tier 1 Level-1 Precision and Terrain (L1TP) of Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI) was acquired for years 2001 to the end of 2017 from the USGS Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) On-Demand Interface (https:// espa.cr.usgs.gov/). We downloaded 13,285 images across 33 unique tiles (WRS-2 Paths 147–155 and Rows 30–33) which were already projected into Albers Conic Equal Area. We masked poor-quality pixels and applied an inter-calibration equation to adjust Landsat 5 TM surface NDVI and Landsat 7 ETM+ surface NDVI to the surface Landsat 8 OLI NDVI, which on average was shown to have higher values (cf., Table 3 in [32], Surface NDVI from OLI = 0.0235 + 0.9723 ETM+). Because of the small differences between the Landsat 5 TM and Landsat 7 ETM+ data [33,34], we used the same equation for both datasets.

For the analyses conducted in this study, we additionally used three other geospatial datasets: (1) digital elevation model, (2) pasture land-use mask, and (3) point coverage of settlements. We downloaded 133 tiles of SRTMGL1, the NASA Shuttle Radar Topography Mission Global 1 arc second (~30 m) V003 elevation product [30] from USGS Earth Explorer (https://earthexplorer.usgs.gov/). Tiles were merged and reprojected into the Albers Conic Equal Area at 30 m spatial resolution using bilinear resampling. We then generated aspect and slope layers. For this study, we masked out pixels from all layers where slope was greater than 30 degrees. Additionally, we created sets of three terrain layers where we left only pixels on the contrasting aspects (northern: NW-N-NE-E, southern: SE-S-SW-W) at four elevation ranges: 1800–2400 m, 2400–2900 m, 2900–3400 m, and 3400–4000 m, and elevation−aspect interactions.

The pasture land-use class (122,405 km2) was obtained from a Soviet-era land use map that was updated in 2008 using Landsat 7 ETM+ and MODIS datasets for the CACILM project [28,29]. We used those data to mask out non-pasture pixels in the terrain and LSP layers to focus only on the pasture land-use pixels.

The settlement point layer was also obtained from the CACILM project where the dataset was collected and consolidated by ECONET WWF Project [35] national teams. Over Kyrgyzstan, the layer has 703 points representing detailed administration levels, including small villages. We conducted a quality check on this dataset and eliminated points that were clearly mislocated, yielding 617 point locations (Table 1) for analysis.


**Table 1.** Number of settlement points after dataset revision and filtering per province (oblast).

#### *2.3. Methods*

2.3.1. Land Surface Phenology

We used a downward-arching convex quadratic (CxQ) function to characterize LSP [36–38]. The model uses a vegetation index—here the NDVI from a Landsat surface reflectance time-series—as proxy for active green vegetation, and thermal time—here accumulated growing degree-days (AGDD) from MODIS LST—as proxy for insolation.

To obtain AGDD, we first transformed two diurnal and nocturnal observations from the MODIS on Terra and on Aqua into a mean LST using the following Equation (1):

$$\text{mean LSTt} = \left[ \text{max(LSTt}\_{\text{TERRA1030}} \text{LSTt}\_{\text{AQULA1330}} \right] + \min(\text{LSTt}\_{\text{TERRA2230}} \text{LSTt}\_{\text{AQULA0130}}) \right]/2 \tag{1}$$

where LSTtTERRA1030 is the LST for period t at the Terra daytime overpass, LSTtAQUA1330 is the LSTt at the Aqua daytime overpass, LSTtTERRA2230 is the LSTt at the Terra nighttime overpass, and LSTtAQUA0130 is the LSTt at the Aqua nighttime overpass.

We filled gaps that resulted from missing or filtering excluded pixels using the Seasonally Decomposed Missing Value Imputation method [39], replaced all negative values with 0 ◦C, and calculated growing degree-days GDD (2) at compositing period t as the maximum of mean LST and Tbase, where Tbase was set to 0 ◦C [37,40].

$$\text{GDD}\_t = \max(\text{mean LSTt} - T\_{\text{base}}, 0) \tag{2}$$

For each year, we produced 46 GDD composites that were multiplied by 8 to account for the 8-day MODIS product composite period and accumulated each year into time series of AGDD (3), with an annual reset in January to 0 ◦C:

$$\text{AGDD}\_{\text{t}} = \text{AGDD}\_{\text{t}\cdot 1} + (\text{GDD}\_{\text{t}} \times 8) \tag{3}$$

Prior the CxQ LSP model fitting, it was necessary to further clean and filter NDVI and AGDD datasets to reduce noise and spurious data. Therefore, we filtered out observations with NDVI < 0.1 and AGDD < 100 to avoid non-vegetated or snow-covered pixels. To account for cloud contamination that may have slipped through the masking process, we looked for unusual abrupt dips in the NDVI time-series. We first calculated the simple average of NDVI observations on either side of the focal observation. We then calculated the percentage difference between the average NDVI and the focal observation and excluded observations that were ≥15% than the average of the two neighboring observations [25,26]. Having NDVI and AGDD datasets prepared for each pixel and each year (from 2001 to 2017), we applied the CxQ LSP model shown in (4):

$$\text{ANDVI} = \alpha + \beta \times \text{AGDD} + \gamma \times \text{AGDD}^2 \tag{4}$$

Using the fitted coefficients (4) for the intercept (α), slope (β), and quadratic (γ) parameters, we calculated two LSP phenometrics: Peak Height [PH=<sup>α</sup> − (β2/4 × <sup>γ</sup>)], which is the maximum modeled NDVI; and Thermal Time to Peak [TTP = −β/2 × γ], which is the quantity of AGDD required to reach PH and corresponds to the duration of green-up phase. To control CxQ LSP model fitting performance, we used a suite of six quality criteria: (i) the fitted quadratic parameter coefficient was less than zero (γ < 0); (ii) TTP greater than the AGDD at the first observation; (iii) adjusted R<sup>2</sup> greater than 0.7; (iv) Root Mean Square Difference (RMSD) less than 0.08; (v) at least seven observations in the time series to be fit, where at least three observations were distributed before and at least three after the PH; and (vi) the PH less than or equal to 1.0.

If any criterion was not fulfilled during the fitting process, then the last observation in the time series was removed and the model fitting procedure was rerun over the reduced time series. We repeated this fitting procedure until either the fitted model passed all criteria or the length of the time series was fewer than seven observations. In the latter case, the model fit for that pixel was labeled as failed and no phenometrics were calculated for that pixel.

Next, for each pixel where model fit was successful and phenometrics were obtained, we calculated the mean values of PH and TTP across 17 years. We also highlighted PH and TTP from the years 2007 and 2009, which were drier and wetter weather conditions, respectively. As mentioned in Section 2.2, we used the pasture land-use layer to mask data.

#### 2.3.2. Settlement Ring Buffer Analyses

From the settlement point coverage, we randomly selected (ArcMap Software 10.6.0.8321, Random Generator Type: default ACM599, seed:0) points with the spatial constraint of 10 km from each other, which resulted in 293 focal points for analysis (Table 1).

Then, we created 10 ring buffers around each settlement focal point from 500 m to 5000 m distance in 500 m intervals (0–500 m, 500–1000 m, 1000–1500 m, 1500–2000 m, 2000–2500 m, 2500–3000 m, 3000–3500 m, 3500–4000 m, 4000–4500, and 4500–5000 m). The spatial constraint of 10 km ensured that the ring buffers would not intersect. Finally, we prepared nine datasets (raster stacks) at four elevation classes—(1) 1800–2400 m, (2) 2400–2900 m, (3) 2900–3400 m, and (4) 3400–4000 m—yielding 36 in total, as follows:


Once these data were prepared, we extracted pixels from each set of settlement point ring buffers across all 36 datasets. During the extraction process, we ensured that only those pixels whose center point was within the ring buffer polygon were included so the same pixel would not spill into the neighboring ring buffers. Because of the prior pixel filtering by pasture land-use map, elevation classes, slope threshold value, and the fact that LSP CxQ modeling did not always succeed (i.e., no phenometrics calculated), the

number of extracted pixels varied within each ring buffer for each dataset, and in some instances, there were no pixels to extract due to the multiple filtering steps. When there were extracted pixels, we calculated the mean values for each ring buffer, which were used for the complementary analyses of the influence of elevation, aspect, growing season weather, and distance from village on the two phenometrics.

#### **3. Results**

#### *3.1. Buffer Mean Values of the Peak NDVI as a Function of Distance and Elevation*

Figure 2 displays the spatial mean ±2SE values of the fitted NDVI Peak Height within each ring buffer calculated from the mean PH values across all years. Not surprisingly, there is considerable variation, but some patterns are evident. First, the PH decreased with increasing elevation. Second, the variation in the PH increased with elevation. Third, PH increased with distance from villages at 2400–2900 m and 2900–3400 m.

**Figure 2.** Elevational gradients in modeled peak NDVI: mean ±2SE of the ten ring buffer mean values calculated from the temporal mean Peak Height. Sequential color scheme starting from left represents four classes of elevation: 1800–2400 m, 2400–2900 m, 2900–3400 m, and 3400–4000 m.

#### *3.2. Contrasting Mean Values of Phenometrics Nearby and Far from Villages*

Figure 3 subsets the distributions to focus on either end of the spatial series: the 0–500 m ring buffer captured those pasture areas closest to villages and the 4500–5000 m ring buffer captured pastures that were far from the focal village and not intruding on the ring buffer of another village. For PH (Figure 3, left panel), there was no difference in the lowest elevation class (1800–2400 m), but the differences between PH values nearby and far from villages were strong at both 2400–2900 m and 2900–3400 m, where the PH values in the distant ring buffer were higher than in the nearby ring buffer. However, the distributions of the mean PH values were significantly different between the 0–500 m and 4500–5000 m buffers only at the 2900–3400 m elevation range, according to the Kolmogorov–Smirnov two-sample test with the Dunn-Šidák correction for post-hoc multiple comparisons [41]. Note that in the highest elevation class (3400–4000 m), there were no pasture areas in the ring buffer nearest (0–500 m) to villages. It is also evident that the variation of PH in the nearby ring buffer appeared larger at higher elevations relative to that in the lowest elevation class, but this difference may result from fewer samples at high elevations.

The Thermal Time to Peak phenometric (Figure 3, right panel) decreased with elevation, as expected, and the distant TTP ring buffer means were consistently—but not significantly—lower than the nearby values.

**Figure 3.** Contrasting values from pasture areas nearby and far from villages: mean ±2SE of the two (0–500 m in light blue, 4500–5000 m in purple) buffer mean values of temporal mean Peak Height (**left**) and mean Thermal Time to Peak (**right**) for the four elevation classes.

#### *3.3. Influence of Weather on Phenometrics*

Figure 4 shows the PH (left) and TTP (right) means for the ring buffers nearby versus far from the villages. For PH, the contrast between years and distance from villages was significant only for the 2900–3400 m elevation class (Figure 4, left). At the lowest elevation class, the nearby ring buffer means were not significantly different, but the distant means were. TTP was clearly lower during 2009, the wetter year, at all elevations, and lower at higher elevations, as expected. Evaluation by the Kolmogorov-Smirnov two-sample test (with the Dunn-Šidák correction) confirms what Figure 4 (right panel) shows: the TTP distributions were significantly different between 2007 and 2009 at every elevation class, but the difference between 0–500 m and 4500–5000 m ring buffer distributions was significant only at the 2900–3400 m elevation class and only in the drier year of 2007.

#### *3.4. Influence of Aspect on Phenometrics*

We tested for differences in the distributions of buffer means as a function of aspect (NW-N-NE-E vs. SE-S-SW-W) by distance from settlement point and elevation class using the Kolmogorov-Smirnov two-sample test with the Dunn-Šidák correction for post-hoc multiple comparisons. Northerly aspects consistently showed higher Peak Height means than in southerly aspects, but these differences were not statistically different at every distance or elevation class (Figure 5). In the lowest elevation class (1800–2400 m), the effect of aspect on PH was significant starting at the 2500–3000 m ring buffer; in the 2400–2900 m elevation class, the difference in PH between aspects appeared significant 500 m closer to the settlement point (i.e., at the 2000–2500 m ring buffer). In the two higher elevation classes, the aspect differences were not significant, but the variation about the means appears lower in the more distant ring buffers.

**Figure 4.** Differences in phenometric values (Peak Height at **left**, and Thermal Time to Peak at **right**) arising from distance and weather: mean ±2SE of the two ring buffer (0–500 m as a triangle, 4500–5000 m as a circle) means of 2007 (in red) and 2009 (in dark cyan) for the four elevation classes.

#### *3.5. Interaction of Weather, Aspect, and Distance on Phenometrics*

Figure 6 shows the distributions of phenometrics for northerly and southerly aspect slopes from pasture areas nearby (0–500 m) and far from (4500–5000 m) villages in the drier year of 2007 and the wetter year of 2009. The Kolmogorov-Smirnov two-sample test with the Dunn-Šidák correction reveals four patterns of interest in the Peak Height distributions: (1) the aspect effect on PH is significant at the two lower elevation classes only far from villages; (2) at the far ring buffer at each elevation class between 1800 m and 3400 m, the distributions of PH means from southerly aspects in the drier year of 2007 are significantly different from the northerly aspect PH means in the wetter year of 2009; (3) at each distance and elevation, distributions of PH means are not significantly different by year for the same aspect, i.e., the distributions of southerly (or northerly) aspect PH means are not different between years; and (4) the distributions of northerly aspect PH means in the drier year of 2007 are not significantly different from the southerly aspect PH means in the wetter year of 2009. This last pattern is clearer in Figure 7 (top panels) and illustrates how the northerly aspect slopes moderate the impact of dry years on pasture LSP.

The patterns in the TTP distributions were not surprising (Figure 6, right): (1) decreasing TTP with increasing elevation; (2) lower TTP values in the cooler, wetter year of 2009; and (3) no apparent aspect effect in the TTP buffer mean distributions. The bottom panels in Figure 7 show more clearly how strongly the TTP mean distributions diverge between years.

**Figure 5.** Contrasting values from pasture aspect: mean ±2SE of the ten ring buffer mean values of temporal mean Peak Height at two contrasting aspects (northern aspects in dark purple; southern aspects in dark orange) for the four elevation classes.

**Figure 6.** Interaction of weather, aspect, and distance: distributions of the phenometrics (Peak Height (**left**) and Thermal Time to Peak (**right**) from two nearby (0–500 m) and distant (4500–5000 m) ring buffer mean values calculated from the hotter, drier year of 2007 (northern aspects [N] in red, southern aspects [S] in orange) and the cooler, wetter year of 2009 (northern aspects [N] in a dark cyan, southern aspects [S] in a light blue)) at the four elevation classes.

**Figure 7.** Differences in phenometric values (Peak Height on top and Thermal Time to Peak at bottom) arising from contrasts in weather and aspect: mean ±2SE of the two ring buffers at 0–500 m (**left**) and 4500–5000 m (**right**) from villages for hotter, drier 2007 (northern aspects [N] in a red triangle, southern aspects [S] in an orange circle) and cooler, wetter 2009 (northern aspects [N] in a dark cyan triangle, southern aspects [S] in a light blue circle) for the four elevation classes.

#### **4. Discussion**

The results of these linked analyses address the questions posed earlier. The Peak Height phenometric tends to increase in pasture areas as distance from the village increases. One interpretation of this pattern is the pasture areas near villages tend to be degraded. Winter pastures are closer to villages than either summer or transitional pastures and have been shown to exhibit significant differences in biogeochemistry and vegetation composition [5]. However, this general pattern can be modulated by the terrain, both elevation and aspect. Of the four elevation classes, significantly larger PH differences were

seen between nearby and distant buffer rings in the 2900–3400 m elevation class (Figure 3, left). TTP decreased with elevation class but showed no clear difference between nearby and distant values (Figure 3, right), due in part to the much coarser spatial resolution of the LST data.

For the same reason, there was no clear influence of aspect on TTP (Figure 6, right panels; Figure 7, bottom panels). In contrast, the influence of aspect on PH was strongest in the lower elevation classes (Figure 5, bottom panels). When distant from a settlement, the PH on a northern aspect slope during a dry year appeared very similar to the PH on a southern aspect slope during a wet year (Figure 7, top right). However, in pasture areas closest to settlements, this aspect influence vanishes (Figure 7, top left).

The advantages of our study are several. First, we integrated long time series of two independent remote sensing products through a simple biometeorological model of land surface phenology that responds to the progress of growing season temperature rather than to the mere passing of days [37,38]. Second, the CxQ model has been shown to capture the initial seasonal peak well in grassland LSP [25,37,38]. Third, the finer spatial resolution of our analysis captured the influence of terrain features on vegetation growth and development as captured by the phenometrics [25].

Four limitations of our study are key. First, even the 30 m Landsat data can miss important landscape features influencing LSP [42]. Second, the 1 km spatial resolution of the MODIS LST product is coarser than would be optimal for use in rugged terrain. However, there are no effective alternatives at this point. Third, the CxQ model captures the initial seasonal peak but not necessarily the post-peak decay, particularly if there is heavy grazing post-peak [25]. This limitation is not as serious as it may first appear: we expect clear evidence of pasture degradation to appear only after several years of overgrazing, particularly if coupled with a severe drought. Fourth, while the pasture land-use mask was critical for omitting non-pastoral land-uses, there is substantial uncertainty associated with its development and accuracy. We have no accuracy assessment associated with it, and we expect that commission error is more likely than omission error, which would have the effect of increasing variance and decreasing the significance of the patterns. Finally, a key caveat for interpretation of these results should be noted: the 5 km extent of analysis was not meant to capture the summer pastures associated with villages. Summer pastures can be 10–50 km or more from a given village, and the spatial allocation and arrangement of summer pastures can be quite complicated.

Given our findings, to what extent do they raise questions about previous remote sensing studies on widespread degradation of pasture resources? The relationship between the scale of observation and the scale of the phenomenon of interest is crucial to the understanding and interpretability of the remote sensing data [43]. It is clear from these results that terrain effects can be obscured by coarser spatial resolution data, yet there is another scale to consider as well. Prior studies have demonstrated significant effects of climate oscillation modes on LSP in Central Asia [12,44,45]. However, more recent work at higher spatial resolution was not able to discern a significant influence from climate modes because the influence on LSP was overridden by local landscape structure [26].

Another facet of the degradation monitoring problem relates to the impact of spatial heterogeneity on the characterization of LSP. An important empirical study [46] explored how spatial heterogeneity in the land surface phenology observed at finer spatial resolution influences the timing of phenophase detection when observed at coarser spatial resolution. They found, when the land cover was homogeneous, that greater than 60% of Landsat 8 OLI 30 m pixels exhibiting start of season (SOS) detections were within 1 day of the SOS detection within a VIIRS pixel of 500 m spatial resolution. In contrast, less than 20% of the 30 m SOS detections were within 1 day if the land cover within the 500 m pixel was heterogeneous. They further reported that the SOS detection timing at the coarser spatial resolution was controlled by the timing when roughly 30% of the 30 m pixels within the 500 m pixel transitioned to SOS [46]. Thus, phenophase detections at coarser spatial resolutions are biased toward the earlier components in the vegetated land surface if the

target landscape exhibits spatial heterogeneity. The study shows that the scaling effect on LSP is not resolvable with simple averaging. The mountain pastures of the Kyrgyz Republic and elsewhere in Central Asia exhibit a higher spatial heterogeneity than the croplands of central Iowa examined in [46]. Using a nonstandard MODIS product at 250 m spatial resolution, a study of pasture degradation in western Kyrgyz Republic found that pastures with a higher abundance of non-palatable vegetation exhibited later timing of peak NDVI during dry years in sub-alpine and mountain steppe ecozones but negligible impact in normal years [8]. The study is notable in that it included an extensive field component and demonstrated that higher NDVI values can accompany pasture degradation.

Several prior studies that have used MODIS data to detect trends in browning or greening in Kyrgyzstan and elsewhere in montane Central Asia may have been biased towards detecting browning trends and interpreting them as evidence of degradation due to using Collection 5 (C5) of the MODIS products. The differences between C5 and Collection 6 (C6 aka V006) are pronounced due to the loss of sensitivity of the red channel in the Terra MODIS [47,48]. Trend comparisons between the two collections have shown that significant negative trends in C5 were no longer significant in C6, and many nonsignificant trends in C5 appeared now as significant positive trends in C6 [49,50]. Two earlier studies [4,14] finding degradation in Kyrgyz pastures clearly used C5 products. Another study finding degradation [8] used data that was unlikely to have included the C6 corrections. A recent study [7] does not state which collection was used, but as the analysis period extended only until 2014, it is likely that it used C5 instead of C6 as well. The key point here is that as remote sensing products undergo upgrades that significantly change observed patterns, it is important to revisit those studies relying on earlier product versions to re-evaluate their findings in light of the new information.

A further concern is the use of simple linear slopes to detect significant trends. Applying ordinary least squares regression to a time series to fit a slope and declare the slope to be the trend has serious limitations from a statistical standpoint [51,52]. Positive autocorrelation reduces residual variance and inflates significance, increasing the risk of a Type I inferential error. Alternatives exist, such as the nonparametric Seasonal Kendall trend test, but it is not resistant to strong interannual autocorrelation. Fitting time series with autoregressive models has not been common in remote sensing studies.

Our findings suggest potential degradation in pasture areas as indicated by lower PHs closer to villages, but our analysis cannot rule out functional degradation arising from an increase in the coverage of nonpalatable species that can enhance NDVI without providing accessible forage [8,13,19]. Furthermore, what constitutes degradation in a particular socio-ecological system is frequently more subjective and culture-bound than is typically acknowledged [13,53–55].

#### **5. Conclusions**

We analyzed the landscape patterns of phenometrics based on fitted parameter coefficients of the land surface phenology model applied to 17 years of Landsat and MODIS seasonal time series across the montane pastures of the Kyrgyz Republic. We found that the Peak Height of NDVI generally decreased closer to villages, but the patterns were modulated—sometimes strongly—by elevation, aspect, and growing season weather. These findings raise questions about reports of pasture degradation based on coarser spatial resolution image time series. Our goal here was not to differentiate the relative contributions of these factors, because they are not susceptible to linear "unmixing". Rather, we sought to recognize their potential influence on the mixed signal observed at coarser spatial resolutions and to urge caution in interpreting, for instance, declines in NDVI trends with pasture degradation and increases with pasture remediation. The situation is more complicated. Due to the spatial heterogeneous distribution of pastures and pasture usage in mountainous landscapes, contextual information should be used to interpret remotely sensed patterns and trends in an appropriately nuanced manner.

**Author Contributions:** Conceptualization, G.M.H. and M.A.T.; methodology, G.M.H. and M.A.T.; software, M.A.T.; resources, G.M.H.; data curation, M.A.T.; writing—original draft preparation, G.M.H. and M.A.T.; writing—review, editing, and revision, G.M.H. and M.A.T.; visualization, M.A.T.; supervision, G.M.H.; project administration, G.M.H.; funding acquisition, G.M.H. Both authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the NASA LCLUC program, grant number 80NSSC20K0411. The APC was funded by the same grant.

**Data Availability Statement:** Part of the data used in this work were supplied from our previous work [25] and can be found in the NASA Land Use/Land Cover project repository: https://lcluc.umd. edu/metadatafiles/LCLUC-2015-PI-Henebry/RSE/, accessed on 20 June 2021. The remaining data presented in this study are available on request from the corresponding author (henebryg@msu.edu).

**Acknowledgments:** We appreciate the assistance and data sharing from Kamilya Kelgenbaeva. We also appreciate discussions with Munavar Zhumanova about pasture degradation. We would like to thank all reviewers and editors for their time and effort to improve the clarity, precision, and relevance of this publication.

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

#### **References**


## *Article* **Mapping of** *Kobresia pygmaea* **Community Based on Umanned Aerial Vehicle Technology and Gaofen Remote Sensing Data in Alpine Meadow Grassland: A Case Study in Eastern of Qinghai–Tibetan Plateau**

**Baoping Meng 1,2, Zhigui Yang 1,2, Hongyan Yu 3, Yu Qin 4, Yi Sun 1,2, Jianguo Zhang 1,2, Jianjun Chen 5, Zhiwei Wang 6, Wei Zhang 4, Meng Li 1,2, Yanyan Lv 1,2 and Shuhua Yi 1,2,\***


**Abstract:** The *Kobresia pygmaea* (KP) community is a key succession stage of alpine meadow degradation on the Qinghai–Tibet Plateau (QTP). However, most of the grassland classification and mapping studies have been performed at the grassland type level. The spatial distribution and impact factors of KP on the QTP are still unclear. In this study, field measurements of the grassland vegetation community in the eastern part of the QTP (Counties of Zeku, Henan and Maqu) from 2015 to 2019 were acquired using unmanned aerial vehicle (UAV) technology. The machine learning algorithms for grassland vegetation community classification were constructed by combining Gaofen satellite images and topographic indices. Then, the spatial distribution of KP community was mapped. The results showed that: (1) For all field observed sites, the alpine meadow vegetation communities demonstrated a considerable spatial heterogeneity. The traditional classification methods can hardly distinguish those communities due to the high similarity of their spectral characteristics. (2) The random forest method based on the combination of satellite vegetation indices, texture feature and topographic indices exhibited the best performance in three counties, with overall accuracy and Kappa coefficient ranged from 74.06% to 83.92% and 0.65 to 0.80, respectively. (3) As a whole, the area of KP community reached 1434.07 km2, and accounted for 7.20% of the study area. We concluded that the combination of satellite remote sensing, UAV surveying and machine learning can be used for KP classification and mapping at community level.

**Keywords:** *Kobresia pygmaea* community; unmanned aerial vehicle; Gaofen satellite; spatial distribution

### **1. Introduction**

Alpine meadow is the major vegetation type on the Qinghai–Tibet Plateau (QTP), China. It is important for animal husbandry, water conservation and biodiversity conservation [1,2]. Since the 1980s, due to the dual effects of climate change and human activities, alpine meadow grassland has experienced different extents of degradation, especially in

**Citation:** Meng, B.; Yang, Z.; Yu, H.; Qin, Y.; Sun, Y.; Zhang, J.; Chen, J.; Wang, Z.; Zhang, W.; Li, M.; et al. Mapping of *Kobresia pygmaea* Community Based on Umanned Aerial Vehicle Technology and Gaofen Remote Sensing Data in Alpine Meadow Grassland: A Case Study in Eastern of Qinghai–Tibetan Plateau. *Remote Sens.* **2021**, *13*, 2483. https://doi.org/10.3390/rs13132483

Academic Editor: Elias Symeonakis

Received: 21 May 2021 Accepted: 22 June 2021 Published: 25 June 2021

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

**Copyright:** © 2021 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/).

the source region of Yellow River, which is on the eastern part of the QTP [3]. The degradation has restricted the sustainable development of animal husbandry and seriously threatened local ecological security [4]. The degradation succession stages of the alpine meadow grassland community include Poaceae, *Kobresia humilis* (KH), *Kobresia pygmaea* (KP) and black soil type (BS). The KP community is the key stage for the management of degraded grassland [5]. In the first two stages, the original community structure and function can be quickly restored under the grazing prohibition and artificial measures [6]. However, further degradation of KP community will cause irreversible degradation until the severest stage of black soil type [7]. Therefore, it is vital to map the current distribution of KP community grassland for mitigation and adaptation measures. However, previous grassland classifications have been performed at the vegetation type level, and few at the community level [8,9]. At present, the spatial distribution and impact factors of KP community on the QTP are still unclear [1,2,10]. Therefore, it is urgent to develop a method for mapping alpine meadow at community level.

Traditional grassland vegetation community samples are mainly obtained with the few field investigation, expert knowledge and literature reviews. Due to the complex distribution and dynamic of grassland vegetation communities, the field investigation cannot meet the accuracy requirement of classification [11–15]. In addition, remote sensing (RS) vegetation indices have been commonly used as classification variables, "the same object with different spectrum" or "the different object with same spectrum" have occurred frequently [16–18]. Successful classifications at the community level requires: (1) the RS images with proper temporal-spatial resolution, coverage, sensitive spectrum band; (2) massive field observations; and (3) effective classification methods.

Compared with traditional multi-spectral remote sensing (e.g., MODIS, Landsat, HJ-1A/1B), the Gao Fen 1 (GF1) and Gao Fen 6 (GF6) satellites have significant advantages in grassland resource monitoring [19]. Each of these satellites has a high resolution of 16 m (wide field view images, WFV), a relatively large detection width of 800 km and a short revisit period of two days (four days for each, two days for combination) [20]. Additionally, GF6 satellite adds the red edge band, which is beneficial to vegetation classification. Thus, it is easier to collect high quality remote sensing images at a regional scale [21].

The massive field observation is the basis of RS classification of grassland communities. However, the resolution of satellite images is insufficient to identify grassland communities, and traditional methods require large amounts of time, labor, cost and resources. In recent years, with the development of unmanned aerial vehicle (UAV) technology, the shortcoming of satellite and traditional methods in grassland resource monitoring are supplied [22–24]. On the one hand, the aerial photographs provided by UAV have high resolution, which can be used to identify the grassland vegetation community effectively [25]. On the other hand, UAV has a large observation range, which can save time and effort. Yi et al. (2017) [26] also developed a set of UAV aerial photography system with fixed-point, multi-site, collaborative observation, which can realize massive observation over large regions [27].

With the development of classification methods, machine learning algorithm has obvious advantages in RS image classification [28,29]. Based on neural network (NN), support vector machine (SVM), random forest (RF) and other machine learning algorithms, the satellite vegetation index, phrenological characteristics, image texture and topography are considered to improve the accuracy of RS classification [24,30–32]. However, RS classification in grassland mainly includes the land use type (e.g., grassland, non-grassland, woodland, etc.) [33], different biophysics characteristics (such as grassland with high, medium and low coverage) [34] and types in different climatic zones (such as class, group and type of grassland) [35].

In this study, we aimed to map the KP community over the eastern of QTP by using the combination of UAV aerial photographing, GF WFV images and machine learning algorithms. We hope this study can be helpful for guiding further mapping of the KP community over the whole QTP and provide scientific basics for restoration and management activities.

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

*2.1. Study Area*

The study area is located at the eastern of the source region of the Yellow River, including Zeku County and Henan County of Qinghai province, and Maqu County of Gansu province (Figure 1). It is one of the most important animal husbandry basis on the QTP and also an important water source conservation area in China. The study area is located at 33◦03 ~35◦33 N, 100◦33 ~102◦33 E, with elevation ranging from 2871 to 4850 m (Figure 1c). The mean annual precipitation ranges from 400~600 mm, mean annual temperature is between −2.4~2.1 ◦C It belongs to the continental plateau temperate monsoon climate. Alpine meadow is one of the main alpine grassland types, accounting for 79.67% of the whole study area. Other than alpine meadow, mountain meadow, swamp meadow and alpine steppe account for 13.22%, 1.78%, and 1.69%, respectively (Figure 1b). The growth period of grassland plants is relatively short, only about 150 days, mainly from May to September. The grasslands are mainly used for yak and sheep grazing.

**Figure 1.** Location, grassland type and topography of study area. (**a**) Location of the study area in Qinghai–Tibetan Plateau and the observation sites; dots of different color represent the vegetation communities. Poaceae: *Elymus nutans + Stipa silena + Festuca ovina* community, KH: *Kobresia humilis* community, KP: *Kobresia pygmaea* community, BS: black soil type, MM: marsh meadow, SM: shrub meadow. (**b**) Grassland type of study area; (**c**) topography of study area.

#### *2.2. Data and Preprocessing*

#### 2.2.1. Field Observation and Preprocess of Aerial Photographs

We carried out the field monitoring for vegetation communities of alpine meadow based on aerial photographs by Phantom 3 professional and Mavic 2 zoom Quad-Rotor intelligent UAVs (manufactured by DJI Innovation Industries; http://www.dji.com (accessed on 1 June 2018). According to grassland growth status and spatial representativeness, an area in the range of 250 × 250 m was selected as an observation site, and four flight routes were designed in each site, including one GRID flight way (200 × 200 m) and three BELT flight ways (40 × 40 m) (Figure 2a). The flight way of UAVs was designed by FragMAP [22], Phantom 3 professional was used to perform the GRID flight way at a height of 20 m (red dot in Figure 2a,b), Mavic 2 zoom was used to perform the BELT flight way at a height

of 2 m (green dot in Figure 2a,c). The positional accuracy of two UAVs was ±1.5 m horizontally and ±0.5 m vertically. A photograph was then taken vertically downward at each way point automatically, the photograph resolutions of GRID and BELT were 1 and 0.09 cm, and the ground coverages were 26 × 35 m and 2.57 × 3.43 m, respectively.

**Figure 2.** Strategy of field observation and data collection: (**a**) strategy of observation site; (**b**) Phantom 3 professional and (**c**) Mavic 2 zoom Quad-Rotor intelligent UAVs.

To better identify the vegetation species, about 9~15 aerial photographs were collected randomly by operating Mavic 2 zoom manually at a height of 0.5 m in each sample site. The number of photographs was determined by the uniformity of community growth status. These aerial photographs could clearly identify plant species, which was corresponding to the traditional ground observation quadrat (Supplementary Figure S1).

According to the dominant species of grass vegetation, grassland coverage, texture features and plateau pika (*Ochotona curzoniae*, hereafter pika) activities, the aerial photographs were divided into six types, including four alpine meadow vegetation communities of Poaceae, KH, KP and BS (Figure 3 and Table 1), two land covers of shrub meadow (SM) and marsh meadow (MM). Additionally, the forest and others (bare land, construction use and waters) were acquired based on the Google Earth images and GF WFV images. Field observation was carried out at the peak time of grassland growth, and 751 sample sites were observed from 2015 to 2019 in total (Figure 1c). About 30 sample sites were acquired for forest and others.

#### 2.2.2. Region of Interest Construction

According to the GPS information recorded in FragMAP and stored in aerial photographs property files, the names of photographs were renamed by the number of 1 to 16 by the DJI Locator software [22] in each site. Then the region of interest was built based on photograph location information in the same observation site in ArcGIS and ENVI software (Figure 3d,h,l,p). Additionally, about 30 samples (region of interest, ROI) for forest and others (the water, bare land, and construction land) were selected in ENVI software, according to the GF WFV images and Google Earth images.

**Figure 3.** Classification criteria for aerial photographs of alpine meadow vegetation communities. (**a**,**e**,**i**,**m**) were aerial photographs taken at the height of 20 m by Phantom 3. (**b**,**c**,**f**,**g**,**j**,**k**,**n**,**o**) were aerial photographs taken at a height of 2 m by Mavic 2. (**c**,**g**,**k**,**o**) acquired with 2× wide-angle zoom lenses; (**d**,**h**,**l**,**p**) were Gaofen images of four types of vegetation communities corresponded; Poaceae, KH, KP and BS represented communities of *Elymus nutans + Stipa silena + Festuca ovina*, *Kobresia humilis*, *Kobresia pygmaea* and black soil type, respectively.

**Table 1.** Characteristics of vegetation communities in alpine meadow.


Poaceae, KH, KP and BS represented *Elymus nutans + Stipa silena + Festuca ovina*, Kobresia humilis, *Kobresia pygmaea* and black soil type, respectively.

#### 2.2.3. Acquisition and Preprocessing of Remote Sensing Data

The remote sensing data, including GF1 and GF6 WFV imager images, were downloaded from the China Centre for Resources Satellite Data and Application (http://www. cresda.com/EN/ (accessed on 20 September 2019)). The WFV imager was carried by GF1 and GF2 satellites, with four multi-spectral bands (800 km of swath width) and eight multi-spectral bands (850 km of swath width), respectively. The resolution of WFV image was 16 m, and the revisit period for each satellite was 4 days. Together, the revisit period could be reached up to 2 days (Table 2). Three scenes of WFV images with no cloud cover in Zeku, Henan and Maqu, during the peak of grassland growth of 2019 and 2020 were downloaded (Table 3). The GF WFV data were preprocessed using ENVI 5.3 software, and the Radiometric Calibration module, FLAASH Atmospheric Correction module and RPC (Rational Polynomial Coefficient) Orthorectification module was used for converting the original DN value to atmospheric surface reflectance, atmospheric correction and precise geometric correction of WFV images, respectively. Then, the Band Math module was used to calculate the vegetation indices of NDVI, NDWI and NDMI. The Co-occurrence measures module was used to extract image texture features of WFV images based on a sliding window with 3 × 3 pixels, and the texture indices mainly included Mean, Variance, Homogeneity, Contrast, Dissimilarity, Entropy, Second Moment and Correlation.


**Table 2.** characterization of Gao Fen (GF) wide field view (WFV) cameras.

**Table 3.** List of GF1/6 WFV images used in this study.


The DEM data were 90 m shuttle radar topography mission (SRTM) images (version V004) (http://srtm.csi.cgiar.org/ (accessed on 1 June 2018) in Geo-TIFF format. The Slope, topographic position index (TPI) and aspect were calculated based on the DEM. Then, all indices above mentioned were uniformly projected as UTM\_Zone\_47N (same as GF WFV).

#### *2.3. Vegetation Community Classification and Accuracy Evaluation*

#### 2.3.1. Classification Method

The maximum likelihood estimate (MLE), NN, SVM and RF classification methods were employed. MLE assuming each statistic of different types in every band was normally distributed, the likelihood of each pixel belonging to a certain training sample was calculated. Finally, the type of pixel was determined based on the highest likelihood [36]. NN (also called artificial neural network, ANN) referred to a multi-layer network structure, the Levenberg–Marquardt function algorithm was selected for NN training. The number of neurons and hidden layers were determined based on a trial-and-error process [37]. SVM was constructed by a set of hyperplanes in high- or infinite-dimensional space, the higher the functional margin, the lower the generalization error of the classifier. The radial basis function (RBF) was used as the kernel function, and the optimal cost and gamma values were obtained for final classification [38,39]. The RF algorithm was constructed by the classification tree, which applied a set of decision trees to improve prediction accuracy. The bootstrap sample was employed to construct a decision tree. The training samples were constantly selected to minimize the sum of the squared residuals until a complete tree was formed. Multiple decision trees were formed, and voting was used to obtain the final prediction [40,41]. MLE, NN and SVM methods were performed in ENVI supervised classification toolboxes of Maximum Likelihood Classification, Neural Net Classification and Support Vector Machine Classification, respectively. RF method was performed in ENVI Extensions toolbox of Random Forest Classification [42].

#### 2.3.2. Classification and Accuracy Evaluation

Given the classification accuracy and efficiency, three input datasets were used: (1) GF1/6 WFV spectral band (band1 to band8); (2) vegetation and texture indices (NDVI, NDWI, SAVI, Contrast, Correlation, Dissimilarity, Entropy, Homogeneity, Mean, Second moment and Variance); (3) vegetation, texture, and topography indices (DEM, Slope, Asp and TPI). About 70% of observation sites were selected randomly as a training set, and the rest were used to validate classification accuracy in each county. The standard confusion matrix was employed to evaluate the classification accuracy of images, and the overall accuracy (OA), Kappa coefficient (Kappa), user's accuracy (UA) and producer' accuracy (PA) based on the validation datasets were used to assess the precision of classification results.

#### **3. Results**

#### *3.1. Characteristics of Field Observation and Its Corresponding Multi-Indices*

The distribution of observed sites was shown in Figure 1a. The vegetation communities of alpine meadow showed a considerable spatial heterogeneity. Among the 751 observed sites, the proportion of KH community is highest, with 56.32% of all observed sites. Followed by KP community (17.04% of all observed sites), the number of KP community observation sites were 68, 37 and 22 for Maqu, Zeku and Henan, respectively. The proportion of SM, MM, BS and Poaceae only accounted for 3.33~9.85% of all observed sites.

For the four types of alpine meadow grass communities and four types of land cover, the statistics of GF1/GF6 WFV image bands, vegetation indices, topography indices and texture indices were calculated in the study area (in Supplementary Materials). The result showed that the characteristics of multi-indices in alpine meadow vegetation communities were very similar, and it was difficult to distinguish with commonly used indices (Figure 4a,b,e). Even though eight land covers could be coarsely distinguished between each other in red edge bands (band5 and band6 of WFV image) and DEM, there was relatively large error in classification (with little difference in mean values and wide range in variation) (Figure 4c,d,f).

**Figure 4.** Statistical analysis results of band3 to band6 of GF1/GF6 images (**a**–**d**), NDVI (**e**) and DEM (**f**), respectively; Poaceae, KH, KP and BS represent *Elymus nutans + Stipa silena + Festuca ovina*, *Kobresia humilis*, *Kobresia pygmaea* and black soil type, respectively.

#### *3.2. Accuracy Evaluation of the Different Classification Methods*

Accuracy assessment of classification was performed with the validation samples listed in Table 4. Among the four classification methods, the RF method performed best, with the highest overall accuracy and Kappa coefficient ranged from 74.06% to 83.92% and from 0.65 to 0.80 in three counties, respectively. This was followed by the SVM method, with an overall accuracy that ranged from 69.39% to 78.53% and Kappa coefficient that ranged from 0.60 to 0.73. The accuracies of the NN and MLE method were the worst (overall accuracy ranged from 40.78% to 73.89%; Kappa coefficient ranged from 0.24 to 0.67). Among the three classifications input, in general, the MLE, SVM and RF methods based on the input data set of vegetation indices + texture + topography exhibited the best performance, followed by the spectrum and vegetation indices + texture. However, the performance of the NN method based on the above input data set showed contrary results to the MLE, SVM and RF methods.

**Table 4.** Overall accuracy and Kappa coefficient of eight land covers based on maximum likelihood estimate (MLE), neural network (NN), support vector machine (SVM) and random forest (RF) and difference input data set in County of Zeku, Henan and Maqu.


Note: overall accuracy, OA; maximum likelihood estimate, MLE; neural network, NN; support vector machine, SVM; random forest, RF.

The results of the standard confusion matrix were shown in Table 5. The PA and UA based on the RF method were highest in three counties, with 60.84% to 97.23% and 60.73% to 78.09%, respectively. The PA and UA based on other methods showed lower value, and the classification results of the KP community were easily confused with other grass communities and land cover types.

**Table 5.** Producer's accuracy and user's accuracy of *Kobresia pygmaea* community based on maximum likelihood estimate (MLE), neural network (NN), support vector machine (SVM) and random forest (RF) in Zeku, Henan and Maqu County.


Note: producer's accuracy, PA; user's accuracy, UA; maximum likelihood estimate, MLE; neural network, NN; support vector machine, SVM; random forest, RF.

#### *3.3. Distribution and Area of KP Community*

According to the vegetation community distribution map acquired by the RF method, the spatial distribution of the KP community was fragmented with large spatial heterogeneity and small area (Figure 5). Among the three counties, the distribution of the KP community was mainly located in: the north, east and around the county urban area of Zeku County (around the town of Zequ, Qiakeri and Xipusha), with an area of 445.60 km<sup>2</sup> (6.82% of Zeku County); the northeast and central part of Henan County (east of county urban area, towns of Tuoyema and Duosun, and north of Saierlong), with an area of 176.76 km<sup>2</sup> (4.48% of Henan County); the part of county urban area, towns of Oulaxiuma, Muxihe and Awancang in Maqu County, with an area of 811.70 km2 (8.59% of Maqu County). As a whole, the area of KP community reached 1434.07 km2, and accounted for 7.20% of the study area.

**Figure 5.** Distribution of *Kobresia pygmaea* community in Counties of Zeku, Henan and Maqu.

#### **4. Discussion**

#### *4.1. Influence Factors of KP Community in the Qinghai–Tibet Plateau*

Generally, the KP community builds almost closed, non-specific, golf-course like the lawn with a felty root mat. This characteristic mat not only protects soil against intensive trampling by herbivores, but also helps to cope with nutrient limitations enabling medium-term nutrient storage and increasing productivity and competitive ability of roots against leaching and other losses [43–45]. However, with browning (patchwise dieback of lawns), crack, collapse, fragmentation of KP community turf, the water budget [46], carbon cycle [47,48] and soil nutrition [44,45,49] have been significantly changed [10].

Pastoralism may have promoted the dominance of KP community and is a major driver for felty root mat formation [10]. However, the degradation of KP grassland may be caused by both human activities and climate change [9]. The mean annual precipitation in the northern and western parts of the QTP (the elevations ranged from 4400–4800 m) was less than 450 mm, with an increase of inter-annual variability towards the west [2,10]. Grassland suffered from co-limitation of summer rainfall and nutrient shortage [10,50–53]. The types of grassland were diverse, but the species richness was low [10,15]. Hence, the ten distinct plant communities were described in this area [2]. The grassland is dominated by KP community in closed lawns with covers of 98%, and companion species less than 10 [10,43].

Our study area is located at the eastern edge of the QTP (including three counties), the mean elevation is 3758 m (Figure 1c) and mean annual precipitation ≥ 450 mm. The alpine meadow in study area consists of four types of vegetation communities, including (a) the Poaceae community (*Elymus nutans + Stipa silena + Festuca ovina*), (b) the *Kobresia humilis* community, (c) the *Kobresia pygmaea* community (KP), and (d) the denuded black soil ecosystem. Those communities consist of more than 40 species, with mosaics of KP community patches and grasses, other sedges and perennial forbs growing as rosettes and cushions [54,55]. Overgrazing is the main inducing factor for grassland vegetation community variation [5,10,56], but effect of climate still cannot be eliminated. Although we have mapped the distribution of KP community, the relative contributions from climatic and anthropogenic forces require further investigation. The main effect factor can be distinguished by combining the potential distribution based on the ecological niche model [57] and realistic distribution based on remote sensing, which is very important for alpine meadow protection.

#### *4.2. Challenges and Prospects for Alpine Meadow Grass Communities Classification* 4.2.1. Field Observation

KP community plays a vital role in alpine meadow degradation succession in QTP. However, its spatial distribution is difficult to map: on the one hand, the field observation data is lacking; on the other hand, the distribution of the KP community is under a dynamic variation with different disturbances [5,43]. The massive field observation is the basis of RS classification for grassland community. Traditional grassland vegetation community samples were obtained with the few field investigation, expert knowledge and literature reviews [11–14,58]. Field observation is mainly carried out at quadrat, plot and belt transection scales [15,59,60]. Due to the complex distribution of grassland vegetation communities, the field investigation is difficult and time-consuming. Meanwhile, the expert knowledge and literature reviews cannot meet the accuracy requirement of classification, because of the subjective bias, the dynamic climate and anthropogenic activities [15,61].

In this study, the field observation was performed by UAV based on FragMAP [14]. The resolution of each aerial photo is ~0.87 cm and covers ~35 × 26 m of ground at the height of 20 m, which is close to the traditional ground observation plot [25,62]. Moreover, the UAV is efficient and easy to operate (about 15 min to finish each observation site), which provides the possibility for rapid observation in large regions [24]. Most importantly, the waypoints, once established, can be repeatedly used (the error of two flights of the same waypoint is 1–2 m, and two photos on the same waypoint from two different flights

are almost overlapped). It is suitable to monitor the dynamic variation of grassland communities in a long-term period [25,62].

Limited by the UAV control range and battery life, the size of ROI was only 250 × 250 m, and the proportion of image raster used for training classification is relatively small. Besides, most of field observation sites were located in the flat area, which was near major traffic roads. Therefore, the spatial distribution of KP community still had some uncertainty in other regions of the study area. Moreover, the vegetation communities were distinguished by manual visual interpretation, and it requires good knowledge of plant taxonomy and time-consuming. Hence, the automatic identification of vegetation community based on aerial photograph and deep learning algorithm requires further exploration.

#### 4.2.2. Classification Variables

NDVI, NDWI and SAVI have been commonly used as the classification variables for grassland classification [20,33,59]. The vertical variation of grassland vegetation is significantly changed with topographic features in the QTP [63], hence, topographical factor is an important classification basis in alpine vegetation communities classification [64]. Additionally, texture features are also essential variables in object-based classification, which usually reflect local spatial information relating to the change of image tone [16,17]. The common method in texture feature extraction is the grey level co-occurrence matrix (GLCM). The texture metric includes angular second moment, contrast correlation, entropy, homogeneity, difference, average and standard degrees [18]. Incorporating texture feature information usually enhances the recognition of "the same object with different spectrum" or "the different object with same spectrum" [16–18].

Our results showed that, the threshold range of these RS indices for identifying the alpine meadow communities are commonly confused during extraction and identification. According to the descriptive statistical value of those RS indices corresponding to the four alpine meadow grass communities, the threshold range of KH was close to KP, and that of Poaceae was close to BS among the NDVI, NDWI and SAVI (Figure 6a–c). Although four grass communities could be distinguished in topography and texture metrics, there were relatively few differences and large errors (with little difference in mean values and wide range in variation) (Figure 6d–i). Therefore, it was difficult to distinguish the alpine meadow grass communities based on single variable and simple combinations [33–35]. RS classification accuracy can be improved by combining the RS, topographic and texture indices (Tables 3 and 4).

Due to large errors in spatial quantification of some variables (such as texture indices), the classification still has some limitations and uncertainties [29]. Hence, we consider using high spatiotemporal resolution images in future research, such as the Sentinel- 2A/B satellite images, to reduce the effects of spatial heterogeneity on spectral reflectance and acquire more detailed texture features. Secondly, screening and reconstructing the remote sensing vegetation index: combining existing vegetation index, screening out indices that are more suitable for alpine meadow vegetation community classification.

#### 4.2.3. Classification Method

Limited by the low temporal-spatial resolution, few spectrum band of RS images and field observations, most of natural grassland classification were applied in land use types (such as non-grassland, grassland, woodland, etc.) [33], different biophysics characteristics (for example, grassland with high, medium and low coverage) [34] and types with different climatic zones (e.g., groups and types of grassland) [35]. The most frequently used classification methods are visual interpretation, maximum likelihood classifiers, k-nearest neighbor and decision tree classification, and so on [65–67]. With the development of classification methods, the machine learning algorithm has obvious advantages in RS image classification [28,29]. However, the previous grassland classifications have been done at the vegetation type level, and few at the community level [8,9].

**Figure 6.** Characteristics of RS indices (**a**–**c**), topography (**d**–**f**) and texture metrics (**g**–**i**) in eight types of land covers: Poaceae, KH, KP and BS represent *Elymus nutans + Stipa silena + Festuca ovina, Kobresia humilis, Kobresia pygmaea* and black soil type, respectively.

Referenced with previously classification methods [20,58,68,69], the ANN, AVM and RF were used to distinguish the alpine meadow grass communities based on RS, texture and topographic indices in the QTP. Our results demonstrated that the RF algorithm had higher overall accuracy than other algorithms by using the same training samples (with 74.06% to 83.92%). Compared with other methods, RF is a data-driven algorithm. With the increase of input dataset, classification accuracy is improved correspondingly [66,70,71]. The RF algorithm can estimate complex nonlinear relationship and all the quantitative and qualitative information distributed within the models better; thus, these models are robust and fault-tolerant [69,70]. Moreover, the input classification indices can be acquired by different multi-spectral remote sensing images, and it helps to integrate multi-source remote sensing data [66,70,71]. However, it is difficult to train the RF model effectively with a small sample dataset. RF algorithm composes a large sample decision tree, and classification is performed based on the voting results of each decision tree, thus, has a strong tolerance for data error [40,41]. Constructing decision trees consumes more time while performing random forest classification [70].

#### **5. Conclusions**

Based on the band spectral, vegetation indices, texture feature of GaoFen 1/6 wide field view images, topographic indices and UAV field observation, this study examined four classification methods and evaluated their accuracy. Our results showed that the characteristics of RS indices in alpine meadow vegetation communities were very similar, and it was difficult to distinguish the alpine meadow grass communities based on single variable or simple combinations. The KP community could be distinguished through the RF method based on combination of RS, texture and topographic indices. The spatial distribution of KP community was fragmented with large spatial heterogeneity and small area in three counties. The area was 1434.07 km2, which accounted for 7.20% of the whole study area. Our study demonstrated it was feasible to map at the community level using the satellite remote sensing, UAV surveying and machine learning methods. In future work, more detailed texture features derived from the high spatiotemporal resolution images are required to improve the grassland vegetation community classification.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/rs13132483/s1. Figure S1: Aerial photographs of alpine meadow vegetation communities. a, b, c and d were photographs taken at the height of 2 m. e, f, g and h were photographs taken at a height of 0.5 m; Poaceae, KH, KP and BS represent communities of Elymus nutans + Stipa silena + Festuca ovina, Kobresia humilis, Kobresia pygmaea and black soil type, respectively. Figure S2: Statistical analysis results of band1 to band8 of GF1/GF6 images; Poaceae, KH, KP and BS represent Elymus nutans + Stipa silena + Festuca ovina, Kobresia humilis, Kobresia pygmaea and black soil type, respectively. Figure S3: Statistical analysis results of texture indices of GF1/GF6 images; Poaceae, KH, KP and BS represent Elymus nutans + Stipa silena + Festuca ovina, Kobresia humilis, Kobresia pygmaea and black soil type, respectively. Figure S4: Statistical analysis results of vegetation and topography indices of GF1/GF6 images; Poaceae, KH, KP and BS represent Elymus nutans + Stipa silena + Festuca ovina, Kobresia humilis, Kobresia pygmaea and black soil type, respectively.

**Author Contributions:** All authors contributed significantly to this manuscript. B.M. and S.Y. designed this study; Z.Y., B.M., H.Y., Y.Q., Y.S., J.Z., J.C., Z.W., W.Z., M.L. and Y.L. were responsible for the field observation, data processing, analysis, and writing of the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key R&D Program of China (2017YFA0604801), the National Nature Science Foundation of China (42071056, 31901393, 41861016). The APC was funded by National Key R&D Program of China (2017YFA0604801).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** GF1 and GF6 WFV imager images were downloaded from the China Centre for Resources Satellite Data and Application (http://www.cresda.com/EN/, accessed on 5 May 2020); The DEM data were 90 m shuttle radar topography mission (SRTM) images (version V004) (http://srtm.csi.cgiar.org/, accessed on 5 May 2020).

**Acknowledgments:** The authors thanks the Three-River-Source National Park Administration for their help in field sampling. We also appreciate the editor's and reviewers' constructive suggestions to greatly improve the paper.

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

#### **References**


## *Article* **Assessment of Rangeland Degradation in New Mexico Using Time Series Segmentation and Residual Trend Analysis (TSS-RESTREND)**

**Melakeneh G. Gedefaw 1, Hatim M. E. Geli 2,3,\* and Temesgen Alemayehu Abera <sup>4</sup>**


**\*** Correspondence: hgeli@nmsu.edu; Tel.: +1-575-646-1640

**Abstract:** Rangelands provide significant socioeconomic and environmental benefits to humans. However, climate variability and anthropogenic drivers can negatively impact rangeland productivity. The main goal of this study was to investigate structural and productivity changes in rangeland ecosystems in New Mexico (NM), in the southwestern United States of America during the 1984–2015 period. This goal was achieved by applying the time series segmented residual trend analysis (TSS-RESTREND) method, using datasets of the normalized difference vegetation index (NDVI) from the Global Inventory Modeling and Mapping Studies and precipitation from Parameter elevation Regressions on Independent Slopes Model (PRISM), and developing an assessment framework. The results indicated that about 17.6% and 12.8% of NM experienced a decrease and an increase in productivity, respectively. More than half of the state (55.6%) had insignificant change productivity, 10.8% was classified as indeterminant, and 3.2% was considered as agriculture. A decrease in productivity was observed in 2.2%, 4.5%, and 1.7% of NM's grassland, shrubland, and ever green forest land cover classes, respectively. Significant decrease in productivity was observed in the northeastern and southeastern quadrants of NM while significant increase was observed in northwestern, southwestern, and a small portion of the southeastern quadrants. The timing of detected breakpoints coincided with some of NM's drought events as indicated by the self-calibrated Palmar Drought Severity Index as their number increased since 2000s following a similar increase in drought severity. Some breakpoints were concurrent with some fire events. The combination of these two types of disturbances can partly explain the emergence of breakpoints with degradation in productivity. Using the breakpoint assessment framework developed in this study, the observed degradation based on the TSS-RESTREND showed only 55% agreement with the Rangeland Productivity Monitoring Service (RPMS) data. There was an agreement between the TSS-RESTREND and RPMS on the occurrence of significant degradation in productivity over the grasslands and shrublands within the Arizona/NM Tablelands and in the Chihuahua Desert ecoregions, respectively. This assessment of NM's vegetation productivity is critical to support the decision-making process for rangeland management; address challenges related to the sustainability of forage supply and livestock production; conserve the biodiversity of rangelands ecosystems; and increase their resilience. Future analysis should consider the effects of rising temperatures and drought on rangeland degradation and productivity.

**Keywords:** NDVI; precipitation; drought; breakpoints and timeseries analysis; ecosystem structural change; BFAST

**Citation:** Gedefaw, M.G.; Geli, H.M.E.; Abera, T.A. Assessment of Rangeland Degradation in New Mexico Using Time Series Segmentation and Residual Trend Analysis (TSS-RESTREND). *Remote Sens.* **2021**, *13*, 1618. https://doi.org/ 10.3390/rs13091618

Academic Editor: Elias Symeonakis

Received: 11 March 2021 Accepted: 14 April 2021 Published: 21 April 2021

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

**Copyright:** © 2021 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/).

temesgen.abera@helsinki.fi

#### **1. Introduction**

Land degradation affects ecosystem productivity and threatens its capacity to sustain human, livestock, and wildlife population specially in dryland environments. Drylands that are susceptible to desertification occupy 39.7% (~5.2 billion ha) of the global terrestrial ecosystems (~13 billion ha) [1,2]. Of this, sever land degradation is prevalent in over 10–20% of the dryland ecosystems [1–3]. For these reasons, land degradation in dryland ecosystems is recognized as one of the major environmental and socioeconomic challenges that can alter ecosystem services and human wellbeing [4–7]. Thus, understanding the rate, expansion, and severity of drylands degradation has received (and will continue to receive) considerable attention, due to their pivotal role in food production and water availability for more than 2 billion people in the world [3,8–11].

The main causes of drylands degradation are principally associated with population growth, overgrazing, inappropriate land and water use practices, and climate change impacts [1,11,12]. As a result, noticeable and persistent loss of vegetation cover and biomass productivity, reduction in forage and crop production, and decline in carrying capacity of rangelands are becoming common in these ecosystem [3,13,14]. Future predictions showed that climate change impacts (i.e., increase in surface air temperature and evapotranspiration, and decrease in precipitation) are expected to worsen poverty and inequality in developing countries [15,16].

Changes in above-ground biomass in drylands on which forage production and other life securing ecosystem services depend on can be measured by the net primary production (NPP) [13,17]. Several studies used temporal changes in NPP as an indicator of land degradation [2,18–20]. Nevertheless, the measurements and estimation of land degradation have been characterized by arbitrary assumptions, qualitative and inconsistent judgment, unreliable and spurious estimations [7,21,22]. Several studies have been criticized for underestimating the extent and severity of rangeland degradation [3,8,23]. The main reasons being the climate of dryland ecosystems (i.e., low annual precipitation and its increased interannual variability) [19,24], the successive occurrence of degradation for many decades at regional or continental scales [20], and the lack of objective measurements that quantify the extent and severity of all forms of degradation [4]. These reasons further complicate the identification of changes in dryland productivity either due to human activities (grazing and cropping) or natural variability [13,17]. Hence, the assessment of drylands degradation requires techniques that encompass spatial and temporal properties that strongly adhere to the measurement principles of repetitiveness, objectivity, and consistency [21,25]. In this regard, earth observations (OE) proved to be the only feasible means for long time and large-scale monitoring of vegetation productivity [8,9,18,24].

Vegetation indices (VIs) such as the normalized difference vegetation index (NDVI) have been used as proxy to detect and quantify long-term land degradation in drylands [6,7,11,17,18,21,25]. Two main categories of methodology have been employed to assess land degradation in drylands [7,22]—the first one observes changes in relationships between climate variables (i.e., precipitation and temperature) and VIs (i.e., measure of vegetation greenness and productivity) [8,22,23]. The second category analyzes vegetation phenology to detect structural change caused by processes such as deforestation and long-term trends [22,24]. Within each category, various methods have been applied. For example, the methods in the first category aim to control the influence of precipitation and/or temperature trends to detect permanent degradation [7,23]. The rain use efficiency (RUE) as proposed by [19]—the ratio of NPP to rainfall—has been used as an indicator of degradation [7,17,20,26,27] and ecosystem functioning [18,28]. The residual trend analysis (RESTREND) method as proposed by [11] was used frequently to estimate changes in productivity by relating annual maximum NDVI and precipitation [11,13,25,29]. In the second category, the breaks for additive and seasonal and trend (BFAST) developed by [30,31], have been employed to assess change in vegetation phenology with time (e.g., [32]).

Both RESTREND and BFAST, however, have limitations in accurately detecting land degradation when the relationship between vegetation and climate breaks and due to over

sensitiveness in areas where natural climate variabilities are prevalent, respectively [5,8,22]. Over sensitiveness is represented by the detection of false breakpoints using BFAST when dryland vegetation skips phonological cycles in response to drought while the ecosystem is healthy. The time series segmented residual trend analysis (TSS-RESTREND) method was developed by [8,22,23] to address this limitation by combining BFAST and RESTREND to detect the breakpoints where the relationship between vegetation and precipitation or temperature changes [8,23]. The TSS-RESTREND leverages the ability of BFAST in detecting breakpoints in NDVI long timeseries [5,8]. The TSS-RESTREND was applied to detect land degradation in Australia and its results were qualitatively evaluated using wildfire data [8]. However, the TSS-RESTREND lacked in providing a quantitative assessment of detected breakpoints [8].

About 40% of the United States of America (USA) land is classified as drylands [33]. The largest portion of these drylands is rangelands—comprising 31% of the USA land [34] which are mostly situated in the Western USA [35]. Rangelands support livestock production particularly beef cattle and major local economies in the Great Plains, the Intermountain West, and the Southwest rely heavily on this industry [36]. Since the pre-settlement era, about 34% of the USA rangelands have been permanently modified by human activities [37]. Woody plants encroachment, invasive species, prolonged drought, low resilience and management of rangelands could be the most critical factors that can affect current and future rangeland productivity [37,38]. Further, urban development, energy extraction (i.e., oil and gas), and expansion in agricultural lands are expected to contribute to rangelands fragmentation [39]. The degree to which climate variables (e.g., interannual variability of precipitation) can cause structural changes and affect forage productivity in New Mexico's (NM) rangeland ecosystems has not been adequately studied. Structural change in an ecosystems can be described as the significant change in the pattern of organization of the ecosystem necessary for functioning [35]. A structural change of rangeland was considered as when a pixel has a significant break in its vegetation- precipitation relationships (VPR) that leads to irreversible degradation or significant increase in productivity. Moreover, quantitative assessment of rangeland degradation and its impacts on NM's food-water-energy systems is critical for the sustainable use of its resources.

The goal of this study was to identify long-term structural and productivity changes in rangeland ecosystems due to interannual climate variability based on precipitation in New Mexico during the 1982–2015 period. Taking advantage of the consistent NDVI timeseries of Global Inventory Modeling and Mapping Studies (GIMMS), the study tested the hypothesis that there is no significant change in rangeland productivity in NM between 1982–2015 by evaluating the prevalence of significant difference in NM's rangelands mean productivity before and after the years where the VPR breaks were occurred. The specific objectives were to: (1) characterize changes in rangeland productivity in terms of (a) direction of, and (b) type of structural changes in NPP as represented by NDVI using TSS-RESTREND; (2) develop an assessment framework for the identified changes in productivity; (3) use the assessment framework along with independent productivity data to identify areas affected by significant structural changes in NM's rangeland ecosystems.

#### **2. Data**

#### *2.1. Study Area*

The study area was the state of New Mexico, USA, which covers a total land area of 314,918 km2 (Figure 1). NM's climate is dominated by arid and semi-arid conditions with an average annual precipitation ranging from less than 254 mm in the southern desert to more than 500 mm in higher elevations in the northern part of the state [40]. The mean annual air temperature ranges from less than 4.4 ◦C in high mountains and valleys in the northern to 18 ◦C in the southern parts of the state.

**Figure 1.** A map showing the location of the New Mexico, the contiguous USA (red polygon), and the nine ecoregions within New Mexico [41].

NM's land cover encompasses nine ecoregions that include the Chihuahua Deserts (22%) in south, southwestern part of NM, Southwestern Tablelands (22%) and Western High Plains (10%) in central, West and in the northwest parts of NM, Arizona/New Mexico Mountains (14%), Arizona/New Mexico Plateau (19%), Madrean Archipelago and southern Rockies [41] (Figure 1). Three major types of vegetation biomes (i.e., forest, shrubland and grassland) comprise the respective ecoregions [42].

#### *2.2. Vegetation Cover*

NM's grassland biomes, particularly in the Desert Grassland Association—desert plain grassland and mixed grassland or mixed prairie—are dominated by black grama grass (*Bouteloua eriopoda*) and tobosa grass (*Hilaria mutica*) in the desert plains. Bluestem (*Andropogon scoparius*), san blue stem (*A. halli*), and Indian Grass (*Sorghastrum nutans*) are parts of mixed grassland or mixed prairie. Woodland biome can be found within a range of 1371 m to 2286 m elevation (amsl) and consists of one-seed juniper (*Juniperus monosperma*) and pinon pine (*Pinus edulis*) sometimes with oak (Quercus spp.) with an understory of grassland, forbs, and shrubs. Coniferous forest biomes can be found within a range of 2590 to 3658 m amsl in Petran subalpine and petran montane dominated by Enngelman spruce (*Picea engel- mannii*) and subalpine fir (*Abies lasiocarpa*). Petran montane forest association (2591.8 to 2896 m a.m.s.l) covers an extensive area of the state dominated Douglas fir, and white fir species. The land use/land cover data used in this analysis included the National Land Cover Dataset of 2011 [43] as well as the state's ecoregions [41] and quadrants. These datasets were overlayed with the identified breakpoints to be related to changes in productivity over these different land cover types.

#### *2.3. GIMMS NDVI*

To study temporal changes in rangelands' structure (i.e., significant changes in NPP trend), long-term NDVI timeseries was used as a proxy to NPP. Specifically, the GIMMS NDVI product based on the third generation Global Inventory Modeling and Mapping Studies NDVI (GIMM NDVIv3.1g) data was used in this study. These data was based on the Advanced Very High Resolution Radiometer (AVHRR) sensors [5,8] and it is of the most accurate datasets currently available. This version of the data corrected for calibration errors in a previous one [22]. The GIMMS NDVI data is available for the 1981–2015 period only at spatial and temporal resolutions of 8 km and 16 days, respectively. The GIMMS NDVI dataset was obtained from ECOCAST [44].

#### *2.4. Rangeland Productivity*

Rangeland productivity for the USA was obtained from the Rangeland Production Monitoring Service (RPMS) dataset developed by the United State Forest Service [45]. The dataset was prepared using NDVI from the Thematic Mapper for the 1984–2020 period at 250-meter pixels. The dataset provides estimates of annual production of rangeland vegetation in pounds per acre, which is useful in understanding trends and variability of rangeland forage resources [45]. The RPMS dataset was coupled with the assessment framework (Section 3.2) to evaluate the significance of rangeland productivity changes compared to those identified by the TSS-RESTREND method.

#### *2.5. Precipitation, Drought, and Fire*

Gridded precipitation data developed by Parameter-elevation Regressions on Independent Slopes Model (PRISM) was obtained from PRISM Climate Group [46]. Monthly precipitation at ~4 km pixels was used. To evaluate the accuracy of breakpoints in terms of timing, extent, and direction of change relative to potential disturbances, historical drought and fire events were compared with the breakpoints. Mean monthly and annual self-calibrated Palmar Drought Severity Index (PDSI) was acquired from [47]. Fire data were obtained from New Mexico Resources Geographic Information System (RGIS) [48].

#### **3. Methods**

This analysis followed three main steps—data preparation, the application of TSS-RESTREND, and the development and application of an assessment framework to evaluate the detected breakpoints and changes in productivity (Figure 2). The first step involved data acquisition, projection, resampling, and extraction of pixel values of GIMMS NDVI, PRISM precipitation, and RPMS productivity. The second step (Section 3.1) involved the application of the TSS-RESTREND method to identify breakpoints in space and time, their significance, and their structural changes. The last step provided a framework (Section 3.2) to evaluate the detected breakpoints using the RPMS—an approach that was lacking in the work by [8,23].

#### *3.1. Characterization of Change in Productivity Using TSS-RESTREND* 3.1.1. NDVI and Precipitation Relationships

The Bimonthly GIMMS NDVI data were filtered using a quality control (QC) procedure to remove non-reliable values based on a quality flag [49,50]. Pixels with at least 75% reliable values were used. Non-vegetated pixels were excluded based on a threshold of a median NDVI of less than 0.1 [50,51]. Out of 4474 pixels, a total of 4454 pixel were used in the analysis. Complete mean monthly timeseries of NDVI (ctsNDVI) was assessed over each pixel. The PRISM precipitation data was resampled using a bilinear method to match the spatial resolution of the GIMMS NDVI and assessed at monthly time scale.

An ordinary Least Square regression (OLS) was used to develop the relationships between the ctsNDVI and the complete timeseries of optimal accumulated precipitation for each pixel for the 1982–2015 period (referred to as ctsVPR) [8]. The optimal accumulated precipitation was determined by identifying the pairs of ctsNDVI and accumulated precipitation that provide the highest correlation coefficient. The correlation coefficients of all possible OLS relationships between ctsNDVI and a matrix of precipitation accumulation period (1–12 months) with offset period (0–3 months) were evaluated following [8]. The relationships with the highest correlation coefficients were used to calculate the residual between the observed and predicted ctsNDVI (referred to as ctsVPR-residual).

**Figure 2.** Depreciation of TSS-RESTRND analysis and validation of breakpoints. In step 1, Data preparation which includes (i.e., acquisition of NPP data, precipitation, and productivity data), projection, resampling, and extraction of pixel level values. In Step II (TSS-RESTREND): method of residual fit, *p*-vector, trend of productivity change (i.e., decreasing, increasing, non-significance change and indeterminant), and ecosystem structure change were detected. In Step III (Evaluation), include stratified random sampling for validation, Welch's *t* test on randomly selected pixels (i.e., to test the significant change of productivity before and after the break years).

#### 3.1.2. Identification of Breakpoints

The BFAST method was applied on the ctsVPR-residual to list potential breakpoints. Briefly, the BFAST method decomposes the timeseries into season, trend, and reminder components—an approach that allows to detect changes in the season and trend components [52,53]. The list of potential breakpoints identified by the BFAST method [8] based on the ctsVPR-Residuals were further evaluated for their significance in the VPR—allowing to assess their impact on NPP as represented by the maximum NDVI timeseries. A Chow test was applied on the VPR-Residuals on all pixels with significant VPR (α = 0.05). Based on this test, all pixels that have no significant breakpoints in the VPR-Residual (α = 0.05) but have significant VPR (α = 0.05) were further assessed using the standard RESTREND by developing a regression between the VPR-Residuals and time (Equation (1)) [8].

$$y\_i = \beta\_0 + \beta\_1 \ge (\text{RESTREND})\tag{1}$$

where *β*<sup>0</sup> is intercept, *β*<sup>1</sup> is slope, and *x* is year

Pixels that showed significant breakpoints in VPR-Residuals (based on the Chow test above) and had significant VPR were further evaluated for the significance of these breakpoints but in the VPR also using Chow test. Pixels with significant breakpoints in VPR-Residuals but not in the VPR were further assessed using the Segmented RESTREND by developing a multivariate regression between VPR-Residual, time, and a dummy variable (Equation (2))

$$y\_i = \beta\_0 + \beta\_1 \mathbf{x}\_i + \beta\_2 z\_i + \beta\_3 \mathbf{x}\_i z\_i \text{ (Segmentted RESTRENID)} \tag{2}$$

where *β*<sup>0</sup> is intercept, *x* is year, *z* is value of dummy variable (0 or 1), *β*<sup>1</sup> is slope, *β*<sup>2</sup> is offset at the breakpoint, and *β*<sup>3</sup> is the change in the slope at the breakpoint.

Pixels that showed significant breakpoints in VPR-Residuals and in VPR may indicate the presence of significant structural changes and thus the assumption of stationarity of the accumulation and offset periods used in the calculations of optimal accumulated precipitation before and after a breakpoint [5,54]. The set if the NDVImax and precipitation timeseries before and after a breakpoint was separated to recalculate new and independent VPR on either side of the breakpoints [8]. The precipitation data were standardized to account for the differences in the accumulation and offset periods among the breakpoints (Equation (3)).

$$z\_i = \frac{(\mathbf{x}\_i - \mathbf{u})}{\delta} \tag{3}$$

where *z* is the standard score, *xi* is observed values, *μ* is the mean, and *δ* is the standard deviation. The NDVImax and standard score timeseries were further evaluated by fitting a multivariant regression (Equation (4)) [5,8].

$$\text{Yi} = \beta \mathbf{o} + \beta \mathbf{x}\_1 \mathbf{x}\_i + \beta \mathbf{z} z\_i + \beta \mathbf{y} \mathbf{x}\_i z\_i \text{ (Segmentted VPR)}\tag{4}$$

where *x* is the standardized precipitation for year *i*, *z* the value of the dummy variable (0 or 1), *β*<sup>0</sup> is intercept, *β*<sup>1</sup> is slope, *β*<sup>2</sup> is the offset at the breakpoint, and *β*<sup>3</sup> the change in the slope at the breakpoint.

Pixels that did not meet any of the above conditions were classified as indeterminant met the following conditions: (1) had no significant VPR and no significant breakpoints in VPR; or (2) no significant VPR, significant breakpoints in VPR, and no significant breakpoints in segmented VPR.

#### 3.1.3. Identification of Structural Changes

Structural changes of each pixel within NM ecosystems were identified based on three properties that include the significance of the breakpoints, direction of change (i.e., increasing and decreasing in productivity), and method of detection (as described in Section 3.1.2). With the regard to the significance level and direction of change, all pixels were classified into nine categories as shown in (Table 1) following [55] and similar to other dryland degradation studies in Australia [8] and China [5].



INC = Increase No Change, DNC = Decrease No Change, I1 = increasing trend in productivity level 1, I2 = increasing trend in productivity level 2, I3 = increasing trend in productivity level 3, D1 = Decreasing trend in productivity level 1, Decreasing trend in productivity level 2, D3 = Decreasing trend in productivity level 3, NSC = Non-significant Change.

Consequently, a pixel can be considered having:

• Non-reversable degradation or initiation of degradation if it exhibited a breakpoint with *p* < 0.05 and a negative change (Table 1) in productivity as detected by Segmented VPR or Segmented RESTREND, respectively,


A summary of the pixels based on these structural change categories was presented relative to NM's total area, quadrant, ecoregions, and major land use/land cover classes.

#### *3.2. Breakpoints Assessment Framework*

The TSS-RESTREND method by [8] provided only a qualitative assessment of breakpoints with other independent data in terms of significance and direction of changes. These two properties are important to properly characterize changes in productivity. A framework was developed in this study to address this gap by proposed means to quantitatively assess these properties. This framework consisted of four steps: (1) Develop random samples within the identified significant breakpoints; (2) Select and use independent productivity data; (3) Evaluate the random samples at the pixel level; and (4) Group and evaluate all random samples that fall within identified ecoregions—allowing to identify whether the changes in productivity at the individual pixels are reflective of consistent regional changes.

Random Samples: A set of random samples in terms of size and location can be developed following [56,57]. To estimate the size of random samples, a prior knowledge of image accuracy/variability is required [57]. The degree of variability in the RPMS data is unknown, thus it was assumed that the maximum variability would be about 50% with 95% confidence level and ±5% precision [56]. This criterion helps to determine representative sample pixels using a stratified random sampling approach [58]. The strata were developed based on ecoregions, land use/land cover, and direction of change. Based on the formula from [57], the sample size was calculated as (Equation (5)).

$$m\_{o-} \frac{Z^2 \text{ pq}}{c^2} \tag{5}$$

where *no* is the sample size, *Z* is the selected critical value of desired confidence level (1.96), *p* is the estimated proportion of an attribute that is present in the population (0.5), and *q* = 1 − *p* (0.5) with *p* and *e* (0.5) represent the desired level of precision.

The allocation of the random sample was based on the majority of the identified breakpoints in each method (Section 3.1.2). The random samples can be broadly allocated based on the direction of change and land cover within the identified significant breakpoints in the two main categories—the Segmented VPR and Segment RESTREND. Ecoregions and land use/land cover (NLCD 2011) [43] maps were overlaid to identify the locations of the pixels within these regions for further assessment.

Selection of independent data: The productivity data from the RPMS [45] was used to evaluate the accuracy of the identified breakpoints based on the TSS-RESTREND in terms of direction and significance of change in productivity. The data was resampled to match that of GIMMS NDVI [44,58]. Using the coordinates of GIMMS NDVI pixels, the corresponding RPMS productivity was extracted. The rangeland productivity was converted from pound per acre per year to kilogram per hectare per year. this process allowed to prepare a dataset that include RPMS productivity and with the corresponding TSS-RESTREND estimates of timing and direction of change, significance of breakpoints, and method of detection for each random sample.

Pixel level assessment: Using the location of the randomly selected samples (breakpoint pixels) and their identified year of breaks, the mean annual productivity for each pixel based on the RPMS data can be calculated for before and after the break years. Over each pixel, the number of years before and after the breaks were different and thus these mean values had different sample sizes. These mean values before and after the breaks over each pixel were then statistically compared using the Welch's student-test for the significance in their difference and direction of change. The Welch's *t* test was conducted assuming that the variances were not equal before and after the breaks [58,59]. The obtained results based on RPMS data using the Welch's test were compared with those from the TSS-RESTREND method. A summary of the agreement of this comparison was provided over these pixels as well as over the representative land cover classes.

Ecoregion level assessment: a similar approach was followed in this step of the framework except in this case the randomly selected pixels were grouped to represent ecoregions—all sampled pixels within each ecoregion represent a single entity (analysis unit). The mean values of productivity based on RPMS of each group before their individual breaks were further averaged and compared to the corresponding average after the break. The averages were assessed for their significance and direction of change using the Welch's student-test. The results were then compared with those obtained based on the TSS-RESTREND over the equivalent ecoregions.

#### **4. Results**

#### *4.1. Characteristics of Change*

4.1.1. Breakpoints and Direction of Change

The number of significant breakpoints detected (increased and decreased productivity) and precipitation anomalies between 1982 and 2015 were shown in Figure 3. Out of all the analyzed pixels (i.e., 4454), 814 (18.3% or ~50,000 km2) had significant breakpoints—with either increasing or decreasing trends—were detected over different land cover types. The remaining pixels that showed insignificant change (NSC), indeterminant, or identified as Agricultural represented about 55.6% (158,591 km2), 10.8% (30,656 km2) and 3.2% (9122 km2) of NM's area (~315,900 km2). The distribution of these pixels is shown in Figure 4. Out of the 814 significant breakpoints, 52.3% (26,176 km2) and 46.5% (23,232 km2) had negative (decreasing) and positive (increasing) change in productivity, respectively (Table 2). The areas that showed significant increasing trends in productivity as I1, I2 and I3 were about 10,752 km<sup>2</sup> (3.8%), 9408 km<sup>2</sup> (3.3%), and 3072 km2 (1.1%) respectively. The areas that showed decreasing trend in productivity with D1, D2, and D3 level of significance were about 9088 km<sup>2</sup> (3.2%), 13,056 km<sup>2</sup> (4.6%), and 3776 km<sup>2</sup> (1.3%), respectively.

**Table 2.** A summary of the identified direction of change in productivity based on the TSS-RESTREND method in New Mexico during the 1982–2015 period.


\* NSC = Non-Significant Change.

From all significant breakpoints, 58%, 20%, and 15.7% were observed in shrubland, grassland, and evergreen forest ecosystems, respectively (Figure 5). The highest number of detected breakpoints was 250 (or 31% out of 814) in 2005. Among them, shrublands and grasslands combined accounted for 87.6% (219), out of which 72.8% and 14.8% were over shrublands and grasslands, respectively.

**Figure 3.** A summary of the number of detected significant breakpoints along with corresponding precipitation anomaly (mm) based on a 33 year mean from PRISM. Top and bottom panels show pixels with decreasing (red line) and increasing (green line) trends in productivity, respectively.

**Figure 4.** The spatial distribution, significance, and direction of change in productivity using TSS-RESTREND (1982–2015). Bands of red (D1, D2, D3, D4) and green (I1, I2, I3, and INC) pixels indicate those with decreasing and increasing trends in productivity, respectively. NSC (yellow) and indeterminant (gray) pixels were those with non-significant and unidentified change, respectively.

**Figure 5.** The percentages of breakpoints that were identified within shrubland, grassland, evergreen forest, and other land cover types between 1982 and 2015 in NM categorized by the Segmented RESTREND and Segmented VPR methods.

#### 4.1.2. Observed Types of Structural Changes

In general, identified breakpoints by the Segmented VPR method indicate significant ecosystem structural change—either decrease (i.e., irreversible degradation) or increase in productivity. Those identified with the Segmented RESTREND method indicate changes either decrease (initiation of degradation) or increase (reversal from degradation)—in productivity that are not significant enough to alter the ecosystem structure.

From all obtained 814 significant breakpoints, 62.7% and 37.3% were identified using Segmented VPR and Segmented RESTREND methods, respectively. From those identified by Segmented RESTREND method, 20.8%, 10.5%, 3.4%, and 2.6% were over shrubland, grassland, evergreen forest, and other land cover classes, respectively. From those identified by Segmented VPR method, 37.3%, 12.3%, 9.5%, and 3.6% were over shrubland, evergreen forest, grassland, and other land cover classes, respectively.

The total number of pixels identified by the different methods and their direction of change were presented in Table 3. Those showed decreased (significant or insignificant) productivity were 786 or 17.6% (46,976 km2) out of 4454. About 4.9% (~14,016 km2) were identified using the Segmented VPR method; 4.3% (~12,160 km2) were identified using the Segmented RESTREND method; and the remaining 8.5% were identified using the RESTREND method (insignificant decrease or increase in productivity).

**Table 3.** A summary of the pixels (in % relative to 4454 pixels) and methods used to identify the direction of change in productivity in New Mexico during the 1982–2015 period.


NSC = Non-Significant Change.

From all pixels that showed increased productivity (12.8% or 570 pixels, 34,816 km2), 5.7% (16,320 km2) were identified using the Segmented VPR method (i.e., significant gradual increase), 4.6% (12,992 km2) were identified using the RESTREND method (i.e., insignificant increase), and the remaining 2.5% (7168 km2) were identified using the Segmented RESTREND method (i.e., reversal of degradation).

More than half of NM's area (55.6%~158,592 km2) show insignificant change in productivity based on the total number of pixels that fell within the NSC category (*p* ≥ 0.1). The pixels identified with the RESTREND method that accounted for 25.8% (~73,472 km2), and the remaining 29% were considered indeterminate.

Regionally, non-reversible decrease in productivity (i.e., degradation) was mostly observed in the northeastern (1.8%~5184 km2), northwestern (1.2%~3456 km2), and southeastern (0.99%~2816 km2) quadrants NM (Figure 6). Moreover, the initiation of degradation (i.e., decreasing trend based on the Segmented RESTREND method) was mostly detected in northeastern (1.8%~4992 km2), and southeastern (2.1%~5888 km2) quadrants, with a negligible initiation of degradation in the northwestern and southwestern NM (Figure 6).

**Figure 6.** Area and direction of change in productivity in northeastern, southeastern, northwestern, and southwestern quadrants of New Mexico as identified using the Segmented VPR (degraded vegetation or gradual increase) and Segmented RESTREND (reversal and initiation of degradation) methods.

From all analyzed pixels (i.e., 4454 pixels), the ones that showed increasing trends in productivity based on the Segmented VPR method (i.e., 5.7% or ~16,320 km2) were detected mostly in the southwestern (2.27%~6464 km2), southeastern (2%~5824 km2), and northwestern (1.3%~3712 km2) quadrants (Figure 6). Furthermore, 2.51% of the pixels (7168 km2) revealed reversal from degradation—increased productivity based on the Segmented RE-STREND method (Figure 6). From which, 0.45% and 0.92% were in the northeastern and southeastern quadrants, respectively (Figure 6). The largest number of pixels that showed NSC was detected in northwestern (46,144 km2), while southeastern quadrant exhibited the least (32,256 km2~11.1%). Northeastern (40,896 km2) and southwestern (39,296 km2) NM revealed an equal number of NSC pixels during the study period.

#### *4.2. Dominant Land Cover Class Changes*

Significant trends (increasing or decreasing) in productivity (irreversible degradation, initiation in degradation, reversal in degradation, and initiation in reversal of degradation) were observed on NM's dominant land cover classes that include shrubland, grassland, and evergreen forest (Figure 7).

From all analyzed pixels (i.e., 4454), 2.2% of NM's grassland (6336 km2), 4.5% of the shrubland (12,800 km2), and 1.7% of the evergreen forest pixels (4800 km2) showed significant decreasing trends in productivity–either with a complete or initiation of degradation. From all analyzed pixels, ~1.4% (3840 km2) and 0.88% (2496 km2) of NM's grasslands that showed significant decrease in productivity were attributed to initiation of degradation (i.e., the Segmented RESTREND method) and non-reversal degradation (i.e., the Segmented VPR method), respectively (Figure 7).

Similarly, from all analyzed pixels, 2.2% (6336 km2) and 2.3% (6464 km2) of NM's shrublands that showed significant decrease in productivity were attributed to initiation of degradation and non-reversal degradation, respectively. The degradation of shrublands was dominant in the northwestern (2048 km2) and southeastern (2112 km2) quadrants, while the degradation in grasslands was notably observed in northwestern (512 km2) and northeastern (1664 km2) quadrants.

**Figure 7.** Area and direction of change in productivity in grasslands, shrublands, and evergreen forests in New Mexico as identified using the Segmented VPR (degraded vegetation or gradual increase) and Segmented RESTREND (reversal and initiation of degradation) methods.

On the other hand, 5.7% (out of the 4454 pixels) of NM's shrubland pixels (16,248 km2), 1.3% of the grassland pixels (3648 km2), and 0.92% of the evergreen forest pixels (2624 km2) experienced either significant gradual increase in productivity or reversal of degradation (Figure 8). From all analyzed pixels, 1600 km<sup>2</sup> (0.56%) and 2048 km<sup>2</sup> (0.72%) of the grasslands that showed significant increase in productivity were attributed to reversal in degradation (using the Segmented RESTREND method) and significant increase in productivity (using the Segmented VPR method), respectively.

#### *4.3. Assessment of Breakpoints*

This section provides a summary of the results obtained based on the breakpoints assessment framework described in Section 3.2.

#### 4.3.1. Identified Random Samples

The total number of breakpoints based on the Segmented RESTREND and Segmented VPR methods were about 304 and 510, respectively. From which, 384 samples were randomly selected with 165 and 219 were based on the Segmented RESTREND and Segmented VPR methods, respectively (Tables A1 and A2 in Appendix A). Since the Segmented VPR had more significant breakpoints with a noticeable change in productivity, only those pixels were subjected to the random selection. Only 155 samples out of the 219 (or 71%) showed significant difference in productivity before and after the break years (either decreasing or increasing) (Figure 9). The remaining 64 samples did not meet the criteria for significance and were not considered. From the 155 samples, 24% and 76% were obtained over grasslands and shrublands, respectively. The distribution of the grassland samples (i.e., 24%) represented the Southwestern (SW) Tablelands (46%), Arizona/New Mexico (AZ/NM) Plateau (35%), and Chihuahua Desert (16%) ecoregions. Similarly, 58%, 25%, and 11% of the shrubland samples (i.e., 76%) represented the Chihuahua Desert, AZ/NM Plateau, and SW Tablelands ecoregions, respectively.

#### 4.3.2. Changes in Productivity at Pixel Level

A summary of the comparison of the significance of the differences in mean productivity before and after the breakpoints between the Segmented VPR and the RPMS data was shown in Table 4 along with the direction (increase or decrease) of change in productivity. The results were presented as the percent of pixels that fell within each category relative to the number of the random samples (i.e., 155).

**Figure 8.** Direction of change in productivity in New Mexico's quadrants during the 1982–2015 period indicating (**a**) non reversal degradation, (**b**) gradual increase, (**c**) initiation of degradation, and (**d**) initiation of reversal from degradation along with the corresponding methods used.

Based on the Segmented VPR method, all the 155 randomly selected samples indicated significant difference in productivity before and after the break years. However, based on the RPMS only 55% of them showed significant differences in mean productivity before and after the break years until 2019, respectively (Table 4). From the 155 samples, 37% and 18% showed persistent decreasing and increasing in mean productivity after the break years until 2019 on RPMS data, respectively. About 36% and 9% showed increasing and decreasing trends before and after the break years on the RPMS data, respectively.

From the 37% of the sampled pixels that exhibited persistent decrease in mean productivity after the break years, 13% were in the Chihuahua Desert, and 18% in the SW Tablelands ecoregions. From the 18% of the sampled pixels that exhibited significant and consistent increase in mean productivity after break years, AZ/NM Mountains and AZ/NM Plateau ecoregions accounted for 3% and 15%, respectively. From the 36% of the sampled pixels with increased but insignificant difference in mean productivity before and after the break years, 35% were over the Chihuahua Desert ecoregion, and the remaining were equally obtained over AZ/NM Mountains and AZ/NM Plateau (Table 4). From the 9% of the sampled pixels with decreasing trend but insignificant difference in productivity before and after the break years, 5% were over AZ/NM Plateau.

**Figure 9.** The distribution of the randomly selected samples from the identified breakpoints based on the Segmented VPR method along with the ecoregions in New Mexico.

**Table 4.** The percentages of pixels with decreasing and increasing trends in productivity as estimated from the RPMS data relative to randomly selected pixels (155) based on the Segmented VPR method categorized based on their significance of difference in mean productivity before and after the break years over major ecoregions in New Mexico.


In the sampled grasslands and shrublands pixels (i.e., 24% and 76% of the 155 samples, respectively) with either decreased or increased productivity, the difference in mean productivity before and after the break years was insignificant in 5% and 40%, respectively on RPMS (Table A3 in Appendix B). The grassland (12% of the samples) and shrubland pixels (26% of the samples) with decreased productivity had significant lower mean productivity before the break years on RPMS data. Similarly, out of the pixels that showed increased trend on productivity, 7% from grasslands and 11% from shrublands had significantly higher mean production than that before the break years on RPMS (Table A3).

#### 4.3.3. Changes in Productivity at Ecoregion Level

A summary of the comparison between the Segmented VPR and RPMS at the ecoregion level (Section 3.2) including grassland and shrubland cover classes based on the

sampled pixels (i.e., 155) was provided in Table A5 (Appendix C)—allows to highlight whether the changes at the individual pixels are reflective of those at the regional level.

Based on the RPMS data, a continuous and significant decrease in shrublands' mean productivity after the break years in the Chihuahua (Welch's test *p* ≤ 0.0001) and the AZ/NM Plateau (*p* = 0.0397) ecoregions was observed. Similarly, significant decline in mean productivity of grasslands of the SW Tablelands (*p* ≤ 0.0001) and AZ/NM (*p* = 0.019) ecoregions were observed after break years. The shrublands within the AZ/NM Mountains and the SW Tableland ecoregions and the grassland in Chihuahua and the AZ/NM Plateau ecoregions exhibited insignificant differences in mean productivity between before and after the break years—stable ecosystem productivity during the study period.

In contrary, the sampled pixels over the shrublands in the Chihuahua Desert (*p* = 0.0194) and the AZ/NM (*p* = 0.00765). Mountain ecoregions showed significant increase in mean productivity after the break years (Table A5). There was a significant increase in mean productivity in the grasslands within AZ/NM Plateau ecoregions after the break years. However, sampled shrubland and grassland pixels in the AZ/NM Plateau and the Chihuahua Desert, respectively, showed insignificant difference in mean productivity before and after the break years, suggesting a negligible increase in mean productivity after the break years during the study period.

#### **5. Discussion**

#### *5.1. Characteristics of Change*

The significance of the breakpoints as identified by the TSS-RESTREND methods can be interpreted relative to observed ecosystem structural changes [8,13]. Out of 67.9% of the pixels that met the criteria of RESTREND [13], 8.5% and 4.6% showed decreased and increased productivity, respectively. These pixels exhibited gradual change as their VPR remined consistent over time with no major ecosystem structural changes [13,60]. The behavior of the pixels that were identified by the Segmented VPR method (11.6%) experiencing irreversible degradation (4.9%) or increased productivity (5.7%) can partially be attributed to abrupt land use changes (decreeing or increasing) induced by human activities and climate variability [13,61] such as overgrazing or easing of drought conditions [62]. More than half of NM that did not experience significant change in productivity (NSC = 55.6%—Table 2) was dominantly in the western part of the state. This can partially be explained by the fact that western NM is the driest region in the state. Thus, it experiences weaker interactions related to climate variability and human activities—thus minimal effects on productivity. The pixels that were identified as indeterminant (10.8%) were those that none of the methods was able to fit their observed behavior [8,55] and there was no clear explanation for such behavior.

There were relatively higher human activities in northeastern, northwestern and southeastern NM represented by crude oil and natural gas production and livestock grazing practices. The prevalence of irreversible degradation and initiation of degradation were apparent in these regions [63]. On the other hand, significant changes in productivity (i.e., increase) was observed mostly in southwestern and northwestern NM, where forest were minimally influenced by human activities. On landscapes where human activities were dominant, i.e., southeastern NM, significant increase in productivity can be associated with parallel restoration activities and land use management practices [55,64].

Overall, there was a consistent increase in the number of detected significant breakpoints during the 1980s, 1990s, and 2000s, with 9.5%, 19.2%, and 71.4%, respectively (Figure 10). Again, the combined effects of climate change and anthropogenic factors can explain this increase [65]. For example, NM's precipitation remained close to the long-term average (only recently showed increased variability) since 1920 while air temperature showed an increasing trend since 1970s [63,66,67] with a simultaneous increase in fossil fuel production activities since 1980s.

**Figure 10.** A summary of breakpoints with (**a**) decreasing and (**b**) increasing trends in productivity along with the weighted average of self-calibrated Palmer Severity Drought Index (sc-PSDI).

It is challenging to identify a single factor that can be considered as the direct cause of these breakpoints and ecosystem structural changes [8,68]. However, it was possible to compare these breakpoints with some observed ecosystem structural changes whether gradual (e.g., land use management and climate change) or abrupt (e.g., wildfire). Such a comparison can help in evaluating the accuracy of these breakpoints in terms of timing, distribution, and direction of change [8,53].

#### *5.2. Land Cover Changes Relative to Drought and Wildfire*

Some of the major drivers of change in NM's dryland ecosystems were identified as climate (e.g., drought) which is influenced by increased concentration of atmospheric greenhouse gases (GHG) (e.g., CO2); wildfire; grazing practices (i.e., livestock density); and land use conversion [21,35,68] were used here to highlight their effects on the breakpoints.

#### 5.2.1. Detected Changes Compared to Previous Studies and Current Restoration Activities

The significant decrease in grasslands and shrublands was mostly observed in northwestern, northeastern, and southeastern NM—consistent with [4] that indicated that degradation was mostly over grassland-savanna. Degradation of grasslands, in some cases, was directly linked to increased productivity of shrublands. In NM, increased productivity in shrublands was noticed southwestern deserts and plains. This increase was attributed to land cover conversion from black grama and other valuable grasses-dominated areas to bushes due to over grazing and drought [69]. The evident climate warming, expected increase interannual precipitation variability, and projected increase in aridity can further enhance the growth of shrubs (i.e., encroachments of woody plants) over grasses since the later are heavily depended on transient surface moisture [70,71]. It was also argued that increased CO2 concentration in the atmosphere can improve plant CO2 uptake and reduce water loss through plant's stomata thus increasing the photosynthesis process and this mechanism actually favor woody vegetation over non-woody ones [72,73].

The general notion that was suggested in some studies was that encroachments of shrublands to grasslands can be considered as one stage of degradation that can threaten the integrity of the rangelands [74]. However, it was argued that, shrublands are not necessarily degraded, nor do they necessarily represent "degradation" due to their ability to support valued ecosystem services [71,75], and they also have a long-term mean annual above ground NPP—equivalent to that of grasslands [76]. In NM, there have been a number of restoration and brush management efforts undertaken to control invasive species, and retain favored shrub species [64]. This suggested that the attribution of drivers to changes in productivity in environmental studies remains challenging owing to the diversity of driving forces and limited sources of ground truth data for validation [77].

The increased productivity in grasslands that was observed in western NM can partially be attributed to local scale successful restoration efforts by the Bureau of Land Management [64]. These efforts targeted the replacement of Creosote and mesquite by healthy grasslands, and reclamation of surfaces resulted from oil and gas extraction operations. According to Powell [78], gradual improvements in range conditions (increased productivity) in southwestern NM was not only related to the results of better moisture, but also cumulative efforts on rangeland management, such as proper stocking, vastly improved grazing distribution, and brush management.

#### 5.2.2. Breakpoints and Drought

The sc-PSDI for NM and the detected significant breakpoints with decreasing and increasing trends were shown in Figure 10. A weighted average of the sc-PDSI was calculated over the breakpoints with decreasing and decreasing trends separately to evaluate their timing against drought events. Some previous findings suggested that extended periods of drought (or relatively dry conditions) can introduce lasting negative impacts (degradation) on rangelands and other ecosystems [8,79]. Based on Figure 10, it was noticed that from all the detected breakpoints during the study period (1982–2015), 67.9% were observed only in 2000s from which 38.2% exhibited decreased trend in productivity (Figure 10a). Coincidently, these breakpoints overlapped with frequent and extended periods of drought events in NM that were also observed regionally in the southwest USA during 2000s. The impact of this regional drought had resulted in a decrease in productivity in more than 30% of the coterminous USA, out of which ~15% (equivalent to more than 41 million ha) was rangelands [37].

The increased number of breakpoints after 2000 suggested a permanent degradation or damage to ecosystems which mechanistically can occur when these ecosystems have little to no time to recover from a previous consecutive drought event. A recovery period from a drought event as defined by [79] is the return of an ecosystem to pre-drought values of GPP and can vary from immediate to multiple years depending on vegetation, climate, disturbance, and drought. It was indicated that dryland ecosystems such as those in NM, had experienced increased recovery periods that were highly sensitive to precipitation and temperature [79]. With the expected future projections of rising temperature and drier conditions, the drought recovery of dry ecosystems would even be longer. When droughts have shorter return period (or more frequent) this makes ecosystems more susceptible to drought or ecosystem degradation continues and builds up until a threshold is reached and the degradation become irreversible (or permanent), a condition that was referred to in [79] as a tipping point. This study suggested that these thresholds may have been reached as represented by the detected breakpoints. This explains, to some extent, the appearance of breakpoints compared to drought events. The timing of these breakpoints suggested the accumulation of dry conditions that impact ecosystem functions until reaching tipping points. Further analysis is needed to understand the timing or emergence of the breakpoints and drought accumulation periods and the response of ecosystems relative to vegetation type and climate.

Moreover, these recent exceptionally drought conditions in NM (during 2000–2015) also varied in spatial extent, duration, and intensity [80]. These variations may have contributed to the occurrence of 38% breakpoints with decreasing trends (i.e., initiation and irreversible degradation) in southeastern (Eddy, Lea, Chaves, and Otero counties); northeastern (Colfax, De Bacca, and San Miguel counties); and northwestern NM (San Juan county). From all pixels that experienced significant decrease in productivity, 35% were

observed in De bacca, San Miguel, San Juan, and Otero counties as drought might have a profound influence in reducing vegetation productivity (e.g., plant mortality) [81,82].

On the other hand, reversal from degradation and significant increase in productivity was dominantly observed in the northwestern, southeastern, and southwestern NM (Figure 8). In NM, wetter years were observed from 1986 to the end of 1990s, followed by drier years from 1999–2003 and 2005 (Figure 10). Significant breakpoints with increase in productivity were identified during 2000s dry years in Otero, Dona Anna, Luna, San Juan and Socorro counties. About 13.5% of these breakpoints were observed in the first four of counties. This can partially be attributed to the fact that precipitation remained near the long-term mean after the drier years over these counties that resulted in reversal from degradation or significant increase in productivity, respectively (Figure 10).

#### 5.2.3. Breakpoints and Wildfire

Wildfire generally reduces plant cover, alters habitat structures, decrease rangeland conditions, and requires much longer recovery period [83]. The significant breakpoints with decreased productivity (i.e., irreversible degradation based on Segmented VPR) that coincided with fire events are shown in Figures 11 and 12. Major wildfire incidents occurred in 1989 and 1994 followed by a consistent increase in frequency each year since 2001. This increase was noticeably concurrent with drought events. The number and time of some breakpoints mostly followed these of the fire incidents. The alignment of fire incidents with breakpoints can indicate the ability of TSS-RESTREND to detect the timing of ecosystem changes as suggested by [8]. However, based on Figures 11 and 12, it appeared that only a small number of breakpoints coincided with fire incidents during the study period—suggested that fire may have a limited contribution to the development of breakpoints. While, the timing of these fire incidents can partly explain the occurrence of the breakpoints, it was not rationale to state that these fire incidents were the only direct cause of the breakpoints. Other factors need to be considered to provide a rational explanation of the remaining breakpoints such as vegetation types and their response to environmental and climate disturbance (drought). The accumulation of successive dry conditions can push an ecosystem to pass a threshold of permanent damage of vegetation ideal fire-prone conditions. Further analysis is needed to identify the causes of breakpoints relative the fire and was out of the scope of this analysis.

**Figure 11.** Timeseries of breakpoints concurrent different types of fire incidents in New Mexico during the 1984–2015 period.

#### *5.3. Breakpoints and the RPMS*

Based on the RPMS data, a decrease in the mean annual productivity was detected over the pixels with significant breakpoints (i.e., identified by Segmented VPR method) during the 1984–2019 period (Figure 13). The mean productivity of these pixels was lower than that of the long-term of the 36 years which was about 613.2 kg/ha. Moreover, the variability in the annual productivity of these pixels ranged from 588 kg/ha in 1995 to 641 kg/ha in 2006.

**Figure 12.** A summary of breakpoints and fire incidents indicated with the timing of some major incidents.

**Figure 13.** Productivity anomaly based on the Range land Production Monitoring Service (RPMS) [45] averaged over the pixels that were identified with the Segmented VPR method with decreasing trend during the 1984–2019 period.

Of the sampled pixels that were evaluated for changes in productivity before and after the break years, 38% with decreasing trend showed a significant difference in mean productivity after the break years with 12% in grasslands (Table A3) and 26% in shrublands (Table A4). Similar to this finding, [84,85] indicated that there was an apparent decline in NPP in the southwestern USA that was attributed to the response of these ecosystems to the combination of a warmer temperature and a decline in precipitation. The pixels that showed significant decrease in productivity in NM's rangeland (i.e., in Chihuahua Desert and AZ/NM plateau ecoregions) coincided with those from [45] that showed a similar behavior (Table A5). The difference in the productivity before and after the break years might be substantiated by the xeric nature of the southwestern USA and southeastern Great Plains in response to interannual variability of the precipitation [86]. On the other hand, 8% and 36% of sampled pixels with decreased and increased productivity (i.e., based on the Segmented VPR), respectively showed insignificant difference in mean productivity before and after the break years based on the RPMS data [45]. These findings were in agreement with those by [83] which suggested that the structure and composition of semiarid and arid regions of the southwestern USA have undergone noticeable changes over the last two decades.

#### **6. Limitation and Future Work**

As the study aimed to detect degradation of vegetation, rangelands (i.e., shrubland and grasslands), the authors acknowledged some limitations that need to be considered when using these findings. The TSS-RESTREND was able to identify the breakpoints, their location, time, and type and direction of change in productivity at the state level. Due to the limited availability of ground biomass data, pixels that exhibited significant changes were validated using derived productivity from the RPMS [45] that showed only 55% agreement. Ground-based data can be more accurate means in assessing the identified breakpoints. Another important factor that can contribute to degradation is the temperature which can significantly affect NM's dryland ecosystems.

For future, the authors would address some of these limitations using ground observations (when available) for example from the long-term rangeland monitoring observatories at the Jornada Experimental Range. The obtained breakpoints seemed consistent with some of the disturbances (drought and fire). However, the authors would consider the combined use of temperature and precipitation to improve the detection of breakpoints such as the pixels that have been classified as indeterminant. Additional quantitative analysis would be conducted to determine how drought may affect the timing of breakpoints, rangeland condition, and productivity, which can be useful in developing rangeland management practices to adapt and mitigate drought impacts.

#### **7. Conclusions**

This study evaluated the degradation of New Mexico's rangelands during the 1984 2015 period with respect to climate using an NDVI timeseries—as a surrogate of NPP and precipitation to represent climate variability. These datasets were evaluated using the TSS-RESTREN method to detect breakpoints in the NDVI timeseries, direction and significance of change in productivity at each pixel. The study developed a breakpoint assessment framework that allowed to quantitively evaluate the identified changes that used an independent productivity data (i.e., RPMS). A qualitative assessment of the changes against land management activities, drought, and wildfire was also conducted.

The study indicated that about 17.6% of New Mexico experienced a decrease in productivity while 12.8% of the state experienced an increase in productivity. More than half of the state (55.6%) had insignificant change productivity, 10.8% was classified as indeterminant, and the remaining 3.2% was considered as agriculture was not evaluated.

The degradation in productivity was observed in about 2.2%, 4.5%, and 1.7% of NM's grassland, shrubland, and evergreen forest land cover classes, respectively. Simultaneously, about 5.7%, 1.3%, and 0.92% of NM's shrublands, grasslands, and evergreen forests were characterized by an increase in productivity, respectively. Regionally, significant decrease in productivity was observed in the northeaster and southeastern quadrants of the state while significant increase was observed in northwestern, southwestern, and a small portion in southeastern quadrants. The timing and number of detected breakpoints coincided with NM's drought frequency and severity and some fire incidents.

The TSS-RESTREN showed 55% agreement with the RPMS data over areas with significant decrease in productivity based on randomly selected pixels. Regionally, there was an agreement between the TSS-RESTREND and RPMS on the occurrence of significant degradation in productivity over the grasslands and shrublands within the Arizona/New Mexico Tablelands and in the Chihuahua Desert ecoregions, respectively.

A long timeseries assessment of rangeland productivity in New Mexico is critical to support decisions related to ecosystems management and conservation. The findings of this study can be used to address some of the rangeland degradation challenges that directly impact rangeland conditions, forage supply in the region, and New Mexican's livelihood as these systems support livestock production. In areas where degradation was prevalent (i.e., northeastern, and southeastern NM), special intervention would be needed to conserve the biodiversity of rangeland and increase the resilience of these ecosystems. As this study considered only precipitation, it was noticed that the rising temperature can play a significant role in vegetation. Thus, future analysis should consider its effects on rangeland productivity. Future research should also pay more attention to the association of degradation and productivity with recurring droughts.

**Author Contributions:** Conceptualization, M.G.G. and H.M.E.G.; data curation, M.G.G. and T.A.A.; formal analysis, M.G.G.; funding acquisition, H.M.E.G.; investigation, M.G.G., H.M.E. and T.A.A.; methodology, M.G.G., H.M.E. and T.A.A.; project administration, H.M.E.G.; resources, H.M.E.G.; software, M.G.G. supervision, H.M.E.G.; validation, M.G.G., H.M.E. and T.A.A.; visualization, M.G.G., H.M.E.G. and T.A.A.; writing—original draft, M.G.G.; writing—review & editing, H.M.E.G. and T.A.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially funded by the National Science Foundation (NSF), awards #1739835 and#IIA-1301346, and New Mexico State University.

**Data Availability Statement:** The data used in this study are openly available in data [43–48].

**Acknowledgments:** The authors would like to thank two anonymous reviewers for their time and effort in providing valuable comments and suggestions that helped in improving the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

#### **Appendix A**

Proportional allocation of the randomly sampled pixels as identified by Segmented VPR and Segmented RESTREND is shown in Table A1.

**Table A1.** Allocation of the randomly selected pixels as identified with Segmented VPR method over the different land cover classes based on the National Land Cover Dataset of 2011 (NLCD 2011) along the corresponding ecoregions.


31 = Bareland, 42 = Evergreen Forest, 52 = Shrub/Scrub, 71 = Grassland/Herbaceous, 82 = Cultivated Crops.

**Table A2.** Allocation of the randomly selected pixels as identified with Segmented RESTREND method over the different land cover classes based on the National Land Cover Dataset of 2011 (NLCD 2011) along with the corresponding ecoregions.



**Table A2.** *Cont*.

22 = Low Intensity Developed, 52 = Shrub/Scrub, 71 = Grassland/Herbaceous, 90 = Woody Wetlands, 95 = Emergent Herbaceous Wetlands.

#### **Appendix B**

The percentages of the randomly selected pixels (i.e., 155 identified by Segmented-VPR method) over grasslands and shrublands that were evaluated for their significance of the difference in mean productivity between before and after the break years at the pixel level using the Rangeland Production Monitoring Service (RPMS) data.

**Table A3.** The percentages of the randomly selected grassland pixels with respect to their significance in the difference in mean productivity before and after the break years based on the RPMS data along with their identified direction of change over the major ecoregions in New Mexico.


**Table A4.** The percentages of the randomly selected shrublands pixels with respect to their significance in the difference in mean productivity before and after the break years based on the RPMS data along with their identified direction of change over the major ecoregions in New Mexico.


#### **Appendix C**

A summary of the statistical significance test result of the difference in mean in productivity before and after the break years using the randomly selected pixels (i.e., 155 identified by Segmented-VPR method) based on the Rangeland Production Monitoring Service (RPMS) data at the ecoregion level.


**Table A5.** Significance of the difference in mean productivity before and after the break years based on the on increasing and decreasing shrublands and grassland randomly selected pixels using RPMS data at the ecoregion level.

\*\*\* *p* values less than or equal to 0.0001, \*\* *p* values less than 0.001, \* *p* values less than 0.05.

#### **References**


## *Article* **Persistent Vegetation Greening and Browning Trends Related to Natural and Human Activities in the Mount Elgon Ecosystem**

### **Dan Wanyama 1,2,\*, Nathan J. Moore <sup>1</sup> and Kyla M. Dahlin 1,3**


Received: 5 June 2020; Accepted: 30 June 2020; Published: 1 July 2020

**Abstract:** Many developing nations are facing severe food insecurity partly because of their dependence on rainfed agriculture. Climate variability threatens agriculture-based community livelihoods. With booming population growth, agricultural land expands, and natural resource extraction increases, leading to changes in land use and land cover characterized by persistent vegetation greening and browning. This can modify local climate variability due to changing land–atmosphere interactions. Yet, for landscapes with significant interannual variability, such as the Mount Elgon ecosystem in Kenya and Uganda, characterizing these changes is a difficult task and more robust methods have been recommended. The current study combined trend (Mann–Kendall and Sen's slope) and breakpoint (bfast) analysis methods to comprehensively examine recent vegetation greening and browning in Mount Elgon at multiple time scales. The study used both Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI) and Climate Hazards group Infrared Precipitation with Stations (CHIRPS) data and attempted to disentangle natureversus human-driven vegetation greening and browning. Inferences from a 2019 field study were valuable in explaining some of the observed patterns. The results indicate that Mount Elgon vegetation is highly variable with both greening and browning observable at all time scales. Mann–Kendall and Sen's slope revealed major changes (including deforestation and reforestation), while bfast detected most of the subtle vegetation changes (such as vegetation degradation), especially in the savanna and grasslands in the northeastern parts of Mount Elgon. Precipitation in the area had significantly changed (increased) in the post-2000 era than before, particularly in 2006–2010, thus influencing greening and browning during this period. The greenness–precipitation relationship was weak in other periods. The integration of Mann–Kendall and bfast proved useful in comprehensively characterizing vegetation greenness. Such a comprehensive description of Mount Elgon vegetation dynamics is an important first step to instigate policy changes for simultaneously conserving the environment and improving livelihoods that are dependent on it.

**Keywords:** bfast; Mann–Kendall; Sen's slope; East Africa; NDVI; breakpoint analysis; vegetation trends; greening; browning; Kenya; Uganda; trend analysis; land use; land cover

#### **1. Introduction**

Vegetation plays very important roles in ecosystem processes, including the mitigation of climate change effects [1] and the regulation of land surface temperatures, carbon and energy cycles [2–4]. At the same time, significant changes in terrestrial vegetation have been reported in the recent past [5,6]. Previous studies have concluded that these changes are driven by (1) slowly-changing natural processes, such as regional climate change (e.g., changes in temperature, precipitation, etc.), nitrogen disposition and increasing atmospheric CO2 concentrations [7], and (2) more rapid anthropogenic activities, including land-use and land-cover (LULC) change (e.g., deforestation [8,9], overcultivation and overgrazing [8], afforestation, expanding green areas in cities [1] among others). These processes do not operate in isolation [7] but rather interact at multiple scales, with global-scale drivers interacting with processes at the regional and local scales [10], thus making vegetation change dynamics a complex phenomenon to examine. Understanding vegetation dynamics has attracted substantial attention in the past years, specifically in the wake of climate change and variability [11]. As such, knowledge of vegetation response to climate change is critically important in the effort to maintain the supported processes as well as the human livelihoods derived from them [12]. Such an analysis has been difficult in the past, due in part to limited access to consistent data both in space and time. However, with recent developments in Earth observation (EO) technologies, spatio-temporally contiguous remote sensing (RS) data have been collected, making it possible to investigate vegetation dynamics accurately and comprehensively in terms of other environmental processes, at multiple scales.

The Normalized Difference Vegetation Index (NDVI), which is the normalized sum of the difference in reflectance between near infrared and red bands, has extensively been used [5,10,13] as an indicator of photosynthetic activity and vegetation amount [14]. Previous studies used this vegetation index (VI) to characterize vegetation processes such as productivity decline [6], phenology [15,16] and greening and browning [1,5,10,12,17,18]. Variability in NDVI generally follows patterns of climatic conditions, mainly precipitation [11,19–22]. The NDVI–precipitation relationship is largely a strong positive one, particularly in locations where total annual precipitation is less than 1000 mm [19,22]. This relationship is further compounded in space and time by site-specific factors such as existing LULC, vegetation structure and composition [22], topography [10] and soil type [23]. As a result, vegetation greenness response to climate has shown significant spatial and temporal variability. For instance, contrasting natures and strengths of the NDVI–precipitation relationship have been reported: a strong linear relationship in the Sahel [19] and a log-linear one in East Africa [22]. In other studies, seasonal precipitation with different lag times was found to correlate strongly with NDVI [13,22]. This complexity implies that (i) results from a specific regional analysis are not transferable to another region; (ii) any slight spatial and/or temporal misspecification may lead to misleading results about vegetation greenness–precipitation relationships and patterns; and (iii) understanding vegetation change requires multiple images (time series, TS) as opposed to the extensively used traditional classification and change detection methods.

East Africa depends heavily on rainfed agriculture, which puts food security and rural livelihoods at risk [24]. The importance of land in this region cannot be overemphasized [24], yet land holdings are small and steadily declining [24,25]. At the same time, populations have been increasing thus necessitating expansion of food production while also navigating the effects of climate change in the area. Over recent decades, therefore, East Africa has witnessed significant landscape transformation due to both human and natural drivers [6,26]. Generally, natural vegetation has been converted to farmlands, grazing lands and human settlements [25]. This LULC transformation has been reported in the Mount Elgon ecosystem (MEE), a major water catchment tower supplying water to three major lakes in East Africa (Lake Turkana, Kenya, Lake Kyoga, Uganda and Lake Victoria, Kenya, Uganda and Tanzania [27]). The MEE is dominated by croplands in most locations, mixed vegetation (primarily savanna, grasslands, and shrubs) in the northern portion and the Afromontane forest (Figure 1). The high population growth and densities in the area have translated into need for more land [28]. Coupled with political interference and corruption among park and reserve staff, the need for more land has resulted in forest encroachment and deforestation as ecologically fragile land is cleared for agriculture and settlement [28–32]. The changes in LULC have, in part, altered the functioning of the ecosystem [29] and, as a result, the mountain area has experienced more frequent landslides,

prolonged droughts, and flooding [30,31]. Evidence of a changing climate has been reported [33] and this may result in increased frequency of these events. As such, the livelihoods of at least 2 million people [27,31] are threatened. Despite this, the complex MEE landscape is currently understudied. A key study of LULC change was conducted by Petersson, Vedeld, and Sassen [27] and employed institutional theory in analyzing processes that led to deforestation within protected areas (PAs) in the transboundary MEE. It was found that, especially on the Kenyan side, it was challenging to correctly quantify LULC change due to the overlap between bamboo, plantation and Shamba system farms. Other studies have been conducted on relatively smaller spatial scales, mostly on the Ugandan side. Such studies have focused on the effect of LULC change on landslide occurrence [28], soil organic carbon, food security and climate change vulnerability [34], carbon stocks and climate variability [35] among others. However, some of the studies have reported contradictory results especially about the magnitude of LULC change within the agricultural land-use class. This may be due to the complex LULC orientation, which leads to persistent greening and browning in the MEE, thus making it difficult to correctly characterize vegetation dynamics especially using traditional classification and change detection methods. More robust methods are therefore needed, and TS analysis has been applied to comprehensively examine spatio-temporal landscape changes, particularly for constantly variable landscapes like the MEE.

**Figure 1.** Mount Elgon ecosystem (MEE) land-use and land-cover (LULC) in 2018. This map was created by reclassifying Moderate Resolution Imaging Spectroradiometer (MODIS) land-cover data (MCD12Q1) created by the University of Maryland [36]. Mixed vegetation includes savanna, grasslands, and shrubs. Cropland includes cropland and the cropland-vegetation mosaic classes. The black-gray line represents the boundary between Kenya and Uganda and the red triangle is Wagagai Peak, the highest point on Mount Elgon (4321 m above sea level).

TS analysis of RS data has been used to characterize environmental phenomena by describing both trends and discrete change events [37]. In recent years, application of TS methods has increased and this has been driven by improved access to RS imagery (for instance, due to opening up of the Landsat archive in 2008 [38]), improvements in the integration of RS and GIS (Geographic Information Systems) and general advancements in computing power [39]. As such, analysis of LULC change, including

vegetation greening and browning, has significantly evolved from the traditional bitemporal image analysis to using multiple and continuous observations. Common methods for TS analysis include Fourier analysis [40,41], principal components analysis, [42,43] and the Mann–Kendall statistic [44]. The Mann–Kendall statistic has been used to identify the presence and nature of monotonic trends in vegetation time series [5,10,45]. It is a non-parametric statistic; thus, data does not have to conform to any specific distribution [46]. Besides, Mann–Kendall compares relative magnitudes of sample data instead of raw data values and therefore missing values are allowed in the analysis. The Mann–Kendall analysis is often followed by the Sen's slope estimator [47], which quantifies the strength of the monotonic trends [12,45]. The Sen's slope estimator calculates the median of the set of slopes generated from Mann–Kendall [47]. Mann–Kendall and Sen's slope have been found to be more robust for TS with outlying values as compared to parametric statistics like ordinary least squares [12]. These two statistics identify and quantify any overall trends in a time series and are therefore well suited for examining vegetation greening and browning. Previous studies have successfully used these methods in assessing the consistency of greening and browning patterns across spatio-temporal scales in northern India [12] and assessing variability in greening and browning patterns caused by use of different RS imagery in the boreal forest of central Canada [48]. However, the assumption that the vegetation trend preserves its change rate throughout the period of study means that some greening and browning changes are masked [5]. For instance, later greening in a consistently browning vegetation may not be detected using Mann–Kendall and Sen's slope. To counter this issue, more TS decomposition methods have been proposed, including *Breaks for Additive Season and Trend*, bfast [16]. Using bfast, even subtle changes in vegetation can be monitored. This algorithm has successfully been used in delineating anthropogenic, fire and elephant damage within forest ecosystems in Kenya [49]. Complementing Mann–Kendall and Sen's slope with bfast can therefore be valuable in examining vegetation greening and browning trends, especially in a dynamic environment like the MEE. In such an analysis, the latter can be used to characterize changes detected by the former. Vegetation trends can then be characterized in more detail and this can be the basis for understanding effects of climate on terrestrial ecosystems [5] and the development of ultimate strategies for the sustainable management of ecosystems [10].

The goal of this study was to characterize comprehensively, over multiple time scales, recent patterns and trends in MEE vegetation greening and browning. Here, the main objectives were; (1) to assess and quantify the nature and magnitude of change in MEE greenness for the period 2001–2018; and (2) to characterize trends and variability in MEE precipitation as a way to disentangle nature- versus human-driven vegetation greening and browning. It is hypothesized that (1) changes in climate have forced local communities in the MEE to expand croplands at the expense of the natural vegetation thus leading to deforestation and degradation; and (2) the high variability exhibited in the MEE landscape requires integration of both general (such as Mann–Kendall and Sen's slope) and sequential (such as bfast) TS analysis methods to be fully characterized. To achieve these goals, the study analyzed spatio-temporal trends and patterns in Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI (2001–2018) and Climate Hazards group Infrared Precipitation with Stations (CHIRPS) precipitation (1986–2018) at multiple temporal scales (dekadal (10-day), 16-day, seasonal), using an integration of Mann–Kendall, Sen's slope and bfast algorithms. The study also incorporates inferences from a field study conducted in the MEE in 2019 to explain some of the vegetation change dynamics in the area. This analysis thus produces a more comprehensive characterization of vegetation dynamics in the MEE, which would not be possible using traditional classification and change detection methods. Results would help fill in some of the existing gaps in literature about nature and magnitude of LULC change and the stability of LULC in the MEE. Such a comprehensive description of MEE vegetation dynamics is an important first step to initiate dialogue aimed to instigate policy changes for simultaneously conserving the environment and improving livelihoods dependent on it.

#### **2. Study Area Description**

The MEE is located in western Kenya and eastern Uganda (Figures 1 and 2). The studied area covers approximately 15,000 km<sup>2</sup> and extends from 1◦37 42.82 N, 33◦55 45.07 E and 0◦42 15.76 N, 35◦14 18.84 E. Mount Elgon is a solitary volcano and is among the oldest in East Africa [28,32]. The highest point, Wagagai Peak, is 4321 m [28,50] and is found on the Ugandan side. This area rises from a plateau that lies 1850–2000 m in the east and 1050–1350 m to the west with a caldera that extends 8 km wide [50]. Vegetation in this area is zoned by altitude and mountain forest, farmland and Afro-Alpine heath and moorland are the common land covers [27]. Declared a protected area in 1968 and 1992 in Kenya and Uganda, respectively [31,51], Mount Elgon Forest, a montane rainforest [52], is home to many important indigenous tree species [27].

**Figure 2.** Map of the MEE in eastern Uganda and western Kenya showing long-term (1986–2018) mean annual total precipitation (CHIRPS [53]). Generally, the driest parts are the grasslands in northeastern MEE. The study area is wettest in the south and around the mountain region. Major protected areas and towns are shown for reference. FR is shorthand for forest reserve, NR is national reserve, NP is national park, WS is wetland system and WR is wildlife reserve.

The MEE receives rainfall in a bimodal pattern (two rainy seasons) and, according to Mugagga, Kakembo and Buyinza [28], most of the rainfall is received between April and October on the Ugandan side (with mean annual amounts ranging from 1500 to 2000 mm). The Kenyan side receives long rains between March and June and short rains between September and November—average annual rainfall ranging from 1400 to 1800 mm [54,55]. There is minimal temperature variation for the area—an average minimum of 15◦C and a maximum of 23◦C on the Ugandan side [28] and 14 and 24◦C on the Kenyan side [55]. However, temperatures and precipitation have a strong variation with changes in altitude [27].

#### **3. Materials and Methods**

The present study integrated Mann–Kendall, Sen's slope and bfast in the analysis of NDVI and precipitation trends to characterize recent greening and browning patterns and trends within the MEE.

#### *3.1. Data and Sources*

This study utilized the following datasets (Table 1) in the analysis of greening and browning trends in the MEE.


**Table 1.** Description of original datasets used in the study.

#### 3.1.1. MODIS NDVI and CHIRPS Precipitation

This study used MODIS NDVI data in the TS analysis of spatio-temporal changes in MEE vegetation greenness signal. The 16-day NDVI composite MOD13Q1.V6 [56] with a 250-meter spatial resolution was obtained through AppEEARS (https://lpdaacsvc.cr.usgs.gov/appeears/) [57]. Data for the period 2001–2018 were used. While Landsat data [58] have a better spatial resolution and historical coverage, and are therefore better suited for this study, they were limited by the many data gaps in the TS over the MEE due to persistent cloud cover and Landsat's long (16-day) revisit time. As such, MODIS NDVI, which has previously been used successfully in similar studies in East Africa [6,22,59], was used in this study. This study also used CHIRPS precipitation data [53]. These data have a spatial resolution of 5 km and provide global daily and pentad records from 1981 to present. For this study, these data for the period 1986–2018 were obtained and preprocessed within Google Earth Engine (GEE) [60]. It is worth noting here that there were no observation data to assess the accuracy of the CHIRPS dataset over the MEE. However, this dataset has been used extensively in similar studies [12,61,62]. Moreover, this dataset was recently validated [63] and the results showed that the CHIRPS data are reasonably accurate in estimating rainfall over east and South Africa.

NDVI preprocessing involved quality assessment using the associated VI quality files. Pixels with low quality, high aerosol content, cloud cover and possible shadows were excluded during this exercise. NDVI and precipitation composites for specific time scales were then generated. First computed were mean NDVI TS during the wet season for the period 2001–2018. Here, imagery in April, May and June was used. The TS of mean NDVI for each 16-day period during the season were also created. As a result, there were generally two TS for each month and were labeled h1 and h2 for TS created from NDVI composites recorded roughly in the first and second half of the month, respectively. For precipitation, the TS of total amount over the wet season were computed for different periods; 1986–2018 (33 years); 1986–1996 (first 11 years); 1997–2007 (median 11 years); 2008–2018 (last 11 years); 1986–2005 (first 20 years); 1999–2018 (last 20 years); and 2001–2018 (18 years, similar to the NDVI TS length). While the present study was aimed at characterizing changes in vegetation greenness from 2001 to 2018, analysis of longer precipitation TS was necessary to understand any longer-term precipitation patterns in the MEE that may influence the vegetation patterns. TS were also generated of total precipitation amounts for each dekad over the wet season. In this analysis, these dekads for the month of April were labelled April d1 (for dates 1–10), April d2 (for dates 11–20) and April d3 (for the rest of the month). The same nomenclature was used for dekad TS in May and June. These were used to assess nature and strength of relationship between vegetation greenness and precipitation variability over a shorter time scale, in which case results from dekad precipitation TS were compared to results from 16-day NDVI. The 16-day NDVI composites were used here due to unavailability of

data to compute dekad NDVI composites. Finally, monthly precipitation composites were generated for the period 1986–2018, for use in analyzing breakpoints in MEE precipitation.

#### 3.1.2. Field-Collected Data

Environmental data from local communities and government officers within the MEE were collected using semi-structured interviews and direct observation in July–September 2019 (IRB Number: STUDY00002404). Most of the interviewees were from significantly changing landscapes (areas showing significant changes in vegetation greenness). The participants in this study were interviewed about perceptions of climate change and land-use change, and historical patterns of agriculture and land-use change. Interviews were conducted and written responses recorded using Qualtrics software [64]. The fusion of such qualitative data with quantitative RS data is important because indigenous and historical accounts of LULC change in such a constantly variable landscape can be used to fill in gaps that may not be fully explained using RS and GIS alone.

#### *3.2. Methods*

This section describes the methods and analyses performed on NDVI and precipitation TS to identify areas and characterize patterns of vegetation greening and browning within the MEE. These two analyses were performed in R [65] and trends were assessed at the 95% significance level. After analysis of monotonic trends in greenness and precipitation, it was necessary to perform breakpoint analysis on these time series, to detect any significant breaks within the data. These two types of analyses are described in more detail below. Figure 3 highlights major methods used in this study.

**Figure 3.** Flowchart of analysis methods used in the study.

#### 3.2.1. TS Analysis: Mann–Kendall and Sen's Slope

To analyze initial spatiotemporal changes in vegetation greening and browning trends in the MEE, the current study borrowed from methodologies presented by Landmann and Dubovyk [6]. However, since NDVI TS often do not meet assumptions for parametric analysis [48], the current study used the nonparametric Mann–Kendall test rather than linear regression. To calculate the Mann–Kendall test statistic, data values are evaluated as an ordered TS [46] and each value is then compared to all subsequent values [46,66]. The Mann–Kendall S statistic is initially assumed to be 0 and a value of 1 is added to (subtracted from) the test statistic if the value of an observation is higher (lower) than that of the previous observation [46,66]. There is no change to the statistic if the values are equal. Equation (1) below shows the Mann–Kendall test equation. High positive and low negative values of S, respectively, indicate increasing and decreasing trends, but the strength of the trend is statistically quantified by computing probability associated with S and the size of the data sample [46].

$$\mathbf{S} = \sum\_{k=1}^{n-1} \sum\_{j=k+1}^{n} \text{sign}(\mathbf{x}\_j - \mathbf{x}\_k) \tag{1}$$

where *sign xj* − *xk* <sup>=</sup> 1 if *xj* <sup>−</sup> *xk* <sup>&</sup>gt; <sup>0</sup> = 0 if *xj* − *xk* = 0

= −1 if *xj* − *xk* < 0 **Source:** Khambhammettu [46]

The magnitude of the trends was quantified using Sen's slope estimator [47]. These algorithms were used in this study first to characterize and quantify any general patterns in MEE greenness during the growing season (April, May, and June). The analysis then investigated more subtle changes in vegetation greenness over shorter (16-day) periods during the growing season. From these, areas of increasing (decreasing) greenness were mapped as areas of vegetation greening (browning).

To disentangle changes due to human and natural forcings, monotonic trends were assessed in the precipitation TS and magnitude of the trends quantified. The analysis was performed on precipitation series at different temporal scales—based on the length of the time series (including 33-, 20-, 18 and 11-year periods) and TS resolution (including dekad totals, and general growing season totals). From these analyses, areas with significant positive (negative) precipitation trends were mapped as areas experiencing significant and consistent wetting (drying) over the analysis period. A cross examination of results from greenness and precipitation trend analysis was conducted to distinguish any humanfrom nature-driven vegetation greening and browning. Using available very-high-resolution imagery from Google Earth Pro [67], the nature of LULC conversion was ascertained. Besides, the statistics of LULC change (including greening, browning pixels) were calculated.

#### 3.2.2. Breakpoint Analysis: bfast

To further understand any temporal patterns in MEE greenness and precipitation, the bfast algorithm was applied using the 'bfastSpatial' package (http://www.loicdutrieux.net/bfastSpatial/). There exists a predictable annual cycle of greening an browning in vegetation, and these coincide with the occurrence of rainy and dry seasons [49]. bfast, created by Verbesselt et al. [16] and Verbesselt and Zeileis [68], creates a best-fit seasonal regression model with a trend component for the time series [69]. This approach follows three main steps (i.) fitting a harmonic model based on a historical ('stable') period, (ii.) testing observations that follow the historical period for any structural breaks from the fitted model, and (iii.) calculating the magnitude of change which is the median residual between observed and expected values in the monitoring period [49,70]. While this approach has been applied on forested landscapes, the current study applied it to the whole of the MEE, whose LULC comprises forested, savanna, grassland, and cropland LULC (Figure 1). This was prompted by the unique LULC orientation in the MEE that makes it difficult to detect LULC changes via satellite image

change detection analysis [27]. bfast breakpoint analysis can sequentially monitor vegetation change on a yearly basis, thus making it a suitable method for assessing changes in such LULC.

The models used in this part of the study were parameterized based primarily on the study by DeVries et al. [70]. First, the time series was divided into historical and monitoring periods. A minimum of two years for the historical period is recommended when using MODIS 16-day data [49,68]. In this study, therefore, the period 2001–2004 was used as the initial historical period and was assumed to be generally stable before the start of the monitoring period. First-, second- and third-order harmonic models were fit on the NDVI data, the results were assessed, and the first-order model was finally selected as the suitable instance to use. The single-order model has been used previously with the assumption that vegetation phenology generally follows a similar trend [70]:

$$y\_t = a + \gamma \sin(\frac{2\pi t}{f} + \delta) + \varepsilon\_t \tag{2}$$

where *yt* dependent variable


As in the study by DeVries et al. [70], the trend component was excluded from the time series to reduce chances of yielding false breakpoints.

The change magnitude was calculated as the median of the residuals in the monitoring period. The median is thought of as a conservative measure, unlike the sum, that minimizes the chances of getting inflated magnitudes and therefore false breakpoints [70]. However, increasing the number of observations before and after a change event, by including long monitoring periods, yielded very high magnitudes and unrealistically numerous breakpoints. As such, this study elected to use sequential non-overlapping one-year-a-time monitoring periods. Here, the TS was trimmed to include the historical period plus one-year monitoring period. For monitoring the period from January to December 2005, for example, the TS for January 2001 to December 2005 was used, and the monitoring period was set to start in January 2005. This sequential monitoring approach is illustrated in DeVries et al. [70]. The approach is advantageous as it enables the assessment of subtle changes in vegetation within the MEE, especially alternating degradation and restoration that would go undetected using other methods. The default values for *h*, the minimal segment size between potentially detected breaks in the trend model given as a fraction relative to the sample size [49], were used in this study.

Breakpoint analysis was also performed on the monthly total precipitation TS. The analysis was performed on both 1986–2018 and 2001–2018 time series, to understand any longer-term patterns as well as recent changes in the precipitation. Here, the third-order harmonic model with the trend term was the best fit. While there was no direct application on record of bfast for precipitation breakpoint analysis, it is noted that the algorithm can be used for this purpose [71,72].

#### *3.3. Validation of Results*

A hybrid validation of analysis results was performed in this study. The collection of reference data borrowed from the methodology used previously by Landmann and Dubovyk [6]. The seasonal greening and browning map from Mann–Kendall and Sen's slope analysis (Figure 4) was linked to very-high-resolution imagery in Google Earth Pro and both qualitative and quantitative accuracy assessment of the results was performed. Locations of greening, browning and no change were investigated by visually interpreting historical imagery in Google Earth Pro. The selection of these testing points was based on the minimum size of detectable change (at least 250 m, the size of a MODIS NDVI pixel) and the availability of sufficient imagery to interpret change. As such, only locations

with very-high-resolution imagery in one of the three initial and final years (2001–2003, 2016–2018, respectively) were used. By visually assessing historical imagery spanning these periods, vegetation change, or lack thereof, could be observed. A total of 153 visually interpreted points (51 browned, 50 greened and 52 no change) were used. Greenness change in pixels at these locations was first recorded. The seasonal greening and browning map from Mann–Kendall and Sen's slope analysis (Figure 4) was also reclassified into greened, browned and no change classes. The testing points were then compared to the reclassified raster and accuracy measures including overall, producer's, and user's accuracies were calculated. The final points were locations with pixels that indicated highly discernible change (like deforestation, reforestation etc.).

It was challenging to evaluate sequentially monitored vegetation change since it was not possible to gather testing points due to many temporal gaps in images in Google Earth Pro. In most instances, one can rarely find images for successive years, thus making it difficult to validate year-to-year change results. Since the bfast algorithm detects even subtle changes in vegetation greenness, changes could not be discerned with high confidence and therefore calculating accuracy assessment statistics would be misleading due to inconsistencies in testing data [73]. As such, these results were only visually and qualitatively assessed, and a general trend of change in the pixels was interpreted rather than year-to-year change. Evaluating such time series results has been found to be challenging previously [37,74]. In this study, inferences drawn from the field interviews were also qualitatively incorporated in this exercise to explain some of the trends detectable in both trend and breakpoint analyses.

#### **4. Results**

This study highlights and characterizes, using the Mann–Kendall, Sen's slope and bfast algorithms, recent vegetation greening and browning trends and patterns in the MEE at multiple time scales. The results highlight portions of the MEE that experienced persistent and significant changes in vegetation greenness, as indicated by changing NDVI. The results from similar analyses of precipitation TS are also presented and, together, attempt to disentangle nature- from human-driven changes observed over the MEE landscape.

#### *4.1. Trend Analysis Results*

#### 4.1.1. Persistent Vegetation Greening and Browning in the MEE

Greening (browning) was defined as any significant increase (decrease) in NDVI as shown by either Kendall τ (Mann–Kendall and Sen's slope) or the magnitude of change (bfast). The results indicate various greening and browning patterns during the growing season (Figure 4) and near-half month (Figure 5) periods. During the growing season, greenness significantly increased in approximately 27% (3400 km2) of the study area. Here, NDVI increased at rates up to 0.025 per year. These locations were concentrated mostly within croplands, grasslands, and savanna (Figure 1). Significant browning was also evident, with NDVI in more than 1400 km2 (about 11% of the total area) decreasing by up to 0.035 per year over the analysis period. These areas were mostly located in the southwestern part of the MEE in Uganda. This location includes the Namatala wetland, which has experienced intensive conversion to agriculture and settlement in the past years [75]. Browning was also evident in other locations around the Mount Elgon ecotone and elsewhere in the MEE. Based on our visual assessment, most of these corresponded to areas where deforestation occurred over the analysis period. No significant trends were found in the rest of the MEE (62%, 7800 km2).

Analysis of greening and browning trends in NDVI for every 16-day period in the months of April, May and June showed similar patterns; most browning occurred in the southwestern MEE and greening elsewhere. Most of the changed areas experienced greening and browning at rates of up to 0.03 and -0.04, respectively. The highest proportions of land where greening occurred was found in the month of May, especially the May h1 period, in which about 20% (2400 km2) of the MEE experienced significant greening (Figure 5). Greenness increased in more than 1600 km<sup>2</sup> (13%) of the study area during May h2. More greening was experienced in the June h2 and April h2 (9% of MEE greened during both periods). The greatest proportion of browned areas was observed in June h1, where about 600 km2 (5%) of the MEE experienced vegetation greenness decline (Figure 5E). Most of the land within the MEE did not experience any significant change during these 16-day periods.

**Figure 4.** Map of significant changed (greened and browned) locations during the growing season. Slope values (Kendall τ), indicative of magnitude of change per time step, are shown here. White pixels indicate no significant change.

**Figure 5.** *Cont.*

**Figure 5.** (**A**–**F**) Greening and browning trends in the Normalized Difference Vegetation Index (NDVI) for every 16-day period in the months of April, May and June for 18 years (2001–2018).

#### 4.1.2. Precipitation Variability in the MEE

The test for monotonic trends in precipitation revealed various patterns of consistently increasing (wetting) and decreasing (drying) precipitation within the MEE (Figure 6). Most of the areas experienced increased precipitation over the years and only negligible proportions of the study area became drier. The greatest proportions of land during which consistently wetter conditions prevailed include 1999–2018, 1986–2018 and 2001–2018 (82%, 58% and 38%, respectively). The areas experiencing change covered approximately 10,300, 7200 and 4800 km2, respectively. These periods also showed the greatest magnitude of change in precipitation amount. Here, precipitation increased by at least 13, 4, and 5 mm per year for 1999–2018, 1986–2018 and 2001–2018 periods. In addition, about 27% of the MEE also experienced wetter conditions during 2008–2018.

There were very small portions of the MEE with wetter conditions during the earlier years in the analysis; 1.61% of land (200 km2) recorded wetter conditions in 1986–1996 while there was no significant increase in precipitation during the 1986–2005 and 1997–2007 periods. Consistently drier conditions were observed in two time periods (1986–1996 and 1986–2005). However, only negligible proportions of land (up to about 2%) experienced this change, although at substantial magnitudes of up to 9 mm per year. Different spatial patterns existed from 2000; significantly increasing precipitation in 2008–2018 was observed mostly on the Ugandan side, and the mountain forest and wetland reserves seemed to be excluded (Figure 6B). In 2001–2018, the changed pixels were mostly found within the mountain area, in the north and some areas in the west (Figure 6E). On the other hand, areas in the northeastern portion of the study area did not experience changes in precipitation during the 1986–2018 period (Figure 6F). Lastly, the 1999–2018 period had significant increases in total seasonal precipitation with an exception of a few eastern and southeastern portions of the MEE (Figure 6D). Overall, there was an increase in precipitation for most areas in the study area, with an increasing magnitude of change, especially in the post-2000 era.

**Figure 6.** (**A**–**F**) Maps of precipitation change for each analyzed time period.

The trend analysis for each dekad in the growing season found that, based on rates of change, the 2008–2018 period recorded the highest change in precipitation in April d1 and May d1, at rates of up to 13 mm per yearly dekad, see Figure 7). Areas where this change was experienced were located mostly in western MEE. An increase in precipitation was also recorded in April d2 in 2001–2018, June d3 in 1986–1996 and June d3 in 2008–2018. Negative trends in MEE precipitation were also observed in some dekads, mostly in 1986–2005. During this period, precipitation decreased in April d1, May d2 and June d1, at the rate of up to 3 mm per annual dekad.

Based on the proportions of land with significant change in precipitation, the study found that most of MEE precipitation significantly changed during May d1 in 1986–2018 (90%, 11,300 km2) and May d1 in 1999–2018 (43%, 5400 km2, Figure 8). From these results, some patterns are clear. For the months of May and June, annual dekad precipitation increased significantly for all time periods, although at varying rates and proportions. Precipitation in most of the MEE was more stable in the earlier years in the study (as in 1986–1996) or generally depicted significant and consistently drier conditions (as in 1986–2005). Based on these results and those in Figure 6 above, it was necessary to conduct a breakpoint analysis to provide a better characterization of any specific spatio-temporal breaks in precipitation and greenness trends.

**Figure 7.** (**A**–**C**) Significantly changed precipitation for dekads in the period 2008 to 2018.

**Figure 8.** (**A**–**E**) Significantly changed precipitation for dekads in the period 1999 to 2018.

#### *4.2. Breakpoint Analysis Results: bfast*

While breakpoint analysis was performed for both 1986–2018 and 2001–2018 precipitation TS, no significant breaks were found in the former. Therefore, results from 2001–2018 are presented in this study.

#### 4.2.1. MEE Precipitation

The breakpoint analysis in bfast revealed no significant breaks in MEE average mean total monthly precipitation. However, breakpoints monitored for each pixel from January 2005 revealed very interesting patterns. In most configurations, the analysis revealed significant breakpoints in 2006 and 2007 for most of the MEE. To reduce the influence of post-change detection observations on change magnitude, the same analysis was performed for the period 2001–2008, with 2001–2004 set as the historical period. Therefore, change maps and statistics provided are based on this adjusted analysis. Precipitation changed significantly in the two years, with about 9800 km2 (66% of the area) and 4800 km<sup>2</sup> (32%) of the MEE experiencing wetter conditions in 2006 and 2007, respectively (Figure 9). Most of the breakpoints in 2006 were detected in the months of July–October, while a great proportion of changes in 2007 were detected in March and April. Similar wetter conditions were detected in some locations in the months of May and June for both years. No drier conditions were detected. The magnitude of precipitation change ranged from approximately 10 to 53 mm during the monitored period (Figure 10A) and the bfast models used could explain up to 70% of the variance in the precipitation breakpoints (Figure 10B). Precipitation was also monitored sequentially for the period 2005–2018 and no breakpoints were detected except for a few pixels in 2005–2007.

**Figure 9.** (**A**,**B**) Months when major breaks were detected in the precipitation time series (2001–2008) when the period 2005–2008 was monitored.

**Figure 10.** (**A**,**B**) Adjusted R2 (greater than 0.5) and magnitude of change in the precipitation time series (2001–2008) using 2005–2008 as the monitoring period.

#### 4.2.2. MEE Greenness

The analysis of breakpoints revealed that some significant breaks existed in the NDVI time series for each year (2005–2018). Having been sequentially monitored, the results include, for each monitored year, months when the breakpoints were observed, the magnitude of change at the breakpoints, the length of historical data used and adjusted R2. In this analysis, only highly statistically significant changes (*p* < 0.05) are reported. This decision was based on two reasons: (i.) There were no ground data for validating results from sequentially monitored bfast. To ensure that only accurate results are reported, breakpoints from models with less than 50% adjusted R<sup>2</sup> were excluded. Thus, only breakpoints from models with average to high explanatory power were reported. This 50% threshold has been used elsewhere by Landmann and Dubovyk [6]. (ii.) Vegetation greenness in the MEE showed significant interannual variability. As such, to ensure a long-enough historical period used by 'bfastSpatial', the study excluded any breakpoints from models using a historical period of less than two years.

The 'bfastSpatial' model was able to detect changes in vegetation greenness, especially in the grasslands of the northeastern MEE (Figures 1 and 11). The observed changes could be due to both natural factors (such as variability in precipitation) and/or human activities (for example, clearing of land for agriculture and settlement, deforestation for charcoal burning and construction etc.). The maximum magnitude of changes ranged from −0.24 and 0.21, thus indicating only subtle changes in the MEE vegetation. The breakpoints were detected in each of the years, but most of them were observed in 2013, 2007 and 2010 (greening) and 2009 and 2017 (browning), as shown in Figure 11. Overall, there were more areas where significant greening trends were detected compared to browning, with about 17% of the MEE (about 2500 km2) and 10% (about 1500 km2) showing greening and browning, respectively, over the monitoring period. These represent annual greening and browning rates of about 1.2% and 0.7%, respectively.

**Figure 11.** *Cont.*

**Figure 11.** Greening and browning in the MEE. The results from bfastSpatial using NDVI time series (2001–2018) while sequentially monitoring each year in 2005–2018. These maps show magnitudes of change for pixels with significant breakpoints.

Based on the yearly changes, differentiated greening and browning patterns were observed during the monitoring period. Cumulatively, the greatest proportion of land with changed greenness was in 2013 in which significant breakpoints were observed in approximately 13% of the MEE (approximately 1850 km2). This was followed by 2009 (1350 km2), 2007 (1200 km2) and 2010 (650 km2) (Figure 12). For 2013, 2010 and 2007, the majority (over 95%) of the changed locations experienced greening. The year 2009 recorded browning in over 99% of the changed locations.

**Figure 12.** Greening and browning in the MEE. The results from bfastSpatial using NDVI time series (2001–2018) while sequentially monitoring each year in 2005–2018. The graphs show the total area of land that experienced greening and browning for each monitored year.

#### 4.2.3. MEE Greenness vs Precipitation

There was an increase in precipitation in 2006 compared to previous years (Figure 13). There followed a steady decrease in 2007–2009 and finally an increase in 2010. Similar changes were observed in greenness during the five years. First, there was greening at most breakpoints in 2007 following the significant wetting in 2006. Significant browning followed, most of which occurred in 2009. This was also observable in an aspatial bfast breakpoint analysis performed on MEE mean NDVI, in which a sudden increase in NDVI was observed towards the end of 2006, followed by a consistently reducing trend until 2009/2010, in which another sudden increase was found (Figure 14). No significant breakpoints were found from a similar analysis using mean MEE precipitation. However, results from 'bfastSpatial' indicate that most of the MEE had very significant increases in precipitation in the second half of 2006 and the first half of 2007. Therefore, the corresponding changes in vegetation greenness can be attributed to this change in precipitation.

**Figure 13.** Violin plot of average total annual precipitation in the MEE. The red box highlights the period 2006–2010.

The significant and extensive browning observable mostly in 2009 (Figure 11) may be explained in terms of vegetation regaining its 'normal' greenness levels following sudden greening due to an above-normal precipitation. However, no significant breaks were found in precipitation around 2013, the year with the highest proportion of land with detected change in vegetation greenness. Thus, such browning and greening trends can be explained with regards to other factors, including temperature and anthropogenic activities that may have altered vegetation greenness during the monitoring period. The 'bfastSpatial' model did not find any significant breakpoints, in precipitation TS, in 2012. This may be due to the unstable historical data prior to the monitoring period. However, the violin plot (Figure 13) indicates that this year had some of the highest precipitation recordings. As such, there is a high chance that the greening following in 2013 was influenced by this increase in precipitation.

A visual inspection of these results in Google Earth Pro indicates that the subtle changes are indicative of various LULC changes. There was evidence of human settlement being introduced into the grasslands in the northeastern MEE. Information gathered from the field corroborates this finding and further explains the implications for landscape greenness. Inhabited by nomadic pastoralist communities, this part of the MEE is susceptible to degradation, especially when these communities move, driven by rainfall patterns, to settle within the grasslands. Based on fieldwork and data from Google Earth Pro, reduced natural vegetation and tree density were evident in these locations. In other instances, the new inhabitants converted natural vegetation to agriculture and, although this would result in environmental degradation, some planted crops were greener than the natural vegetation, and, therefore, these areas were shown to exhibit greening trends. In other locations outside of the grasslands, visible greening trends were indicative of some afforestation practices. Data from the field revealed that this kind of greening was driven by the cultivation of evergreen early maturing blue

gum (*Eucalyptus globulus*) tree species, sometimes together with and other times in the place of the maize crop. This was especially common in the eastern and southeastern parts of the domain and parts of Uganda.

**Figure 14.** bfast mean 16-day NDVI time series decomposition (2001–2018). Breakpoints in trend were found between 2006 and 2010.

#### *4.3. Accuracy Assessment*

Validation of Man-Kendall and Sen's slope results revealed an overall accuracy of 98.04% and user's accuracies of 100%, 98% and 96.2% for browned, greened, and unchanged locations, respectively (Table 2). Producer's accuracies of 98% (for both browned and greened locations) and 98.1% (no change) were also obtained. The visual inspection of bfast results using Google Earth Pro also revealed that most of the detected subtle changes occurred. However, no more detailed information like year of change could be discerned due to the lack of available high-resolution imagery.


**Table 2.** Mann–Kendall and Sen's slope accuracy assessment statistics.

#### **5. Discussion**

#### *5.1. Precipitation and Vegetation Change in the MEE*

Climate change continues to affect economies in most developing nations, especially those relying heavily on natural processes for their livelihoods. Agriculture remains the backbone of many of these nations [76–82], yet agriculture's vulnerability to climate change effects cannot be contested. Climate-related natural hazards, including extensive flooding, extended droughts [83], landslides [28] among others, have significantly impacted agricultural production and endangered lives. The frequency of these events shows an increasing trend [84,85], thus trapping many agriculture-dependent communities in an unending struggle for survival. Therefore, existing food insecurity in these developing nations can be attributed to their overreliance on rainfed agriculture [83]. Moreover, the high population growth and densities in these nations and elsewhere have translated into need for

more agricultural land [28]. Coupled with political interference and corruption among forest park and reserve staff, this need for more land has resulted in forest encroachment and deforestation as ecologically fragile land is cleared for agriculture and settlement [28–32]. Therefore, understanding the major drivers of landscape change is an important first step to inform better decisions to simultaneously conserve the natural environment and improve the livelihoods dependent on it. First, being able to separate nature- from human-induced landscape changes would be valuable in this endeavor.

Changes in the landscape occur at varying rates and magnitudes across space and time, from very subtle tree damages to forest clearings [49]. In this study, two major forms of landscape change browning (areas of declining vegetation) and greening (areas of increasing vegetation) [12]—were studied. The results show that MEE greenness exhibited substantial variability, and some form of the greening and browning change was recorded each year. As expected, these changes varied by scale, with the highest proportions of greened and browned locations observed over the growing season rather than any of the individual 16-day periods. Importantly, there was a lot of activity in areas bordering the mountain forest, as expected. Clearly, both greening and browning trends were observable, particularly on the Kenya side. Significant deforestation occurred as a result of encroaching fertile land on the high slopes of the mountain, for agriculture and settlement (examples in Figure 15). This finding has been reported in previous studies [27,28]. The Shamba system, thought of as a win-win arrangement, enabled local communities to farm in protected areas while tending to the growing trees in their early stages of growth [27]. This was a government effort originally meant to convert native to plantation forests and later to replant trees on harvested forest land. However, the Shamba system farms overlap with plantation forest and bamboo vegetation thus making it more difficult to conclusively identify, characterize and quantify the nature and magnitude of LULC change, especially by the use of traditional classification and change detection methods [27]. The significant and consistent browning in the southwestern MEE are attributed to the conversion of the Namatala wetland to agriculture (mostly paddy rice farming) and settlement (due to the growth of Mbale town) [75]. Since the early 2000s, 80–90% of this wetland has been converted. The wetland area is an Important Bird Area (IBA) [86] and therefore the reported LULC conversion caused many nature–human and human–human conflicts, including, respectively, bird poisoning and competition among people to own and utilize the wetland.

Disentangling nature- from human-induced vegetation change is an important yet complex task. In this study, patterns of change in precipitation varied with the TS duration (33-, 11-, 20- and 18-years) and resolution (dekad, seasonal). Previous studies indicated that precipitation in the area exhibited significant temporal variability, with both positive [63] and negative [80] trends observed over time. It can be inferred that precipitation had changed more (increased) over the period after 2000 than earlier. This implies that some of the vegetation greening and browning should be linked to this MEE wetting. This is true for the period 2006–2010, where bfast reveals that greening and browning events in the MEE follow significant wetting and drying events. However, a visual interpretation of Mann–Kendall analysis maps for both MODIS NDVI and CHIRPS precipitation did not show much similarity. Besides, bfast did not reveal any breakpoints in precipitation around 2013, the year when most breakpoints in greenness were detected. While this may be a fault in the bfast model used (because of, for instance, instability in historical data used), precipitation alone may not be a reliable predictor of vegetation change [49]. Davenport and Nichols [22] concluded that the NDVI–precipitation relationship varies both spatially and temporally, is not linear, generally exhibits a three-month lag, and depends on underlying LULC types. Therefore, it is likely that these complex relationships were not well revealed in both analyses in this study. Besides, the observed greening and browning patterns may be driven by other climate factors (such as temperature), which were not examined in this study, due in part to data unavailability for the MEE. Moreover, the spatial scale mismatch between the CHIRPS and MODIS data make more quantitative comparisons challenging.

**Figure 15.** (**A**–**D**) Some of the greened and browned areas identified by Mann–Kendall and Sen's slope. Images (**A**–**C**) indicate the conversion of natural vegetation (mostly forest) to cropland and settlement. Image (**D**) indicates afforestation. Source: Google Earth Pro [67].

Visual inspections of the results from both analyses revealed the applicability of Mann–Kendall and bfast methods in detecting changes in vegetation greenness. Most people in the MEE are small-scale farmers, owning sometimes less than an acre (0.4 hectares) of land. This represents less than 10 percent of the pixel size used (250 by 250-meter MODIS), meaning that some subtle land conversions were missed. The choice to use these data was necessitated by the frequent availability of MODIS NDVI data. The ability of Mann–Kendall and bfast in delineating areas of change in the MEE using this coarse imagery is important, as free finer RS imagery like Landsat is limited by both cloud cover across tropical regions and their coarse temporal resolutions. Looking forward, combining Sentinel-2 and Landsat imagery may help address this issue, but this will not resolve the problem for historical studies like this one.

The use of Mann–Kendall and bfast algorithms proved to be a valuable integration. Mann–Kendall performed well in mapping locations that had experienced significant change in the entire TS (2001–2018). Consistently greened and browned locations were mapped with high confidence (95% significance level) and similar patterns existed for both 16-day and average growing season periods: significant and consistent browning in southwestern MEE and greening elsewhere. The best and most accurate results were obtained using average growing season NDVI TS. Here, browned locations (especially around mountain boundaries and southwestern MEE) and greening elsewhere were clearly demarcated. This indicates that the growing season period, rather than the shorter 16-day period, is the best temporal scale for monitoring vegetation change in the MEE. bfast on the other hand performed well in mapping changed greenness locations for sequentially monitored periods from 2005 to 2018. This part of the analysis revealed only subtle changes in vegetation greenness in most locations, indicative of vegetation degradation rather than deforestation [70]. Importantly, areas with breakpoints across multiple years

were effectively identified, a finding that would likely be omitted if using traditional post-classification change detection methods. These findings demonstrate that the process of vegetation greening and browning can be studied more thoroughly by fusing these two methods. Using Mann–Kendall and Sen's Slope, one can assess and quantify monotonic trends in their TS to get the general picture in the series. This can be complimented with an assessment of the temporal occurrence of significant breaks in the series using bfast. Using this integrated approach, vegetation greening and browning can be fully characterized and understood to provide more information for better decision making.

#### *5.2. Sources of Uncertainty*

Due to a lack of ground precipitation measurements, the present study did not validate the CHIRPS data over the MEE. The GHCN data [87] would potentially be used, but only one data point existed in the study area and was therefore deemed insufficient for this task. While this is a potential source of uncertainty, these data have been used extensively in similar studies. Moreover, CHIRPS data are based on ground station data since they are created through "smart interpolation" procedures that incorporate both satellite information and gauge data [53] thus making daily CHIRPS imprecise, and reliable only at dekad or higher aggregations. Validation of bfast results was also not possible due to a similar lack of data. Frequent, high-resolution imagery was not freely available to download and process, and the study relied on imagery from Google Earth Pro which had a lot of spatial and temporal gaps. While no accuracy statistics were calculated for these results, the results were interpreted based in part on data collected from the MEE field study.

#### **6. Conclusions**

The MEE of eastern Uganda and western Kenya was found to exhibit significant variability in vegetation dynamics and precipitation regimes. This variability was attributed to the existing LULC orientation especially in eastern MEE and climate change and variability. As such, it is highly probable that analysis of only a few images to ascertain MEE landscape change would yield inconsistent results. In this study, greening and browning in the MEE was examined using both TS trend and breakpoint analysis methods. The MEE had experienced significant and persistent greening and browning at different time scales and this change was attributed to both natural factors (including changing precipitation) and anthropogenic factors (especially the vegetation-to-cropland conversion). The southwestern MEE had consistently browned due to the conversion of the Namatala swamp to paddy rice farming and settlement. A lot of activity was also observed around the mountain forest boundary as people encroached and converted the forest LULC to agriculture and settlement. There were breakpoints in the vegetation greenness TS, particularly in the savanna and grassland land covers in northeastern MEE. The breakpoints were detected in each of the monitored years (2005–2018), but most of them were observed in 2013, 2007 and 2010 (greening) and 2009 and 2017 (browning). The study also concluded that MEE precipitation had significantly changed (increased) in the post-2000 era. More specifically, total precipitation significantly increased in 2006 and 2009–2010 with a consistently decreasing trend in between. We therefore concluded that these precipitation changes influenced significant greening and browning patterns observed in the same period. The greenness–precipitation relationship was weak in other periods as greening and browning changes were not strongly influenced by changing precipitation. This may be attributed to the complex nature of the MEE landscape and/or the spatial and temporal scale mismatch between MODIS NDVI and CHIRPS precipitation data. The integration of Mann–Kendall, Sen's slope and bfast proved useful in comprehensively characterizing recent changes in vegetation greenness within the MEE. Having a comprehensive description of vegetation change is an important first step, especially for such a variable landscape, to effect policy changes aimed at simultaneously conserving the environment and improving livelihoods that are dependent on it.

**Author Contributions:** Conceptualization, D.W. and N.J.M.; methodology, D.W., N.J.M. and K.M.D.; software, D.W.; validation, D.W. and K.M.D.; formal analysis, D.W.; investigation, D.W. and N.J.M.; resources, D.W. and N.J.M.; data curation, D.W.; writing—original draft preparation, D.W.; writing—review and editing, D.W., N.J.M. and K.M.D.; visualization, D.W. and N.J.M.; funding acquisition, D.W. and N.J.M.; supervision, N.J.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by funds from the 2019 Graduate Office Fellowship (GOF) awarded by College of Social Science and Department of Geography, Environment, and Spatial Sciences at Michigan State University.

**Acknowledgments:** Multiple people made this publication a success and authors would like to express their gratitude especially to Eliud Akanga for your valuable support in the fieldwork portion of the study.

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

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Earth Observation-Based Detectability of the Effects of Land Management Programmes to Counter Land Degradation: A Case Study from the Highlands of the Ethiopian Plateau**

**Esther Barvels \* and Rasmus Fensholt**

Department of Geosciences and Natural Resource Management, University of Copenhagen, DK-1350 Copenhagen, Denmark; rf@ign.ku.dk

**\*** Correspondence: lbm692@alumni.ku.dk

**Abstract:** In Ethiopia land degradation through soil erosion is of major concern. Land degradation mainly results from heavy rainfall events and droughts and is associated with a loss of vegetation and a reduction in soil fertility. To counteract land degradation in Ethiopia, initiatives such as the Sustainable Land Management Programme (SLMP) have been implemented. As vegetation condition is a key indicator of land degradation, this study used satellite remote sensing spatiotemporal trend analysis to examine patterns of vegetation between 2002 and 2018 in degraded land areas and studied the associated climate-related and human-induced factors, potentially through interventions of the SLMP. Due to the heterogeneity of the landscapes of the highlands of the Ethiopian Plateau and the small spatial scale at which human-induced changes take place, this study explored the value of using 30 m resolution Landsat data as the basis for time series analysis. The analysis combined Landsat derived Normalised Difference Vegetation Index (NDVI) data with Climate Hazards group Infrared Precipitation with Stations (CHIRPS) derived rainfall estimates and used Theil-Sen regression, Mann-Kendall trend test and LandTrendr to detect changes in NDVI, rainfall and rain-use efficiency. Ordinary Least Squares (OLS) regression analysis was used to relate changes in vegetation directly to SLMP infrastructure. The key findings of the study are a general trend shift from browning between 2002 and 2010 to greening between 2011 and 2018 along with an overall greening trend between 2002 and 2018. Significant improvements in vegetation condition due to human interventions were found only at a small scale, mainly on degraded hillside locations, along streams or in areas affected by gully erosion. Visual inspections (based on Google Earth) and OLS regression results provide evidence that these can partly be attributed to SLMP interventions. Even from the use of detailed Landsat time series analysis, this study underlines the challenge and limitations to remotely sensed detection of changes in vegetation condition caused by land management interventions aiming at countering land degradation.

**Keywords:** developing countries; Google Earth Engine; land degradation; Landsat time series analysis; semi-arid areas; sustainable land management programmes

#### **1. Introduction**

Degradation of land and soil affects approximately one third of the global land area that is used for agriculture [1], involving livelihoods of more than 1.5 million people [2]. Land degradation is of particular concern in developing countries, as this issue poses a threat to food security for a large number of poor people and to local economic activities [3]. In the United Nation's Convention to Combat Desertification (UNCCD) (art. 1f) land degradation is defined as "reduction or loss, in arid, semi-arid and dry sub-humid areas, of the biological or economic productivity and complexity of rainfed cropland, irrigated cropland, or range, pasture, forest and woodlands resulting from land uses or from a process or combination of processes, including processes arising from human activities and habitation patterns". Among those processes are soil erosion or long-term loss of

**Citation:** Barvels, E.; Fensholt, R. Earth Observation-Based Detectability of the Effects of Land Management Programmes to Counter Land Degradation: A Case Study from the Highlands of the Ethiopian Plateau. *Remote Sens.* **2021**, *13*, 1297. https://doi.org/10.3390/rs13071297

Academic Editor: Elias Symeonakis

Received: 30 January 2021 Accepted: 9 March 2021 Published: 29 March 2021

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

**Copyright:** © 2021 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/).

natural vegetation [4]. In Ethiopia, land degradation results mainly from soil erosion by water [5] and occurs particularly in the Ethiopian highlands which are inhabited by 88% of the national population, cover 60% of the national livestock resources and encompass 90% of the area suitable for agriculture [6]. The country's natural physical conditions are one of the underlying drivers of land degradation. Ethiopia has always been prone to soil erosion and droughts due to high rainfall variability, which causes reduced vegetation cover in dry years and soil loss in subsequent wet years [7]. In Ethiopia's highlands, soil erosion is facilitated by steep terrain with slopes in excess of 30% [5]. Moreover, population pressure, increasing livestock (and with it deforestation), overgrazing (due to uncontrolled free grazing) and the expansion of agricultural fields into marginal land are underlying drivers of soil erosion [5,8]. The shortage of fertile cropland led to a shifting of cattle and livestock grazing activities to areas that are specifically vulnerable to soil erosion such as deforested, ecologically fragile hillsides with steep slopes. Gully formation, the removal of soil along drainage lines (channels) by surface water runoff, is one of the apparent consequences of soil erosion in Ethiopia [8].

The restoration of degraded land and soil, the implementation of sustainable land management (SLM) and resilient agricultural practices are targeted in the United Nation's Sustainable Development Goals (SDG). To address in particular target 15.3 which aims to combat desertification and restore degraded and soil, UNCCD adopted the Land Degradation Neutrality (LDN) Target Setting Programme. Achieving LDN will also contribute to reaching other SDGs including those on poverty reduction, food and water [9]. To counter desertification and land degradation UNCCD and affected developing countries have set voluntary targets and implemented National Action Programmes that are supported by international cooperation, including financial and technical resources [4]. In this context, developing countries have implemented land management and land restoration projects on both national and local levels in collaboration with multilateral and bilateral development partners. To rehabilitate degraded landscapes and scale up SLM in Ethiopia, the Ethiopian government launched the Sustainable Land Management Programme (SLMP) in 2009 in collaboration with a range of international donors including the World Bank and the German Development Bank (KfW) [10]. The impact of Sustainable Land Management (SLM) in Ethiopia has been studied through runoff and soil loss measurements [11] and analyses within the economics field, for instance by examining household data to assess the effect on crop yields [12]. A study by Ali et al. used earth observation derived vegetation indices to estimate the impact of SLM in a single watershed in Ethiopia [13].

Monitoring of land surface dynamics, such as land degradation ('browning') or land recovery ('greening') is widely done by implementing earth observation time series analysis. A plethora of studies have examined long term trends by applying ordinary least square (OLS) linear regression models, e.g., regressing vegetation indices with time, based on high temporal resolution data such as from Advanced Very High Resolution Radiometer onboard National Oceanic and Atmospheric Administration (AVHRR-NOAA) or from Moderate Resolution Imaging Spectroradiometer (MODIS). For the Sahel, linear trends of yearly NDVI anomalies derived from AVHRR [14] or linear trends of the seasonal NDVI amplitude and integral [15] have been examined to monitor vegetation. MODIS NDVI data have been used to detect land degradation and regeneration processes in the Sahel [16], Mongolia [17] and Ethiopia [18].

Studies of semi-arid areas [19–22] and of Ethiopia [23] have demonstrated a strong relationship between rainfall and NDVI. Therefore, methods have been developed to disentangle rainfall-related effects from human-induced effects on land degradation. The rain-use efficiency (RUE), defined as the ratio of above-ground net primary production (ANPP) to annual precipitation [24], has been used to detect non-precipitation related land degradation. The basic assumption involved in the use of RUE is the existence of a constant linear relationship between vegetation productivity (or ANPP) and precipitation in areas where land is not affected by human-induced degradation [25]. By normalising for the effect of interannual rainfall variability on ANPP, human-induced changes can then be singled

out [26]. RUE has been used as a measure in several studies, e.g., of the Sahel [25,27], South Africa [19], global drylands [2] and Northern Eurasia [26]. Furthermore, residual trend analysis, i.e., analysing the residuals from a NDVI-rainfall regression model, was found to be effective in disentangling the climate effects from human-induced land degradation [28,29]. In areas with steep terrain land degradation can be driven by climate effects such as high rainfall variability. Hermans-Neumann et al. combined NPP trends, precipitation variability and census data to identify areas in Ethiopia where high in-migration is coupled with land degradation, proposing the latter is likely occurred due to human activities [18].

Since the opening of the Landsat archive by the United States Geological Survey (USGS) in 2008, an increasing amount of studies that exploit medium/high resolution data for time series analysis has been published [30,31]. In parallel, new change detection methods that not only account for linear (in this case gradual trends), but also for abrupt occurrences by separating time series into individual segments, have been developed. Amongst those are Landsat-based detection of Trends in Disturbance and Recovery (LandTrendr) and Breaks For Additive Seasonal and Trend (BFAST) which have been widely used for vegetation monitoring. BFAST has for example been used to detect gradual and abrupt changes in NDVI and rain-use efficiency [26,32–34] and water-use efficiency [35]. LandTrendr has been widely used for monitoring forest disturbances (fire or stand clearing) and forest regrowth [36–38], forest biomass [39] and for agricultural and land abandonment mapping [40].

The aim of this study was to temporally and spatially analyse vegetation dynamics in degraded land areas in Ethiopia between 2002 and 2018 in relation to management programmes implemented to counter land degradation. Due the heterogeneity of the landscapes of the areas examined in this study and the fact that human-induced changes are expected to take place at a small spatial scale, this study aimed at exploring the value of using medium resolution Landsat data derived from different sensors for trend analysis at a spatial and temporal scale compatible with the scale of SLMP interventions. The investigated areas had gone through interventions aiming at avoiding further land degradation and increasing vegetation cover. In this context, associated human-induced and climate-related factors of land degradation and land recovery were examined. The specific objectives of this study were:


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

#### *2.1. Study Area*

The study area consists of 21 major watersheds which are distributed in three different zones of Ethiopia, in Amhara, Oromia and Tigray, and have mean altitudes between 1200 and 3100 m (Figure 1), mean slopes up to 14.4 degrees and mean annual rainfall between approx. 600 and 1900 mm per watershed (calculations based on CHIRPS data). The watersheds are located in semi-arid and sub-humid agro-ecological zones where temperate to cool climate prevails and are surrounded by low-lying tropical warm to hot savannas and semidesert regions [18,41]. They are characterised by heterogeneous landscapes, with croplands and grasslands as the most dominant land cover and hillsides, which have been degraded and closed for farming and grazing. Agriculture is dominated by small-scale subsistence mixed farming systems, i.e., crop production mixed with livestock rearing activities [41]. Crops are mainly grown during the wet period from March through September with harvest taking place mostly in October to December [18].

The watersheds constitute intervention areas of the Sustainable Land Management Programme, a multi-donor supported project, first implemented by the Government of Ethiopia in 2009 and phased out in 2019. SLMP's overall goal was to "reduce land degradation and improve land productivity in selected watersheds in targeted regions in Ethiopia" [10]. It's first component, watershed and landscape management, aimed at reforesting and afforesting degraded communal land, increasing agricultural and livestock productivity, reducing carbon emission, building climate resilience and increasing water availability. To achieve these goals, activities such as hillside communal land treatment, including the prohibition of free grazing, gully rehabilitation and cropland treatment using biophysical measures, promoting agro-forestry and fodder production, and the construction of water harvesting structures were supported in the watersheds [10].

The 21 major watersheds each comprise between 6 to 20 micro-watersheds (Figure 1) with a total of 314 micro-watersheds and an average area of 7 km2. 220 micro-watersheds (1541 km2) received SLMP support from 2011 to 2019 by KfW with the technical assistance of the German Agency for International Coorporation (GIZ). In the following, the supported micro-watersheds are referred to as treatment areas while the remaining 94 microwatersheds (543 km2) that did not receive any support are referred to as control areas.

**Figure 1.** Overview map of the location of the study areas (i.e., major watersheds) and elevation.

#### *2.2. Data*

Landsat Collection 1 atmospherically corrected Surface Reflectance (SR) Tier 1 products were used for the period 2001–2019 including three different sensor systems: Landsat 5 Thematic Mapper (TM) for the epoch 2001–2012, Landsat 7 Enhanced Thematic Mapper Plus (ETM+) for the epoch 2001–2019 and Landsat 8 Operational Land Imager (OLI) for the epoch 2013–2019. Due to the failure of the Landsat-7 ETM+ Scan Line Corrector (SLC) in 2003, the ETM+ data are reduced by about 22% in each scene [42]. Rainfall estimates were derived from Climate Hazards group Infrared Precipitation with Stations (CHIRPS) data. The product is resampled to a spatial resolution of 0.05 degrees [43]. Following Funk et al. CHIRPS data have been largely used to examine rainfall trends and drought patterns in Ethiopia. The dataset is affected by uncertainties due to the inverse distance weighting function that is used for the blending procedure [43].

Polygon shapefiles for the micro-watersheds were provided by GIZ. Furthermore, between 2012 and 2018, the GFA Consulting Group collected georeferenced location data of soil and water conservation (SWC) measures that were implemented through SLMP to monitor the progress in the treated micro-watersheds. These data were used to relate vegetation development directly to SWC measure locations. The types of measures included in the dataset represent combinations of physical SWC constructions and biological activities (e.g., planting). Personal communication with SWC experts and project leaders during field visits in 2019 revealed that two types of measures included in the dataset, hillside terraces and check dams, should have a direct impact on the surrounding vegetation cover, within a radius of approx. 500–1500 m. See Appendix A, Table A1 for details regarding the purpose and the number of geolocations available in the dataset.

#### *2.3. Methods*

All remote sensing data were acquired in Google Earth Engine (GEE), a cloud-based platform for geospatial data processing that stores a large repository of publicly available data [44]. Data processing and analysis were conducted using Python packages such as NumPy and Rasterio and the GEE client library (Figure 2).

**Figure 2.** Workflow of the overall steps of the methodology.

#### 2.3.1. Pre-Processing

For both cloud and cloud shadow masking, the C Function of Mask (CFMask), provided by USGS as pixel quality band as part of the Landsat products [45] was utilised. Visual inspection showed that using only CFMask proved to be insufficient. Therefore, clouds were additionally processed by the Google cloud score algorithm and cloud shadows by the Temporal Dark Outlier Mask (TDOM) [46]. Remote sensing analyses that integrate different sensor data require cross-calibration of the different datasets to ensure consistency [47–49]. As OLI spectral bands widths are narrower compared to ETM/ETM+, OLI NDVI values are on average higher [47] and it was therefore necessary to adjust OLI to ETM/ETM+ NDVI values before temporally aggregating the data. Transformation functions were developed using ordinary least squares (OLS) regression:

$$NDVIETM + = a \times NDVIIOLI + b \tag{1}$$

For this purpose, OLI and ETM+ NDVI images were paired based on the closest available dates. As OLI and ETM+ share the same orbit offset by 8 days, the western and eastern side of a sensor acquisition are overlapped by the eastern and western sides, respectively of the other sensor [47]. In these areas where the sensor acquisitions spatially overlap, image pairs are available with a temporal separation of only one day. To optimise

the inter-comparison, pixels used to produce the OLS models (Figure A1) were therefore extracted for the study area where overlap existed. To minimise problems related to changing surface states and conditions such as different crop cycle stages [47], pairing was done during the dry season (November to February), i.e., for two seasons, 2013/14 and 2017/18, for an area including a large range of vegetation densities (including evergreen vegetation).

NDVI has been used extensively for analysis of dryland vegetation [50] as NDVI saturation that can occur in densely vegetated areas is rarely a concern in drylands [51]. NDVI was here used as a proxy for the vegetation condition during the period of analysis. Annual NDVI maximum value composites (MVC) were produced based on all cloud-free available pixels during the period August through October, which represents in the end of the growing season. As image coverage was insufficient, a three-year (+1/−1 year) moving maximum window was used to fill data gaps (Figure A2). Annual rainfall composites were computed using the seasonal rainfall sum (March to September). As agriculture is primarily rainfed and crop productivity highly dependent on rainfall (MOFED, 2002) [52], the rain-use efficiency (RUE) was used as a proxy for assessing non-climate related changes in vegetation conditions [25,26] inherent to the implementation of the SLMP activities to counter land degradation. As vegetation productivity in the study area is predominantly determined by seasonal rainfall, RUE was calculated as the ratio of the maximum NDVI, as an approximation of ANPP, and the seasonal rainfall sum. For this purpose, the rainfall composites were resampled from 0.05 degrees to Landsat resolution of 30 m using bicubic interpolation.

#### 2.3.2. Theil-Sen Regression and Mann-Kendall (MK) Trend Test

The temporal development of NDVI was used as an indicator of land degradation ('browning') and land recovery ('greening'). Spatiotemporal patterns of NDVI and rainfall were examined using the non-parametric Theil-Sen median slope [53] to analyse changes in vegetation condition (climate-related vs. human-induced). The Mann–Kendall (MK) test [54] was applied to evaluate NDVI trends at the 99% (*p* < 0.01) and rainfall trends at the 95% (*p* < 0.05) significance level. A stricter criterion (*p* < 0.01) for NDVI trends was applied due to the temporal smoothing of the Landsat-based NDVI time series. Pixel-wise slope differences were calculated between the two periods of 2002–2010 and 2011–2018. The split of the two periods is determined by the timing of the SLMP implementation and length of time series.

Aggregated NDVI trends were calculated for treatment and control areas using the median. This was done for the total study area ("regional" scale) as well as for each major watershed that comprises both treatment and control micro-watersheds ("local" scale). To assess the effects of treatments, i.e., SLMP interventions, the significance of the difference in the distribution of per-pixel-NDVI trends between the two sample groups was evaluated using the Mann-Whitney U test, which is a non-parametric test for independent samples, non-normally distributed data and different sample sizes [55]. Spatial patterns were assessed by inspecting trend maps with the additional use of multi-temporal Google Earth VHR images.

In order to investigate the spatiotemporal relationship between vegetation and rainfall, the agreement of NDVI and rainfall trend directions was computed on a pixel-basis for each period. This method is based on a model which, following Horion et al. [26], interprets the nature of changes in ecosystem functioning based on the combination of growing season vegetation and rainfall trends. A decrease in growing season vegetation despite an increase in precipitation or vice versa is likely to be caused by human activities. In contrast, trend combinations with the same direction of change are likely to be caused by climate (see Figure A3 for further explanation). The stronger the magnitudes of change in a given direction, the more likely the cause attribution [26].

#### 2.3.3. LandTrendr

In addition to the trend analysis based on fixed time periods (Section 2.3.2), we conducted an analysis based on the LandTrendr approach in GEE [56]. This was done to further explore the importance of analysing land degradation and SLMP interventions at the level of Landsat pixel resolution rather than at coarser spatial resolution (e.g., MODIS predominantly being used for time series analysis) where subtle vegetation changes at local scale is likely to go unnoticed. The implementation of SLMP interventions was conducted at different times during the eight-year epoch 2011–2018 in the different microwatersheds. Therefore, human-induced changes through SLMP occurred presumably within a time period shorter than eight years. Horion et al. argue that abrupt changes in RUE can indicate significant changes in ecosystem response to precipitation through human activities [26]. To identify trends with a shorter duration than eight years and to obtain more information about the types and timing of the changes, LandTrendr was applied to NDVI and RUE composites. The algorithm was used to fit a model for the period 2002–2018 with a maximum number of three segments (for the sake of simplicity) and a confidence interval of 95% (*p* ≤ 0.05) determining the significance of the fitted segments. For each significant segment, the algorithm returns the magnitude, duration and rate of change as well as the start year in which the segment was detected, the end year and the corresponding NDVI value defined by the identified vertices (Figure A4). The Pearson's correlation coefficient (r) was used to mask pixels where RUE correlated with rainfall over the overall period 2002–2018 using a confidence interval of 95% (*p* < 0.05). This was done, as the use of satellite-based RUE time series to identify non-precipitation related land degradation/recovery is problematic for pixels where RUE remains correlated with NDVI, as this suggests NDVI changes still to be controlled by changes in precipitation [25,57,58].

#### 2.3.4. Effect of Soil and Water Conservation (SWC) Measures on Vegetation Trends

OLS regression was used to estimate the effect of SWC measures on vegetation trends where the trend represents the dependent variable and the distance to the SWC points the explanatory variable. To investigate the influence of SWC distance on the trends, buffers were created based on the geolocations of the different SWC types and implemented based on three different sizes: buffers with a 250 m and 500 m radius, both with a zone width of 50 m (Figure A5), and buffers with a 1000 m radius and a zone width of 100 m. Within each zone the trend results were aggregated in two ways; first, by using the median of all Theil-Sen trends and second, by using the proportion of the significant Theil-Sen as well as LandTrendr increases. The aggregated results were regressed against the corresponding distance value of the zone.

#### **3. Results**

MK trend test revealed that rainfall trends were not significant in any of the subperiods apart from 0.3% of the study area with increasing trends in 2002–2010 (*p* < 0.05). Over the entire study period 2002–2018, monotonic increasing trends were found for 17% and decreasing trends for 1.3% of the total area. When considering all rainfall trends, increasing trends were observed in 91% of the entire study area for the sub-period 2002– 2010, 54% for the sub-period 2011–2018 and 78% for the entire study period 2002–2018.

The increasing rainfall trends during 2002–2010 had a median annual change of 15.2 mm and occurred with mainly decreasing NDVI trends (Figure 3a). During 2011– 2018 rainfall experienced almost no or little changes (median annual change of 1.5 mm) while NDVI was predominately increasing (Figure 3b). Over the entire study period 2002– 2018, minimal increases in rainfall with a median annual change of 4.2 mm occurred with increases in NDVI (Figure 3c). When considering only significant rainfall trends for this period (Figure 3d), the agreement shows a more pronounced pattern of positive significant monotonic NDVI and rainfall increases with an annual change rate of 15.8 mm per year.

**Figure 3.** Spatial agreement of all annual rainfall trends and significant NDVI trends (*p* < 0.01) as density plots for the total study area for (**a**) the period 2002–2010, (**b**) 2011–2018 and (**c**) 2002–2018, and for (**d**) 2002–2018 with significant rainfall trends (*p* < 0.05).

#### *3.1. Treatment and Control Areas* 3.1.1. Theil-Sen Trends

In the period 2002–2010 the median trend in NDVI was negative with an annual decrease of −0.0065, while in 2011–2018 NDVI annually increased by 0.009 (Figure 4a). The proportion of pixels where significant NDVI trends occurred in both periods accounted for 2.6% (treatment) and 2.8% (control). The distributions of the trend differences of these pixels are left-skewed for both treatment and control areas and have a median of 0.019 (treatment) and 0.018 (control) (Figure 4b).

The most common trend was a shift from negative to positive (Figure 4c). The proportion of significant negative-positive trends was slightly higher for treatment areas (1.84%) than for control areas (1.79%) (Figure 4c). The overall trends (all trend types included), the negative-positive trends, and trends with the same sign in both sub-periods all had a positive median trend (Figure 4d). The median trends were similar for both sample groups.

The results of the individual watersheds reveal that cases where treatment areas have larger trends than control areas dominated (Table A2). When treatment and control microwatersheds from all major watersheds were treated as two large sample groups, the overall NDVI trend was larger for treatment (0.019) than for control (0.0183).

**Figure 4.** Theil-Sen NDVI trends. (**a**) Distribution of statistically significant (*p* < 0.01) trends for both sub-periods for the total study area; (**b**) Distribution of NDVI trend differences (slope of 2011–2018 minus slope of 2002–2010) for treatment and control areas; (**c**) Proportion of each type of change for all significant pixels and (**d**) Median of NDVI trend differences (slope of 2011–2018 minus slope of 2002–2010) for the overall trend and for each type of change constellation.

#### 3.1.2. LandTrendr

A LandTrendr-based approach was used to refine the Theil-Sen trend analysis that was based on two fixed time periods defined by the major scheme of SLMP support. LandTrendr detected for 46% of both treatment and control areas statistically significant

(*p* ≤ 0.05) changes in NDVI for the entire study period 2002–2018, that were subsequently analysed as a function of change types. For 92% of the study area pixels showed negative correlation of RUE and rainfall and were therefore not included in the change type analysis (Section 2.3.3). The remaining significant RUE changes (that were not correlated with rainfall changes) accounted for 2% of treatment and 3% of control areas.

The most frequent NDVI change type (amongst the 46% of statistically significant pixels) consisted of one decreasing followed by an increasing trend for 29 and 28% of all significant pixels in treatment and control areas, respectively (Figure 5a). The second most frequent change type were two consecutive increasing trends (24 and 26%, respectively). For RUE, the most common trends were of the same type as for NDVI, however with a larger share of two consecutive increasing trends (44 and 42%, respectively).

To examine the timing of detected trend segments for treatment and control areas, the onset of the decreasing and increasing trend segments with the greatest rate were extracted from each pixel (i.e., from each trend sequence) and aggregated in two groups, respectively. Generally, treatment and control areas showed the same pattern in the timing of trends without any pronounced differences. For both sample groups, the timing of the onset of approx. 90% of both NDVI and RUE decreasing trends was in 2002 (Figure 5b). In contrast, only approx. 40% of the NDVI increasing trends for both groups occurred in this year. For RUE increasing trends the proportions amounted to 53% (treatment) and 51% (control). The remaining NDVI and RUE recoveries occurred proportionally more equally distributed over the years after 2002 (2–10% in the remaining years of the time series).

**Figure 5.** LandTrendr results for treatment and control areas, respectively; (**a**) Change types detected by LandTrendr with "+" indicating increasing and "-"decreasing trends. As the change types resulted from fitting either one, two or three segments into the time series, different combinations of trends existed. For example, while "-" indicates a decreasing trend over the entire study period, "- + +" indicates one decreasing trend followed by two increasing trends; (**b**) Timing of NDVI and RUE trends.

#### *3.2. Visual Inspections of Trends Using Google Earth*

Google Earth imagery was used to inspect all 21 major watersheds and showed that the limited amount of significant RUE trends to a wide extent occurred in areas where tree plantations have been established extensively. The watershed of Banja was selected as an exemplary case of illustration (Figure 6) to link interventions on the ground with the trends captured by the remote sensing time series data. Here, the positive trend differences represent pixels where trees have been planted between 2011 and 2018 and negative trend differences were mostly found in croplands, while grazing land showed no change in NDVI (Figure 6a). NDVI recoveries (an increasing trend segment) detected by LandTrendr were observed to be more dispersed in time (Figure 6c) and space (Figure 6d) as compared the RUE (Figure 6e). The higher uniformity of the RUE results at the pixel level can be explained by the coarser resolution of the CHIRPS rainfall data used to compute RUE time series at the scale of Landsat NDVI data. The strongest changes in RUE were therefore more likely to be found in the same year for different pixels as compared to changes in NDVI. RUE recoveries were mostly detected between 2013 and 2014 (Figure 6c) and represent pixels where trees have been grown on former cropland or grazing land (Figure 6e). This is exemplified in the time series (Figure 6b) at the point of interest (POI) where LandTrendr detected the strongest NDVI recovery in 2012 that continued until 2017 with a rate of 0.07 per year. The rather high NDVI and low rainfall in 2015 led to high RUE in this year; hence the detection of a positive RUE trend in 2013 with a rate of 0.09 per year. The VHR images show that while agricultural fields and grazing lands were largely unterraced in 2005 (Figure 6f), benches of terraces were fully established in 2013 (Figure 6g). From this year on, tree cover expanded until 2020 (Figure 6h).

**Figure 6.** Watershed of Banja. (**a**) NDVI trend differences (slope of 2011–2018 minus slope of 2002–2010); (**b**) Time series at the POI; (**c**) Timing of the largest NDVI and RUE recoveries in the AOI; (**d**) Spatial distribution of the timing of the largest NDVI recoveries and of (**e**) RUE recoveries in the AOI; VHR images showing (**f**) mainly unterraced hillside in 2005, (**g**) the establishment of terraces in 2013 and (**h**) expanded tree cover area in 2020.

Furthermore, an inspection of Google Earth imagery from the watershed Yilmana Densa (Figure 7a) revealed that positive vegetation trends occurred mainly outside agricultural fields. Generally, the timing of NDVI recoveries in this watershed was equally distributed over the entire study period while RUE recoveries occurred mainly in 2009, 2011 and 2013 (Figure 7b). Recovery trends detected after 2010 were mainly detected along riverbanks, on steep hillsides, and on land affected by gully erosion. This is shown in area I (Figure 7c) and area II (Figure 7d) which accordingly showed positive trend differences that were associated with a trend shift of significant monotonic negative-positive trends. The two VHR images of the river in area I show that this development is related to a decline in eroded soil from 2013 to 2016. For area II, the VHR images show that while the hillside

was largely covered by degraded soil in 2014, it was revegetated in 2019 which explains the positive NDVI trends. Furthermore, area III (Figure 7e), characterised by degraded land and affected by gully formation, increased in vegetation cover between 2005 and 2019 which is reflected by the positive NDVI trend differences. This result also coincides with the evaluation through the performance assessment by GFA Consulting Group in which the area of the corresponding micro-watershed was described as "reclaimed land from gully erosion transformed into forage production and other economic activities" (personal communication with GFA's SLMP team leader).

**Figure 7.** Watershed of Yilmana Densa. (**a**) NDVI trend differences (slope of 2011–2018 minus slope of 2002–2010); (**b**) Timing of the largest NDVI and RUE recoveries; Zoom-in areas showing (**c**) a decline of eroded soil at a riverbank (area I), (**d**) an increase in vegetation along a hillside (area II) and (**e**) an increase in vegetation in and nearby a gully (area III) with corresponding trend patterns.

#### *3.3. Effect of SWC Measures on Vegetation Trends*

The regional regression results (results based on the available SWC data in the entire study area) showed strong negative relationships of trends and distance particularly for check dams. The strongest relationships were found between the distance and the median trend differences (slope of 2011–2018 minus slope of 2002–2010) within 250 m (r = −0.97 \*\*, <sup>r</sup><sup>2</sup> = 0.94) and 500 m (r = −0.98 \*\*\*, r<sup>2</sup> = 0.8) (Figure A6a, Table 1) as well as between the distance and the proportion of significant increases 2011–2018 within 250 m (r = −0.97 \*\*, <sup>r</sup><sup>2</sup> = 0.93) and 500 m (r = −0.95 \*\*\*, r<sup>2</sup> = 0.91) (Figure A6b, Table 2). For the latter type of trends, terraces showed strong negative relationships within 500 m (r = 0.87 \*\*, r2 = 0.76) and 1500 m (r = 0.85 \*\*\*, r2 = 0.73) radius (Table 2). Otherwise, linear relationships within 1500 m buffers were for both SWC types predominately weak.


**Table 1.** Regional OLS regression results for median trend differences (all trends).

**Table 2.** Regional OLS regression results for the proportion of significant positive trends in the period 2011–2018.


Differences exist with different significance levels (\*\* *p* < 0.01, \*\*\* *p* < 0.001).

The results showed stronger correlations for a few local regressions models (results of the individual watersheds). An example is the watershed of Tahtay Koraro that experienced increasing trends particularly along the micro-watershed borders where steep slopes exist. Strong negative relationships of trends and the distance to both check dams and terraces were observed. Significant increases between 2011 and 2018 clustered at the location of two check dams and one terrace (Figure 8a). The strongest relationship was found within a 250 m radius (r2 = 0.98) with a decreasing proportion of significant trends from 74% at 50 m, to 36% at 250 m and to 13% at 1500 m (Figure 8b). The VHR images reveal that the hillsides were, in accordance with the SWC data, terraced in 2016 and increased in vegetation cover up to 2019 (Figure 8c).

**Figure 8.** OLS regression of trends and the locations of SWC measures in the watershed of Tahtay Koraro. (**a**) Significant NDVI trends 2011–2018 with SWC points showing the year of completion; (**b**) Combined regression model for all check dams and terraces in the watershed for the proportion of significant increasing NDVI trends 2011–2018; (**c**) Multi-temporal Google Earth VHR images showing the zoom-in area with an increase in vegetation cover from 2004 to 2019 and the construction of terraces in 2016. \*\* *p* < 0.01, \*\*\* *p* < 0.001.

Another example of spatial correlations of increasing NDVI trends and SWC measure locations was observed for the watershed of Gudeyabila where the trend differences (Figure 9a) and recovery trends detected after 2010 (Figure 9b) clustered in close proximity to a check dam. Here, a decrease of the median NDVI from 0.015 at 50 m to 0.005 at 250 m (r<sup>2</sup> = 0.85) was observed (Figure 9c). The VHR images show the existence of a gully as well as mainly unterraced hillside in 2013 (Figure 9d). The construction of the check dam, which was completed at the gully in 2015, was accompanied by treatment through terraces of the surrounding hillside area including grazing land and cropland.

**Figure 9.** OLS regression of trends at a check dam location in the watershed of Gudeyabila; (**a**) Trend differences (slope of 2011–2018 minus slope of 2002–2010); (**b**) Timing of recovery trends; (**c**) Regression model of the median trend difference for the check dam location; (**d**) Google Earth VHR images showing mainly unterraced hillside in 2013 and the construction of terracing in 2015. \* *p* < 0.05, \*\* *p* < 0.01.

#### **4. Discussion**

During the entire period 2002–2018 and during the second sub-period 2011–2018, the study area showed more pixels of NDVI increase than decrease and vice-versa during the sub-period 2002–2010 (Figures 3 and 4a). The browning trends in the latter period coincide with previous findings by Hermans-Neumann et al. who identified declining net primary production between 2000 and 2009 in the highlands of Amhara and Oromia [18]. Besides the greening trends between 2002 and 2018, several results of this study indicate a shift from browning to greening within the same period. On the one hand, this was indicated by the results from the Theil-Sen trend analysis with a shift from decreases to increases in 2011 (Figure 4c), on the other hand, browning to greening was also expressed by the LandTrendr results showing a relatively large number of the negative-positive change types (Figure 5a). For RUE, this change type was the second most common one. The general shift from negative to positive trends was also in agreement with the timing of increases and decreases, as decreases occurred mainly in 2002 and increases mainly after 2010 (Figure 5b).

14% of the study area experienced a significant increase in both NDVI and rainfall in the overall study period 2002–2018. This spatiotemporal pattern indicates climate-related greening [26]. In the first epoch 2002–2010, the study area presented mainly negative rates of change in NDVI despite increasing (but not significant) rainfall trends (Figure 3a). Even though this pattern indicates human-induced browning, evidence of this is not provided due to the insignificance of the rainfall trends and the results rather suggest that declining NDVI was not related to a long-term change in rainfall. Rainfall variability in North, West and central Ethiopia increased during the period 1983 to 2012 [59]. As extreme weather conditions pose major challenges to agricultural activities through reduced crop yields and intensified soil erosion [59], changing intra-annual patterns are overlooked by using seasonally summed rainfall and the results could possibly be improved by examining trends in rainfall intra-annual variability [60].

#### *4.1. Treatment and Control Areas*

The Mann–Whitney U test results showed larger median NDVI trends between 2011 and 2018 for treatment areas than for control areas at the regional level with most local test results (i.e., within each individual major watershed) indicating this as well (Table A2). However, the medians differed only marginally and the significant differences in the distributions of the per-pixel trends were not immediately apparent in the trend maps. The significance of test results may be attributed to large sample sizes facilitating distributions to be significant.

The LandTrendr results did not reveal any substantial differences in the types and timing of changes between treatment and control areas (Figure 5); hence did not provide evidence for an improved development of treatment than control areas. However, it should be considered that this result may be tied to the applied methodology, particularly the use of the maximum NDVI which may not fully capture effects of interventions. This may be the case if the impact of interventions show an increase in crop yield that might be reflected better as an increase in the integral of the phenological crop cycle curve (integrated NDVI) or if the impact is the reduction of eroded soil in the end of the rainy season which could lead to higher NDVI values in the end of the phenological cycle, rather than causing an increase in maximum NDVI. Future research should, provided that the data availability will be sufficient (e.g., Sentinel-2), therefore include seasonal metrics when examining changes in vegetation condition related to SLM interventions.

#### *4.2. Visual Inspections of Trends Using Google Earth*

Whereas the previously discussed results based on the analysis conducted and reported at the level of watersheds did not indicate strong differences in vegetation development between treatment and control areas, closer visual examination of the trend maps showed that human-induced land improvements can be detected from the Landsat-based approach developed, though at localised scales rather than consistently spread throughout entire watersheds. For a few major watersheds, particularly for Banja and Sekela in Amahara region, the VHR images supported that RUE increases between 2012 and 2014 were due to the establishment of tree plantations (Figure 6). Following Mekonnen et al., communal and private lands in rural areas in Amhara region have been extensively used for the expansion of Eucalyptus and Acia plantations due to the demand for wood resource [61]. Furthermore, the government has been promoting tree planting widely through the introduction of campaigns. Consequently, woodlots, home gardens, trees on cropland and farm boundary plantations have become common agroforestry practices [61] which may explain the general emergence of small tree patches in various watersheds.

Furthermore, from visual inspection it was observed that notable changes occurred outside cultivated fields as patches of increasing trends in vegetation cover, particularly at hillside locations and along streams and gullies (Figure 7). These clusters of pixels were characterised by a negative-positive trend shift (MK) and were in agreement with LandTrendr recoveries of which the largest were detected after 2010. Since in these instances VHR images could verify the regrowth of vegetation, these negative-positive trend shifts can be linked to anthropogenic land improvement where degradation previously occurred. Moreover, findings of positive development could be confirmed by results from the GFA Consulting Group performance assessment (Section 3.2).

#### *4.3. Effect of SWC Measures on Vegetation Trends*

The OLS regression results demonstrated that positive changes in NDVI could be attributed to the impact of SLMP infrastructure. Visual inspection showed that this was visible in the spatial patterns for several locations at check dams and terraces, as shown in the case examples (Figures 8 and 9). This coincides with results from a study by Ali et al. who used NDVI and other satellite derived indices to evaluate the impact of SWC measures on different land use (i.e., on cultivated land and non-cultivated land such as degraded hillsides), and found that biophysical measures had a particular high impact on non-cultivated land [13].

In general, the most significant positive changes in NDVI were observed within the smallest buffer zone; i.e., the decrease in the medium change rate or density of change occurrence was mainly observed within 250 m distance from the location, in few cases within 500 m (Tables 1 and 2). At distances larger than this and up to 1500 m, trends levelled off or even increased again. Field visits showed that gullies typically pass through cultivated fields and in these cases revegetation efforts are conducted mainly directly at the gully only. Hence, for land restoration of gullies dense vegetation regrowth does not often take place across larger distances from the restoration activities. Land improvement in proximity to check dams could therefore also occur as a line type pattern in the trend map. In this case, circular buffer zones will not help explaining vegetation trends. Apart from this, it should be considered that rehabilitation activities can occur at a smaller scale than Landsat's spatial resolution of 30 × 30 m, in which cases changes would not easily be detected. The spatial scale of the impact of SWC interventions as implemented in Ethiopia therefore underlines the challenge of detecting changes related to improved land management in developing countries based on the use of traditional remote sensing methods for change detection.

#### **5. Conclusions**

The aim of this study was to examine vegetation dynamics between 2002 to 2018 in degraded areas in the Ethiopian highlands and assess the impact of SLMP interventions. To examine vegetation dynamics in complex landscapes in Ethiopia on a detailed spatial scale, we investigated the potential of combining remote sensing data from different Landsat sensors using cloud-based geospatial processing supporting a high-resolution time series analysis. The vegetation dynamics in the study areas showed a shift from

browning (2002–2010) to greening (2011–2018) along with an overall greening trend over the full period (2002–2018). From the spatiotemporal patterns of NDVI and rainfall it could be concluded that the browning trend was not explained by long-term changes in rainfall. In contrast, the greening trend over the full period could—for 14% of the study area—be explained by increases in rainfall. Overall, no clear patterns of anthropogenic induced changes in vegetation were found when aggregating results at the catchment scale, as NDVI median trends did not clearly indicate better development in SLMP intervention areas than in control areas. Visual inspection based on multi-temporal Google Earth imagery showed that the changes in NDVI and rain-use efficiency did spatially overlap areas of small-scale land improvements related to human management, however, on a smaller scale than a micro-watershed (the smallest aggregation level). The OLS regression results provided evidence of land recovery that could be attributed particularly to SLMP infrastructure (check dams and terraces). Positive impacts on vegetation were found to be contributing to improving the rehabilitation of degraded hillside areas and gullies. These findings underline that the little differences found between treatment and control areas when aggregated to the level of (micro-)watersheds are rooted in a scale issue, and

highlight the need for per-pixel trend analysis using sensor systems like Landsat, or higher spatial resolution, to be able to remotely capture the effect of SLMP interventions.

The ecological improvements through SLMP, identified here at the per-pixel level from the use of Landsat time series, are an important contribution to restore terrestrial ecosystems as targeted in the Sustainable Development Goals. Continuous efforts in developing means for improved monitoring of human-induced vegetation restoration of degraded lands will be essential to maintain rehabilitated land, prevent further land degradation and support environmental sustainability.

**Author Contributions:** Conceptualization, E.B. and R.F.; methodology, E.B. and R.F.; formal analysis, E.B.; writing—original draft preparation, E.B.; writing—review and editing, E.B. and R.F.; visualization, E.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** R.F. acknowledges support by the Villum Foundation through the project 'Deep Learning and Remote Sensing for Unlocking Global Ecosystem Resource Dynamics' (DeReEco).

**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.

**Acknowledgments:** The authors would like to thank GFA Consulting Group GmbH who provided the geolocation data on soil and water conservation measures and the German Agency for International Cooporation (GIZ) GmbH who provided the data on treatment and control areas.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A**

**Table A1.** Description of soil and water conservation (SWC) measure data.


**Figure A1.** Landsat cross-calibration. (**a**) Distributions of ETM+ and OLI NDVI values; (**b**) OLS regression of ETM+ NDVI against OLI NDVI values. For a given OLI NDVI value, the corresponding ETM+ value is usually lower.

**Figure A2.** Interpolation. (**a**) Interpolation of a pixel's annual maximum NDVI time series (blue) through a three-year maximum moving window (green). The data gap in 2003 is interpolated; (**b**) Annual maximum value composite (MVC) with gaps produced by Landsat-7 ETM+ SLC error and clouds; (**c**) Interpolated MVC through applying a three-year maximum moving window.

**Figure A3.** Interpretation of the nature of changes in ecosystem functioning based on the spatial agreement of vegetation (NDVI) and rainfall trends. A decrease in NDVI despite an increase in precipitation or vice versa is likely to be caused by human management. Trend combinations with the same slope direction are likely to be caused by climate change [26].

**Figure A4.** LandTrendr. Fitted segments into an annual maximum NDVI time series.

**Figure A5.** Example of buffer zones used for the OLS regression analysis to examine trends at soil and water conservation measure locations.


**Table A2.** Mann-Whitney U results with the sample size N (number of significant pixels), the median NDVI trend, and U indicating whether treatment areas have significantly <sup>1</sup> larger or smaller trends than control areas. Watersheds that did not include any control micro-watersheds are not included.

<sup>1</sup> Differences exist with different significance levels (\* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001).

**Figure A6.** Regional OLS regression models for check dams with (**a**) the median of all trend differences and (**b**) the proportion of significant positive trends during 2011–2018 within the zones as dependent variable. The red line shows the linear slope when including zones up to 250 m distance, the blue line when including zones up to 500 m. \*\* *p* < 0.01, \*\*\* *p* < 0.001.

#### **References**


### *Article* **Social-Ecological Archetypes of Land Degradation in the Nigerian Guinea Savannah: Insights for Sustainable Land Management**

**Ademola A. Adenle 1,2,\* and Chinwe Ifejika Speranza <sup>1</sup>**


**Abstract:** The Nigerian Guinea Savannah is the most extensive ecoregion in Nigeria, a major food production area, and contains many biodiversity protection areas. However, there is limited understanding of the social-ecological features of its degraded lands and potential insights for sustainable land management and governance. To fill this gap, the self-organizing map method was applied to identify the archetypes of both proximal and underlying drivers of land degradation in this region. Using 12 freely available spatial datasets of drivers of land degradation—4 environmental; 3 socio-economic; and 5 land-use management practices, the identified archetypes were intersected with the Moderate-Resolution Imaging Spectroradiometer (MODIS)-derived land-degradation status of the region, and the state administrative boundaries. Nine archetypes were identified. Archetypes are dominated by: (1) protected areas; (2) very high-density population; (3) moderately high information/knowledge access; (4) low literacy levels and moderate–high poverty levels; (5) rural remoteness; (6) remoteness from a major road; (7) very high livestock density; (8) moderate poverty level and nearly level terrain; and (9) very rugged terrain and remote from a major road. Four archetypes characterized by very high-density population, moderate–high information/knowledge access, and moderate–high poverty level, as well as remoteness from a major town, were associated with 61.3% large-area degradation; and the other five archetypes, covering 38.7% of the area, were responsible for small-area degradation. While different combinations of archetypes exist in all the states, the five states of Niger (40.5%), Oyo (29.6%), Kwara (24.4%), Nassarawa (18.6%), and Ekiti (17.6%), have the largest shares of the archetypes. To deal with these archetypical features, policies and practices that address increasing population in combination with poverty reduction; and that create awareness about land degradation and promote sustainable practices and various forms of land restoration, such as tree planting, are necessary for progressing towards land-degradation neutrality in the Nigerian Guinea Savannah.

**Keywords:** archetypes; self-organizing maps; land degradation; drivers; savannah; Nigeria

### **1. Introduction**

With increasing global population, environmental change, and competing claims on land, the need to maintain land productivity and reduce land degradation has become even more critical. Various global initiatives reflecting this urgency include the United Nations Convention to Combat Desertification (UNCCD) goal of achieving a Land-Degradation Neutral (LDN) World [1]; the African Forest Landscape Restoration Initiative (AFR100) aiming to restore 100 million hectares of land in Africa by 2030 [2]; and the Great Green Wall Initiative across the Sahel [3]. Yet, land degradation (LD), i.e., the persistent reduction (negative trend) or loss of the biological productivity or ecological integrity of land or its value to humans, remains a diverse and complex issue [4,5].

The African Savannah is among the globally threatened landscapes [6], where climatic and edaphic conditions, as well as human activities, constrain vegetation regeneration [7,8].

**Citation:** Adenle, A.A.; Ifejika Speranza, C. Social-Ecological Archetypes of Land Degradation in the Nigerian Guinea Savannah: Insights for Sustainable Land Management. *Remote Sens.* **2021**, *13*, 32. https://dx.doi.org/ 10.3390/rs13010032

Received: 17 November 2020 Accepted: 23 December 2020 Published: 23 December 2020

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

**Copyright:** © 2020 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/).

In West Africa, the savannah ecozones are prone to LD and experience both anthropogenic and non-anthropogenic pressures [7,9]. These expose its 353 million inhabitants and their livelihoods to various impacts such as decline in ecosystem services, food insecurity, migration, and civil conflict [8]. In Nigeria, LD cases across its agroecological zones are mostly triggered by a singular factor or a combination of factors, which include: desertification, deforestation leading to biomass loss, land pollution caused by oil spillage and illegal mining, as well as extensive soil erosion [10]. Thus, agroecological zones like the Nigerian Guinean Savannah (NGS) are experiencing pressures from urbanization, agricultural expansion and an increasing population that is largely dependent on land resources for livelihoods [11,12]. The NGS, covering 49% of Nigeria, is widely degraded, its ecosystem services continue to decline, and livelihoods remain precarious [11,12]. Unsustainable land use and climate stress have been implicated in this widespread degradation [11,13–15]. While the indicators of LD drivers, their interplay, and implications at scale are generally acknowledged [14,16], there is a lack of knowledge of the constellation of factors characterizing specific degraded landscapes, such as the NGS, and their interplay [11]. The explanations of several studies in the NGS on LD drivers are often without an integrative perspective capable of exposing the interactions between potential LD drivers [16,17].

Recent studies thus stress the need for an integrative approach to enable a better understanding of the constellation and interplay of LD drivers in land systems [14,17]. One such integrative approach is archetypical clustering for identifying recurrent patterns in land conditions [18,19]. Archetypes, i.e., patterns or processes that occur repeatedly across space and time [19,20], have been found to provide a more holistic understanding of land system processes [20]. This understanding enables comparability across cases and helps identify strategic policy options to address land management across scales [18–20]. Archetype studies have been conducted on food security in the Peruvian Altiplano [21]; institutional analysis and climate change in the Peruvian Andes [22]; national analysis of ecosystem services in Germany [23]; water governance of river basins [24], and global land resources management [25]. An archetype approach thus helps to illuminate the associative patterns and influence of the complex drivers of global changes (such as LD) that have often been treated in isolation [19,20,25].

This paper thus aims to identify the characteristic patterns of social-ecological factors associated with LD in the NGS and to analyze the implications of the identified LD archetypes for land governance and sustainable land management (SLM) in the region. In the subsections, the description of the social-ecological conditions in the NGS, the study methods and hypothesis, the archetypical patterns of drivers, their thematic, and spatial characterization, including the links to different levels of LD, and the implications of the archetypes for land governance and SLM in the NGS are presented.

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

#### *2.1. Study Area*

The Nigerian Guinea Savannah (NGS) (Figures 1 and 2) is an ecological region found between 6.50◦N and 9.62◦N, 2.77◦E and 13.18◦E. Occupying central Nigeria, it is the country's most extensive ecological zone, referred to as the Middle Belt of Nigeria. The zone consists of parkland savannah, gallery forests, and derived savannah, including distinctive montane vegetation [26], with tropical dry and wet seasons. Rainfall in the wet season (April to October) is about 1240–1440 mm. The dry season lasts from November to March [27]. Average maximum annual temperature varies from 31 ◦C to 35 ◦C while the average minimum ranges from 20 ◦C to 23 ◦C [27]. The region is broadly divided into two sub-regions based on distinctive vegetation types, namely the Southern and Northern Guinea Savannahs. The Southern Guinea Savannah has taller trees, such as locust bean tree (*Parkia biglobosa*), and grasses such as Gamba grass (*Andropogon gayanus*); while the Northern Guinea Savannah is characterized by bush with some trees (e.g., *Isoberlinia* spp.) and relatively shorter grasses (e.g., *Hyparrhenia* spp.) [27]. The NGS has a high level of fauna and flora and hosts major perennial rivers such as River Niger and River Benue. Sub-

sequent years of uncontrolled deforestation and poor land management has transformed the zone largely into a degraded landscape (Figure 2) [12].

**Figure 1.** Map of the Nigerian Guinea Savannah including the administrative boundaries of states (Adapted from [26]).

**Figure 2.** (**a**) Landscape showing a degraded patch of the Nigerian savannah vegetation in Shiroro local government area, Niger state, Nigeria; (**b**) a mix of degraded land with few bushes and trees interspersed in Borgu local government area, Niger state, Nigeria. (Source: own fieldwork, 2019).

#### *2.2. Framing Land-Degradation Drivers*

In this study, LD drivers were understood as determinants, i.e., reasons, factors, or actions, shaping the rate of human activities, leading to a decision to remove or reduce vegetation cover thereby causing a decline or loss of land resources' productivity [14,28]. Factors prompting LD are broadly categorized either as proximal or underlying drivers [16,17,29]. Proximate drivers are human activities or immediate actions, including the decision to directly use or alter the land cover [14]. Underlying drivers include indirect or underpinning factors, such as socio-economic factors (e.g., poverty) and biophysical factors (e.g., topographic variables) that trigger proximate causes of LD [14]. Based on literature [17,28,29] and a report on the LDN target for Nigeria in 2018 [30], three main categories of LD drivers of the Nigerian Guinea Savannah Archetypes (NGSA) namely: environmental, socioeconomic, and land-use management practices were identified. Based on the available spatial and remote sensing data, 12 drivers comprising three environmental, four socioeconomic, and five land-use management practice categories were selected (Table 1a–c). The processes through which the variables (whether as proximate or underlying drivers), drive the LD are explained in Table 1a–c.

#### *2.3. Datasets Selection*

For this study, factors were considered based on their influence on land-use decisions and LD [16,29,31] (see Table 1 for details). Climatic variables such as rainfall were not included because human-induced activities are the prevailing cause of LD, especially in the NGS [11,12]. However, human-induced drivers are better managed than climatic drivers of LD [11,32].

#### 2.3.1. Land-Use Management Practices

Land-use management practices have been linked to degraded land in Nigeria through land clearing and conversion including intensification [13]. Other unsustainable land-use management practices include rapid agricultural expansion, uncontrolled bush burning, deforestation, excessive wood extraction, unplanned infrastructure extension, and urban development, including overgrazing [12,13,33]. For this study and based on spatial data availability, data on land-use management practices were acquired and processed, respectively. The fire-occurrence density was derived from the National Aeronautics and Space Administration (NASA) active fire hotspot data of 2018 (firms.modaps.eosdis.nasa.gov) after running a spatial point-density analysis of the active fire spots at 250 m resolution. The livestock grazing intensity data for 2005, developed by Harvest Choice in 2018 (www.ifpri.org/project/harvestchoice) was used. The generated distance to the major road in 2016 for Nigeria, at a resolution of 3 arc-second (approximately 100 m at the equator) was acquired from (www.worldpop.org/project/categories?id=14) [34]. The Euclidean distance analysis of the extracted major town polygons from Google Earth (www.google.com/earth/) was used to determine the distance to major towns at 250 m resolution. Using the International Union for Conservation of Nature (IUCN) recognized protected area polygon for Nigeria (www.protectedplanet.net/country/NGA), the density of the protected area at 250 m resolution across the NGS was generated through point-density analysis.

#### 2.3.2. Socio-Economic Drivers

Demographic and socio-economic factors such as population density, income, poverty, and illiteracy are drivers of LD in certain contexts [16,29]. In other contexts, they may even be drivers of improved land conditions [35]. However, in the case of Nigeria, its over 200 million inhabitants and population density of 226 km2 place a huge demand on land resources [33,36]. Thus, gridded human population density constructed for 2018 in 2020, from random forest-based dasymetric redistribution at 3 arc-second (approximately 100 m at the equator) spatial resolution, was downloaded from www.worldpop.org [34]. Male and female literacy layers of high resolution at 1 km × 1 km gridded cells developed for 2003 in 2017, based on a geostatistics approach [37], were acquired [34]. These are necessary as they are proxies for access to agricultural extension information, as previous studies show that limited information on SLM drives LD [16,38]. The poverty headcount in percentages for Nigeria at 1 km for 2013 mapped through Bayesian model-based geostatistics analysis was downloaded from www.worldpop.org [34], because poverty can foster practices that cause LD, while LD can foster a poverty trap [38]. These drivers were selected in that they are known factors that influence decisions on land and SLM LD [28,29,38].

#### 2.3.3. Environmental Drivers

Land can be sensitive to degradation due to its environmental and physiographic characteristics, influencing human decisions to use land [31]. Such influential characteristics include soil bulk density (BD), elevation, and slope [13,39]. For this study, the alreadyprocessed NASA Shuttle Radar Topography Mission (SRTM) based slope and elevation for Nigeria by [40] at a resolution of 3 arc-second (approximately 100 m at the equator)

produced for 2020 was acquired from www.worldpop.org [34]. The spatial mapping of the bulk density (BD; soil compaction) in 2018 at 250 m resolution from 0 cm to 30 cm depth were downloaded from www.soilgrids.org and then averaged to give the overview of the BD for the study area [41–43].

#### *2.4. Methods*

#### 2.4.1. Conceptual Framework

The workflow in Figure 3 was applied to the 12 drivers (Table 1a–c), which served as inputs for identifying archetypes. The intermediate outputs from the framework, such as correlation between drivers including cluster features, can be found in supplementary material Figures S1–S5.

**Figure 3.** The conceptual framework.

2.4.2. Identifying Archetypes of Land-Degradation Drivers Using Self-Organizing Maps

To develop the Nigerian Guinea Savannah Archetypes (NGSA) of LD drivers, the 12 driver datasets were clipped to the boundary outline and resampled to a 250 m pixel size, using the nearest neighbor technique (Figure 3, Box 2), and then projected to a uniform coordinate system (i.e., Minna/UTM zone 31N). To enhance the clustering of the datasets, data were normalized (Z-score) to reformate into a common scale (Figure 3, Box 2). Correlation was calculated to identify relationships between the drivers (Figures S5 and S6) and the extent of interdependencies among the dataset that might limit the analysis. Then the Kohonen Self-Organizing Map (SOM) technique was applied to generate a singlelayer map of the archetypes of LD drivers [54]. The SOM approach involves organizing data (in this case, the 12 spatial datasets) into patterns based on their inherent similarities and dissimilarities into different groups [54,55]. By testing the different combinations of clusters and using a performance analysis involving both the Davies–Bouldin index and the mean distance of the classified grid cells [23,25], nine clusters of LD drivers were identified as archetypes (see supplementary material, Figures S2–S5). The Z-score, i.e., standardized score from the SOM, was used to examine values of each driver in terms of distance from the mean and to explain the clusters. A Z-score that is equal to the mean is zero, while positive and negative Z-score, i.e., greater than or less than zero, depict distance way above or below the mean. Through an ordinal scaling, the relative standing of the Z-score from the mean was ranked. If the Z-scores = 0, this implies mean/low influence of a driver; Z-score ≤ ±1 = moderate influence; ±1 < Z-score < ±2 = high influence; and Z-score ≥ ±2 = very high influence, respectively (Figure 3, Box 3). From the absolute values of the Z-scores of each cluster, the percentage dominance categories in the clusters were determined [23] (Figure 3, Box 3).


