Spatio-Temporal Dynamics of Grassland Loss in Sengwa Wildlife Research Area, Sebungwe, Zimbabwe ()
1. Introduction
Savannas represent a major global biome, covering nearly 60% of sub-Saharan Africa, 50% of southern Africa, and approximately 20% of the Earth’s land surface [1] [2]. The savanna biome is characteristically known to have a mixture of forbs, grasses, and woody species, with a discontinuous woody canopy and a continuous grass understory [3]. The balance between woody vegetation and open grasslands dictates habitat quality and biodiversity patterns in savanna ecosystems [4]. Grasslands and woodlands produce a range of ecosystem services, which include food and fibre production, carbon sequestration, maintenance of the genetic library, and recreation purposes [5]. Vegetation shifts from grassland to woodland diminish or completely eliminate habitat for wildlife [6], resulting in localized extirpation of specialised species [7] and loss of ecosystem services [8]. However, grassland habitat loss from bush encroachment threatens the biodiversity of savanna ecosystems [7] [9] and is one of the most critical degradation indicators [10]. Considering the negative impacts associated with grassland habitat loss, monitoring of grasslands in savannas becomes an important issue that requires special attention in the 2020 to 2030 decade of ecosystem restoration [11]. Habitat loss in terrestrial ecosystems is a key indicator in the Global Biodiversity Framework 2050 [12]. Therefore, quantifying the transitions of grassland habitat components in savanna ecosystems is a necessary conservation action and modeling future scenarios provides proactive decision-support tools for adaptive management.
Remote sensing-based land-cover change products are used to monitor spatiotemporal metrics of habitat loss at local and global scales to inform biodiversity assessments and conservation [12]-[14]. The use of high-resolution, multi-temporal transition data enables conservation managers to track loss, recovery, temporary fluctuation, and gains in habitat [15]. The application of predictive metrics, specifically annualized rate of change (ARC) and transition probability matrices (TPMs), provides a scientific basis to forecast the landscape status [13]. Past research in Sengwa Wildlife Research Area (SWRA) has highlighted significant woodland regeneration after many years of elephant population decline and woody bush encroachment [15]-[18]. Little emphasis has been placed on exploring the impact on grassland habitat. This study sought to i) determine the extent of land cover change across seven vegetation communities in Sengwa from 2001 to 2025 using multi-temporal satellite imagery, ii) identify specific transition pathways of grassland loss and potential recovery, and iii) utilize transition probability matrices to model potential future landscape scenarios under current environmental conditions.
2. Materials and Methods
2.1. Study Area
The SWRA is located at the southern end of the Chirisa Safari Area (CSA). Sengwa Wildlife Research Area, where the SWRI is located, covers an area of 373 km2 of the CSA. The area lies between 28˚03' and 28˚20'E and 18˚0' and 18˚13' in Gokwe South District, northwestern Zimbabwe (Figure 1). The research programme for Sengwa Wildlife Research Area has been the long-term monitoring of wildlife populations and vegetation, including the interaction between the two [19]. Annual precipitation ranges from 400 - 700 mm. The mean annual temperature is 22.2˚C. The area has a diverse mammalian community of 7 species of large carnivores and 18 species of large herbivores, 230 species of birds, 15 reptile species, over 50 amphibian species, and 5 fish species [20].
Figure 1. Location of SWRA in Zimbabwe.
2.2. Methods
To determine the spatial extent and structural trajectories of land cover change across the seven vegetation communities, a multi-temporal post-classification change detection analysis was executed spanning a 24-year baseline (2001-2025). Multi-temporal satellite imagery from 2001, 2013, 2020, and 2025 was acquired from the United States Geological Survey (USGS) at 30 m × 30 m spatial resolution and integrated into an ArcGIS Desktop 10.8 GIS environment. The imagery underwent rigorous pre-processing and atmospheric correction prior to supervised image classification using an Interactive Supervised Classification (ISOCLASS) algorithm to partition the landscape into seven vegetation classes, namely, Open Grassland, Mixed Miombo, Mixed Mopane, Mixed Combretum, Riverine, Bareland, and Bare Cliffs. Independent accuracy assessment was performed for each classified image using 350 labelled validation points per year, cross-tabulated against the classified map into a confusion matrix; overall accuracy (sample-based and area-adjusted, with 95% confidence intervals) and Cohen’s kappa were computed for each image. Classification accuracy was high and consistent across all years assessed (2001, 2013, 2020, 2025): sample-based overall accuracy ranged from 97.1% to 98.9%, area-adjusted overall accuracy from 98.3% to 99.7%, and kappa from 0.967 to 0.987 (Table 1). Independent classification and accuracy assessments were performed for each target year, after which a pixel-by-pixel cross-tabulation matrix was computed for each pair of consecutive years via the ArcGIS Combine geoprocessing tool to derive land-cover transitions. The direct mathematical extraction of net areal modifications (measured in hectares) to define the long-term vegetation shifts was done using Python software. Landsat 5 Thematic Mapper (TM) imagery was used for the 2001 classification, and Landsat 8 Operational Land Imager (OLI) imagery for the 2013, 2020 and 2025 classifications. All scenes were acquired in August, within the dry season, ensuring comparable, leaf-off phenological conditions across the four classification years; the 12 - 13-year intervals between epochs (2001-2013, 2013-2020, 2020-2025) also allow enough time for meaningful vegetation change to be detected against classification noise.
Table 1. Classification accuracy assessment summary for each classified image (2001, 2013, 2020, 2025).
Year |
N |
OA (Sample) |
OA (Area-adj.) |
95% CI |
Kappa |
2001 |
350 |
98.9% |
99.1% |
±1.2% |
0.987 |
2013 |
350 |
98.0% |
99.7% |
±0.2% |
0.977 |
2020 |
350 |
97.1% |
98.3% |
±2.1% |
0.967 |
2025 |
350 |
97.7% |
98.7% |
±2.0% |
0.973 |
Overall accuracy and Cohen’s kappa were computed from 350 independent validation points per year; area-adjusted overall accuracy follows Olofsson et al. (2014) and is reported with its 95% confidence interval.
2.3. Data Analysis
To model the specific landscape transition pathways and isolate the underlying ecological drivers governing the contraction or recovery of open grasslands, a programmatic data analysis pipeline was constructed using Python software. Areal datasets extracted from the multi-temporal remote sensing classifications were ingested into Python’s pandas data science library to execute high-fidelity matrix cross-tabulations and longitudinal accounting of absolute structural transformations measured in metric hectares (ha). To evaluate the structural evolution and proportional composition of the seven vegetation communities across all temporal baselines (2001, 2013, 2020, and 2025), a multi-temporal stacked diagram was compiled using Matplotlib and Seaborn to map the shifting dominance thresholds and volumetric allocations over time. Furthermore, to explicitly visualize the directional routing of landscape degradation and restoration, a pixel-by-pixel land-cover transition matrix was computed for both the contraction window (2001-2013) and the recovery window (2013-2025). This quantitative framework was then translated into a topological Sankey flow diagram utilizing the plotly.graph_objects library. The resulting Sankey model mapped the explicit trajectories of land-cover conversions, providing visual and mathematical proof of specific ecological processes, such as the degradation flow of open grassland into bare soil or its replacement via mixed Combretum bush encroachment. This integrated Python-driven computational framework provides a rigorous, reproducible method to audit and visually decouple the complex tree-to-grass equilibrium dynamics within the Sengwa ecosystem.
To simulate future landscape horizons and evaluate the long-term structural equilibrium of the Sengwa savanna ecosystem under contemporary environmental conditions, a comparative predictive modelling architecture was constructed using Python. The framework integrated two distinct mathematical approaches to forecast land-cover distributions toward the 2030 horizon: an average Annual Rate of Change (ARC) model and a stochastic Markov Chain model. The ARC model served as a deterministic baseline, executing a linear extrapolation by calculating the net annual trajectory of each vegetation class across the historical 24-year window and projecting that constant rate forward [21]. The ARC approach was adopted as an ecologically interpretable baseline, even though it cannot capture ecological thresholds or feedbacks [21] [22]. In contrast, the Markov Chain model was deployed to account for localized system resilience, structural inertia, and non-linear land-use feedbacks [23]. Operating as a discrete-time, first-order stochastic process, the Markov model assumes the future state of a land-cover pixel depends exclusively on its current state rather than its historical pathway [24]. This choice is ecologically motivated; savanna tree-grass dynamics are known to exhibit state-dependent feedbacks (e.g. fire-vegetation and herbivory-vegetation interactions) rather than a single constant historical rate [23] [24]. An empirical transition probability matrix (P) was computed from the cross-tabulated 2013-2025 pixel-based transitions (the most recent inter-survey interval, chosen so that the projection reflects the current post-recovery trajectory rather than the earlier 2001-2013 contraction phase), by dividing the number of pixels transitioning from class i to class j by the total number of pixels in class i at the start of the interval; this is mathematically represented as:
where each row element
denotes the localized probability of a vegetation community transitioning from state
to state
, satisfying the probability row-stochastic constraint where
(Turner 1987). The future landscape state vector for the year 2030
was derived by multiplying the contemporary baseline state vector
by the transition probability matrix raised to the power of the required temporal steps
, expressed as
. This dual-mode predictive pipeline allowed for a rigorous comparative analysis between rigid linear trends and matrix-driven trends.
3. Results
Multi-temporal satellite imagery land cover change mapping indicates the highest state of vegetation openness in 2001, with open grasslands maintaining a contiguous and expansive footprint across the study area, particularly anchoring the flat drainage basins and low-lying ecotones. By 2013, the landscape underwent a structural transformation, marking the lowest geographic extent of open grasslands across the entire study period. The spatial mapping indicates a fragmented and heavily encroached grassland matrix, where previous open plains were systematically replaced by closed or open woody canopy cover between 2013 and 2025. Between 2020 and 2025, a notable decrease in bare land is observed, with grasslands emerging as bare areas (Figure 2).
Figure 2. Map of land cover change analysis from 2001 to 2025.
The results show a net loss of open grasslands by 29% between 2001 and 2025, which is a net loss of 414 hectares. Riparian woodland increased by 86% during the same period. Mixed Combretum surged by 1,494 hectares between 2001 and 2013 and decreased by 285 hectares between 2013 and 2025 (Figure 3). In the recovery phase, between 2013 and 2025, grasslands recovered by 188 hectares while bare soil decreased by 152 hectares. The aggressive growth of Combretum woody trees slowed and retracted by 285 hectares. However, the Miombo woodland increased by 79 hectares. The mopane woodland loss was experienced throughout the period from 2001 to 2025, with a net loss of 2.9%. The transitions are shown in the Sankey diagram below (Figure 4).
Table 2. Vegetation changes in hectares between the periods 2001 to 2013 and 2013 to 2025.
Class |
2001 to 2013 Change |
2013 to 2025 Change |
Net Total Change |
Bare |
267 |
−152 |
115 |
Mixed Combretum |
1494 |
−285 |
1209 |
Mixed Miombo |
−1086 |
79 |
−1007 |
Mixed Mopane |
−318 |
−242 |
−560 |
Open Grassland |
−602 |
188 |
−414 |
Riverine |
218 |
335 |
553 |
Cliffs |
29 |
76 |
105 |
All percentage changes reported in the text are calculated relative to each class’s 2001 baseline area (% change = net change ÷ 2001 baseline area × 100). Cliffs were mapped separately from the six vegetation classes because they capture bare rock faces and escarpment exposures with negligible vegetation cover, screened out during classification so that these non-vegetated pixels would not be misclassified into a vegetation class. Its area is included below to complete the land-cover accounting, and its change is comparatively minor and largely driven by classification boundary effects along escarpment edges rather than vegetation dynamics.
Figure 3. Transition pathway of the vegetation classes between 2001 and 2025.
Figure 4. Open grassland to bare patches were experienced between 2001 and 2013, indicating land degradation. During the same period, open grassland also began to transition into Mopane and Combretum shrubland. From 2013 to 2025, riparian cover increased. Open grassland continued to decline sharply after 2013 and showed recovery in 2020 to 2025.
Table 3. Potential future land cover change scenarios.
Vegetation Class |
Hectares (2025) |
Annualized Change (%) |
2030 Projection (ARC) |
2030 Projection (Markov) |
Open Grassland |
1,015 |
1.72 |
1,105.43 |
1,098.71 |
Mixed Miombo |
8,574 |
0.08 |
8,607.13 |
8,609.18 |
Mixed Mopane |
18,563 |
−0.11 |
18,463.09 |
18,463.46 |
Mixed Combretum |
6,435 |
−0.36 |
6,319.85 |
6,321.29 |
Riverine |
1,196 |
2.78 |
1,371.52 |
1,345.17 |
Bare |
977 |
−1.2 |
919.87 |
922.19 |
Cliffs |
413 |
1.71 |
449.4 |
413* |
*Cliffs were excluded from the pixel-transition matrix used for the Markov projections, as they represent static, non-vegetated terrain (rock faces and escarpments) rather than a vegetation state subject to ecological transition; their 2030 Markov value is therefore held at their 2025 mapped area.
The Annualized Rate of Change (ARC, %) values are computed from the 2013-2025 interval (the post-recovery trend), not the full 2001-2025 record; this is why classes with an overall 2001-2025 net decline in Table 2, such as Open Grassland and Mixed Combretum, show a positive short-term ARC here. Open Grassland is projected to continue expanding at roughly 1.7% per year if effective annual rainfall above 600 mm is received. By 2030, it will have recovered significantly, though it may still be below its 2001 baseline. The riverine vegetation has the highest Annualized Rate of Change (+2.78%). This suggests that the Sengwa and Lutope river systems are becoming increasingly lush, potentially providing even better dry-season forage for elephants. Both bare land (−1.2% ARC) and mixed Combretum (−0.36% ARC) are projected to decline. In a Markovian sense, these areas are converting into grassland and riverine classes, indicating both natural restoration and continued localized bush encroachment (Table 3).
4. Discussion
The multi-temporal spatial analysis of land cover change in the SWRA over a 24-year baseline reveals a highly dynamic landscape characterized by non-linear structural transitions between open grasslands, bare land, and woody vegetation cover. The observed baseline trajectory challenges traditional linear views of vegetation change in African savannas. Instead, it demonstrates and reinforces that distinct climate- and disturbance-mediated cyclical phases contribute to observed loss and bush encroachment between 2001 and 2013, followed by a localized stabilization and recovery from 2013 to 2025 [23] [25]. Cumulatively, the net loss of 29% of open grasslands since 2001 and the surge in woodland underscores a long-term trend toward increased woody thickening, a phenomenon widely documented across Southern African savannas [26]. The dramatic reduction of open grasslands from 1,429 ha in 2001 to 827 ha in 2013 represents an intensive period of bush encroachment [27]. The concurrent surge of Combretum species provides structural and compositional evidence that the open savanna matrix has been systematically invaded by hardy, fire-tolerant shrubs [1]. The vegetation change dynamics from 2001 to 2013 are fundamentally tied to an imbalance in the classic sub-Saharan savanna ecological drivers: fire regimes, megaherbivore herbivory, and rainfall dynamics [28].
The initial denudation of grass cover, likely driven by recurrent drought and veld fires, may have stripped the topsoil of competitive perennial root systems [29] [30]. In African savannas, bare soil patches create an ideal micro-environment for woody seedling recruitment; without a continuous fuel load of grass biomass to carry high-intensity canopy fires, woody seedlings easily escape the veld fires [31]. Furthermore, variations in megaherbivore pressure, specifically a reduction in elephant-driven structural canopy damage, reduced the mechanical suppression of trees [32], allowing woody species to outcompete grass species and close the open canopy space [33].
Vegetation transitions also show a net loss of 2.9% in mopane and the conversion of mixed Miombo woodlands into mixed Combretum woodland, which is an interesting observation in this study. Multi-stemmed and highly resilient shrub species are dominating the landscape, altering the structural complexity of the vegetation profile [34]. The period 2020 to 2025 shows stabilisation of the grass layer, potentially due to more effective rainfall. This recovery highlights the inherent structural elasticity and resilience of savanna systems when adaptive management or natural positive feedbacks occur [35] [36]. Reduction in bare areas and increased cover in riparian habitat is further evidence of bush encroachment. This may be a desirable and undesirable development, desirable given that ground cover reduces erosion and increases forage resources for browsers [20]. Undesirable in that unpalatable woody species do not necessarily add value to veld quality in terms of forage but may provide other ecosystem services such as shade from heat, which is key in the face of global warming, as well as cover to evade predators [18].
The structural fluctuations of grassland habitat between 2001 and 2025 have implications for habitat quality and multi-taxa biodiversity patterns within SWRA. African savannas are fundamentally governed by a tree-grass equilibrium; structural changes directly dictate the distribution, behavioural patterns, and niche of wildlife [7]. The net 29% loss of open grasslands severely compresses the available foraging grounds for grazing specialists such as zebra (Equus quagga) and buffalo (Syncerus caffer) [37]. This forces large herbivore herds into smaller, fragmented grassland patches, which accelerates localized tracking, continuous over-utilization, soil compaction, and erosion, potentially creating localized degradation loops. Conversely, bush encroachment expands the feeding niches of browsers and mixed-feeders [37]. For avian communities, these transitions influence habitat quality for birds in many ways. Loss of grassland, for instance, negatively impacts birds that rely on grasslands for food, breeding, and shelter [38].
Predictive landscape modelling is an essential asset for forward-looking conservation planning [39]; however, the drastic divergence between the two 2030 horizons highlights how underlying mathematical frameworks govern spatial predictions. The 2030 ARC projection models a catastrophic collapse of open grasslands alongside an unprecedented surge in woody vegetation. This extreme projection represents a geometric artifact of linear extrapolation. Because the ARC computes a singular, net annual average across the entire 24-year timeline, its long-term average is heavily weighted down by the severe, massive structural losses incurred during the initial 2001-2013 contraction window [22]. By transferring this historical rate directly into the future, the ARC assumes that landscape degradation occurs in a vacuum, entirely blind to the distinct post-2013 recovery, ecological thresholds, or localized management interventions. While mathematically rigid, the ARC projection serves as a valuable indicator of what could occur if the historical forces of bush encroachment are completely unchecked by management.
In stark contrast, the 2030 Markov Chain projection predicts a highly stable, stationary equilibrium that closely mirrors the observed 2025 current state, maintaining a balanced distribution where woody vegetation and open grasslands stabilize [40]. The biological plausibility of the Markov model stems from its operational mechanics; it utilizes a probabilistic transition matrix that calculates the probability of land-cover conversions based on historical pixel-by-pixel changes [41]. Because it accounts for the localized system resilience and structural inertia exhibited during the 2013-2025 recovery phase, it recognizes that the contemporary tree-grass balance has achieved a degree of functional self-regulation [42]. However, a known academic limitation of standard Markovian modelling in ecology is its assumption of stationarity [43], the premise that the probabilities of land-cover transitions will remain constant into the next decade (Baker 1989). It cannot anticipate sudden external shocks, such as an unprecedented multi-year climate drought, shifts in national park funding that cease early-season controlled burning, or rapid localized spikes in mega-herbivore populations [42]. Therefore, while the Markov model represents the most ecologically sound and realistic path forward under contemporary management conditions, both models underscore that the future preservation of vital open ecosystems remains dependent on proactive, adaptive management regimes.
5. Recommendations
Based on the multi-temporal spatial dynamics and predictive landscape modeling, the findings of this study reinforce calls for regular monitoring of wildlife habitat and adaptive conservation management actions. The long-term 29% net loss of open grasslands and the corresponding surge in woody species encroachment underline the need to implement regular monitoring to halt or reverse habitat loss at local scales. Further landscape research is also needed within and outside protected areas to ascertain the rate of grassland loss and its implications for wildlife and livestock production systems.
Given the escalating climate change phenomenon and growing opportunities for using forests as carbon sinks, management actions may shift toward prioritizing woodland over small but ecologically important microhabitats such as open grassland.
Predictive landscape modelling is an essential asset for forward-looking conservation planning that should be encouraged. However, park managers should decide on the best remote sensing and Earth observation data sources as well as modeling tools to use to monitor ecological change alongside traditional field vegetation surveys.
Author Contributions
Mukululi Ndlovu: Conceptualization, Research, and supervision. Victor Badza: Formal analysis and results presentation. Christopher Hungwe: Critical review, correction, and improvement of intellectual content.