Vegetation Trends, Drought Severity and Land Use-Land Cover Change during the Growing Season in Semi-Arid Contexts

: Drought severity and impact assessments are necessary to effectively monitor droughts in semi-arid contexts. However, little is known about the inﬂuence land use-land cover (LULC) has—in terms of the differences in annual sizes and conﬁgurations—on drought effects


Introduction
Drought as a slow-onset event is increasingly an environmental hazard due to its negative impacts on natural and human systems, including livelihoods [1][2][3][4].Drought conditions are initiated by precipitation shortfall in comparison to the climatological normal in the focus context and amplified by concurrent heatwaves and extreme high-temperature events [5,6].Impacts relate to insufficient availability of water to meet human and nature's needs, seasonal moisture deficit including soil moisture and extensive evaporation resulting in a decline in vegetation greenness and health, plant mortality, reduced water levels in dams and other adverse ecological and/or socioeconomic conditions [7][8][9][10].Studies have found and are forecasting increasing drought duration, severity and frequency in different world regions [11][12][13][14][15][16].
Drylands, which by their very nature are water deficient due to limited rainfall and water supplies [17], are particularly vulnerable to droughts.In African drylands, the occurrence of droughts is not unusual.Some studies have found increasing drought events in African drylands [4,18], whereas other studies project future increases in droughts and other high-temperature events [19,20].In the arid regions of South Africa, extreme hightemperature events in parks are increasing in frequency [18].In Botswana, where this study was conducted, there is public awareness on issues relating to drought due to its negative impacts, particularly on agriculture, which is largely rainfed [21].Studies have found observable changes of varying magnitudes in rainfall, temperature trends and drought over Botswana [21][22][23][24].With future global warming, droughts are projected to increase in frequency and severity in this region based on regional climate model simulations [5,25].In the face of climate variability and change, it is increasingly important to assess and monitor droughts because of the need to adapt and minimise impacts on the social-ecological systems in African drylands.
The need to monitor and assess droughts in drylands calls for the use of Remote Sensing (RS)-based data and methods as complementary to meteorological gauge data and climatological indices due to the drawbacks of using only ground-based data, such as inadequate spatial and temporal coverage of meteorological station data.RS image data and vegetation indices are widely used in drought monitoring [26][27][28][29].With the increasing availability of free satellite image datasets of good temporal and spatial resolution, RS affords rapid and cost-effective assessment and monitoring of droughts.Drought effects on vegetation and ecosystem services were examined in the Bobirwa sub-district, Botswana using RS-based indices [30].Although the drought situation is increasingly assessed in Botswana, the use of RS in examining droughts in Botswana is still very limited.Moreover, how drought severity differs between land systems to further exacerbate impacts in this region is not clear.
Using Normalised Difference Vegetation Index (NDVI) image time series datasets, this study contributes to the understanding of the spatial and temporal variations of droughts and the effects across LULC types and change.It integrated remote sensing with spatial statistics in Geographic Information System (GIS) for the analysis and mapping of drought.Variability and vegetation trends over 18 years for the growing seasons between 2000-2001 to 2017-2018 were assessed and afterwards compared to the ongoing growing season (2020-2021).Thus, this study better captured the seasonality of vegetation growth when considering drought severity as this relates to rainfall in dryland contexts.The spatial and temporal evolution of drought severity was assessed and mapped for each year in the time series.We provided an improved methodology for the examination of drought evolution and severity by incorporating indices of LULC change.Transitions in areas of LULC change and persistence, i.e., areas where no changes occurred, are useful indices to understand the processes, especially anthropogenic, driving the observed trends in vegetation productivity and how these relate to drought severity.Most RS-based studies on drought have not evaluated the effects according to LULC types and change.For the examination of drought severity by LULC and change to be meaningful, we further considered the differences in annual LULC configurations and sizes.Considering that both configuration and size differ from year to year, we utilised the annual LULC time series for analysis in each year identified as drought-stricken.Moreover, drought impacts on land-based resources upon which much of livelihoods are dependent would either be exacerbated or ameliorated depending on how the land is put to use and the land management practices.To better demonstrate RS capabilities for assessing drought severity at a finer, sub-national scale, the assessment was conducted in 17 constituencies (a constituency is the fourth administrative level in Botswana).

Description of the Study Location
Seventeen constituencies in the Central District of Botswana (CDB) in the eastern part of the country were used as case studies (Figure 1).The CDB is the largest amongst the nine districts of Botswana both in terms of population and geographic size (Table S1 details the population and geographic area of the 17 constituencies).The district has 576,064 inhabitants (29% of Botswana's population) as of the 2011 census.With 26% of Botswana's land area, it covers an area of approximately 147,730 km 2 .With a semi-arid, hot steppe climate (Koppen's BSh classification), rains occur in the summer months with peak rainfall in January (71-142 mm).Annual average rainfall at the constituencies ranged from 321 mm to 430 mm.Temperature ranges from 32 to 39 • C and can occasionally exceed 40 • C [21].As in most parts of Botswana, the annual evaporation rate of about 2000 mm year −1 far exceeds that of the rainfall (475-525 mm year −1 ) [4].
based resources upon which much of livelihoods are dependent would either be exacerbated or ameliorated depending on how the land is put to use and the land management practices.To better demonstrate RS capabilities for assessing drought severity at a finer, sub-national scale, the assessment was conducted in 17 constituencies (a constituency is the fourth administrative level in Botswana).

Description of the Study Location
Seventeen constituencies in the Central District of Botswana (CDB) in the eastern part of the country were used as case studies (Figure 1).The CDB is the largest amongst the nine districts of Botswana both in terms of population and geographic size (Table S1 details the population and geographic area of the 17 constituencies).The district has 576,064 inhabitants (29% of Botswana's population) as of the 2011 census.With 26% of Botswana's land area, it covers an area of approximately 147,730 km 2 .With a semi-arid, hot steppe climate (Koppen's BSh classification), rains occur in the summer months with peak rainfall in January (71-142 mm).Annual average rainfall at the constituencies ranged from 321 mm to 430 mm.Temperature ranges from 32 to 39 °C and can occasionally exceed 40 °C [21].As in most parts of Botswana, the annual evaporation rate of about 2000 mm year −1 far exceeds that of the rainfall (475-525 mm year −1 ) [4].
The district is of great economic importance to Botswana, with 23% (31,634 holdings) of all traditional agricultural holdings in Botswana [31].Moreover, the majority of the mines in Botswana are located in the CDB, such as the Morupule coal mine and the diamond mines in Lerala, Orapa and Letlhakane.The district is of great economic importance to Botswana, with 23% (31,634 holdings) of all traditional agricultural holdings in Botswana [31].Moreover, the majority of the mines in Botswana are located in the CDB, such as the Morupule coal mine and the diamond mines in Lerala, Orapa and Letlhakane.

Data Sources
Variability in vegetation condition and drought severity in the CDB were examined over 18 years using the 1 km Normalised Difference Vegetation Index (NDVI) decadal (i.e., 10-day composite) image time series from the Copernicus Land Monitoring Service (https://land.copernicus.eu/global/products/ndvi,accessed on 8 January 2021).These images were made available through the European Union-African Union-funded project on Monitoring for environment and security in Africa (MESA).The MESA was implemented for the Southern African Development Community (SADC) region comprising 15 countries and included Madagascar and the Democratic Republic of Congo.The dekadal NDVI datasets from October 2000 to 2014 were derived from SPOT VGT, and data from June 2014 to 2018 are from the PROBA-V [32].These long-term, 1-km NDVI datasets from the two sensors have been pre-processed at the source to ensure compatibility and continuity [33].
The land cover datasets were from the European Space Agency (ESA-LC) Climate Change Initiative (CCI-LC v.2.0.7)ESA CCI and Copernicus Climate Change Service (C3S-LC Mv52 https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-land-cover?tab= form, accessed on 8 January 2021).These datasets are a consistent series of multi-sensor annual maps from 1992 to 2018 [34].The land categories in the ESA-LC were aggregated into six categories for compatibility with previous studies in Botswana [35]-tree-covered areas, grassland, cropland, water bodies, artificial surfaces (settlement including infrastructure) and otherland (Table S2).Tree-covered areas comprise forests and woodlands.In Botswana, forested areas are defined as comprising multi-layered tree canopies with over 10% cover, a minimum size of 0.5 ha and heights of more than 5 m [36].Grassland incorporated shrubland and other sparse vegetation, whereas the otherland category combines bareland, rock outcrops and dunes [35].

Indicators of Vegetation Variability
Indices used for measuring vegetation variability and drought severity are based on the NDVI.NDVI is widely utilised for assessing and monitoring vegetation greenness, net primary productivity, plant phenology and land degradation in natural and human systems [35].NDVI is calculated as in Equation (1), where NIR and R are the near infrared and visible red portions, respectively, of the electromagnetic spectrum.Photosynthetically active plants absorb more incident radiation in the visible red portion but reflect more in the near infrared portion [37].
Three indicators of vegetation variability utilised in this study were derived from NDVI-NDVI difference, NDVI anomaly and NDVI trends.Other derived metrics to gauge seasonal vegetation productivity include the NDVI mean, maximum and cumulative values computed for each month in the growing season within the time series.

NDVI Difference
To analyse the variability of vegetation during the vegetation growing season over the 18-year study period, the NDVI Difference (NDVI diff ) function implemented in the MESA Drought Monitoring Services (DMS) software was utilised.NDVI diff is widely used to get an indication of vegetation state over a specific period by comparing vegetation productivity between two dekads or relative to the long-term average for the same period.This indicator highlights areas where vegetation is under stress as well as those performing well.For this study, seasonal NDVI diff was calculated for every growing season (i.e., annually) as the difference between the start dekad (i.e., first 10-day period) in October (D1) in a certain year i to the end dekad (i.e., last 10-day period) in March (D3) of the following year i + 1 in the time series data of 2000 to 2018 (Equation (2)): Standardised NDVI Anomalies NDVI anomaly captures how vegetation productivity for a certain period deviates from the long-term average dynamics.It is calculated by subtracting the considered month NDVI from the month's long-term average and dividing it by the monthly standard deviation for that period.By distinguishing areas that are normal from those that are above or below normal vegetation productivity, NDVI anomaly is helpful to identify outliers, isolate the variability in the vegetation signal and consider the reviewed period within a meaningful historical context [26].

NDVI Trend
Using the NDVI time series of 2000 to 2020 as input, we computed vegetation change as trends and their significance based on the Mann-Kendall (MK) non-parametric test.Non-parametric approaches estimate trends in a time series by quantifying the rate of change in vegetation greenness for each pixel and characterises trends in the data using the median slope [38].MK tau (τ) coefficient ranges from −1 to +1 with values greater than 0 indicating a continually increasing (greening) trend, and values less than 0 indicating a continually decreasing (browning) trend [39].NDVI trends were reclassified into increase (>0), decrease (<0) and stable (0).The MK is useful to determine the significance of changes in vegetation productivity over time and is robust to outliers [40].The reclassified NDVI trends were afterwards assessed by LULC type.In conjunction with LULC types, NDVI trends are useful in detecting land areas that are potentially degraded if significant negative trends are found over time [35].

Drought Indicators Vegetation Condition Index
Drought severity during the growing season was measured using the NDVI-based Vegetation Condition Index (VCI) as in Equation ( 3) [27].VCI compares the NDVI of a given period j (NDVI j ) with the long-term minimum NDVI (NDVI ltmin ) and long-term maximum NDVI (NDVI ltmax ) computed over a 10-year time series for the same period.
VCI is particularly useful for agriculture, as it assesses changes in NDVI through time since vegetation is water-stressed due to water deficiency such as during drought.VCI for the vegetation growing season was calculated between 2000-2001 and 2017-2018 starting with the first dekad of October in the previous year to the third dekad in March of the following year.VCI is measured as a percentage with values ranging between 0 (lowest) and 100 (highest), with values equal to or below 40% considered as drought to varying degrees of severity (Table 1).The two VCI-based indicators used for characterising drought are drought intensity and drought frequency.

Drought Intensity
Drought intensities for each growing season were calculated through the 18-year study period.These are percentages of pixels (a proxy for the surface land area) whose VCI values fell within the different drought severity and non-drought categories [41].The evolution of the drought hazard was examined for each year in the time series based on the VCI.

Drought Frequency
To relate the pixel-level VCI to each of the 17 constituencies in the CDB, the zonal statistics function in GIS was used.The median VCI value within each constituency was utilised because there is the tendency for more years in drylands to have rainfall below the mean, with the median value deviating more from the mean.Due to the sensitivity of the mean to outliers, the median is a better data distribution measure to gain further insight regarding the frequency of droughts in each constituency.Drought severity was compared to the official declarations of drought years by the government.Drought frequency was computed as the count of the number of drought years in the entire time series (growing seasons 2000-2001 to 2017-2018).A year with an annual VCI value of 40 or less is identified as drought-stricken and the sum of such drought years per constituency in the entire time series amounts to the frequency of drought [26].

Land Use-Land Cover Change
In addition to vegetation trends, drought severity was examined according to LULC types.It was also examined at the constituency level, to take into cognisance the environmental and administrative basis of drought impacts, respectively.For each drought year, drought severity effects on each LULC type were assessed based on the intersection of VCI and LULC values for that particular year.Drought intensity is expressed as a percentage of the area under each LULC type affected by varying drought severity and non-drought conditions.Percentages of the surface land area affected were then normalised by the size of the LULC type for each year.Since LULC configurations and sizes vary from one year to another, the annual LULC map for each identified drought year was utilised.Only for the land change analysis aspect was the maps of the years 2000 and 2018 used for post-classification change detection.

Vegetation Variability and Trend
The spatial and temporal variations in vegetation productivity were presented based on analyses during the growing seasons, i.e., October to March, of the 18 years.Providing insights regarding each growing season throughout the entire time series, NDVI diff maps (Figure 2) depict limited vegetation productivity in the CDB in years 2000-2001, 2003-2004, 2005-2006, 2007-2009, 2013-2014     The NDVI trend was analysed as an indicator of vegetation productivity change.Areas of negative trends amounted to about 90% of the land area in the CDB (Figure 4).Decreasing NDVI trends mostly occurred in the north-east (Nkange, Shashe west, Tonota) and to the south (e.g., Sefhare-Ramokgonami, Mahalapye east).Decreasing vegetation trends imply that these areas experienced an overall decline in vegetation cover and biomass.Increasing NDVI trends (4% of CDB's land area), implying an improvement in vegetation productivity, were most pronounced in the north-west (e.g., Boteti west) and the eastern tip along the Motloutse River (Bobonong).Areas of stable vegetation productivity (6%) were mostly in the western parts around Serowe west, Shoshong, Palapye and Serowe north.The direction of these trends and the extent of land area affected are generally in line with the overall national vegetation productivity dynamics [35].The NDVI trend was analysed as an indicator of vegetation productivity change.Areas of negative trends amounted to about 90% of the land area in the CDB (Figure 4).Decreasing NDVI trends mostly occurred in the north-east (Nkange, Shashe west, Tonota) and to the south (e.g., Sefhare-Ramokgonami, Mahalapye east).Decreasing vegetation trends imply that these areas experienced an overall decline in vegetation cover and biomass.Increasing NDVI trends (4% of CDB's land area), implying an improvement in vegetation productivity, were most pronounced in the north-west (e.g., Boteti west) and the eastern tip along the Motloutse River (Bobonong).Areas of stable vegetation productivity (6%) were mostly in the western parts around Serowe west, Shoshong, Palapye and Serowe north.The direction of these trends and the extent of land area affected are generally in line with the overall national vegetation productivity dynamics [35].

Spatio-Temporal Evolution of Drought Severity during Vegetation Growing Seasons
Drought negatively impacts vegetation growth, the supply of water for nature's needs (e.g., to sustain wildlife) and human needs (e.g., for livelihoods and food security).Seasonal maps depicting spatial and temporal variations in drought severity for the CDB were produced (Figure 5a).The growing seasons of the years 2002-2004 and 2015-2016 were the worst drought periods in the entire series.The growing seasons of 2004-2007, 2012-2013 and 2016-2018 were also affected but to a lesser extent.During droughts, the most affected areas were towards the west, south and the eastern tip of the district except for 2002-2004, when the entire district was affected by drought.Land areas most affected by drought were in Mahalapye west and east, Boteti west, Shoshong, Shashe west, Palapye and Bobonong.The magnitude of the drought hazard varied between the years considered in the time series (Figure 5b).The highest percentages of extreme, severe, moderate and mild droughts were recorded during 2002-2005.These years corresponded to drought years declared by the government [31].Drought severity ranged from extreme to mild as lower than normal, erratic rainfall amounts were recorded in the CDB, similar to most parts of Botswana during these years.

Spatio-Temporal Evolution of Drought Severity during Vegetation Growing Seasons
Drought negatively impacts vegetation growth, the supply of water for nature's needs (e.g., to sustain wildlife) and human needs (e.g., for livelihoods and food security).Seasonal maps depicting spatial and temporal variations in drought severity for the CDB were produced (Figure 5a).The growing seasons of the years 2002-2004 and 2015-2016 were the worst drought periods in the entire series.The growing seasons of 2004-2007, 2012-2013 and 2016-2018 were also affected but to a lesser extent.During droughts, the most affected areas were towards the west, south and the eastern tip of the district except for 2002-2004, when the entire district was affected by drought.Land areas most affected by drought were in Mahalapye west and east, Boteti west, Shoshong, Shashe west, Palapye and Bobonong.The magnitude of the drought hazard varied between the years considered in the time series (Figure 5b).The highest percentages of extreme, severe, moderate and mild droughts were recorded during 2002-2005.These years corresponded to drought years declared by the government [31].Drought severity ranged from extreme to mild as lower than normal, erratic rainfall amounts were recorded in the CDB, similar to most parts of Botswana during these years.

Area of Land Use-Land Cover Change and Persistence
In the CDB, areas of LULC persistence, i.e., unchanged, between 2000 to 2018 amounted to 88% (130,166 km 2 ). Figure 6a depicts the spatial distribution of areas of persistence as well as losses where LULC types have been displaced, i.e., transitioned to other land uses.Figure 6b shows the share of land under each persistent and transitioned land uses between the year 2000 and 2018 as a percentage of the surface land area (Table S3

Area of Land Use-Land Cover Change and Persistence
In the CDB, areas of LULC persistence, i.e., unchanged, between 2000 to 2018 amounted to 88% (130,166 km 2 ). Figure 6a depicts the spatial distribution of areas of persistence as well as losses where LULC types have been displaced, i.e., transitioned to other land uses.Figure 6b shows the share of land under each persistent and transitioned land uses between the year 2000 and 2018 as a percentage of the surface land area (Table S3 provides the LULC transition matrix in km 2 between 2000-2018).Land transitions in the matrix are to be interpreted as 'from-to' changes, whereby a particular LULC type in 2000 (initial year) transitions to another LULC type in 2018 (target year).For example, 37% and 2% of tree-covered areas were derived from grasslands and croplands, respectively.Thus, tree-covered areas increased from 11.5% of the total land area of the CDB in 2000 to 14.5% in 2018.Other notable transitions are the expansion of artificial surfaces such as settlements, with 60% derived from grassland, 1.6% from tree-covered, 1.4% from cropland and 1.8% from otherland areas.Thus, artificial surfaces increased from 0.05% in 2000 to 0.13% in 2018.The main gains by grassland were from tree-covered areas (2.8%) and cropland (1%).The main gains by cropland were derived from tree-covered (8%) and grassland (~28%).Although cropland expanded over time (increased from 6.3% of the total land area in 2000 to 8% in 2018), it lost 2.3% through its conversion to tree-covered areas, 1% to grassland and 1.4% to artificial surfaces.

Land Change and Associated Vegetation Trends
With about 90% of the land area in the CDB experiencing negative vegetation trends, we investigated how changes in vegetation productivity (i.e., increasing, stable and decreasing) are associated with land change.The focus is on major LULC types (tree-covered area, cropland, otherland and grassland), as these made up over 95% of the study area as of 2018. Figure 7a-d depict the spatial distribution of vegetation trends found in these major LULC types.Figure 7e shows the percentage of land area by vegetation productivity change (direction and magnitude).By LULC categories as of 2018, vegetation productivity decreased in about 98% and 94% of tree-covered areas (such as forests and woodlands) and croplands, respectively, and 94% of grasslands.In other words, the majority of tree-covered, cropland and grassland areas experienced decreasing vegetation productivity as of 2018.Seventeen percent (17%) of wetlands and settlement areas, respectively, and

Land Change and Associated Vegetation Trends
With about 90% of the land area in the CDB experiencing negative vegetation trends, we investigated how changes in vegetation productivity (i.e., increasing, stable and decreasing) are associated with land change.The focus is on major LULC types (tree-covered area, cropland, otherland and grassland), as these made up over 95% of the study area as of 2018. Figure 7a-d depict the spatial distribution of vegetation trends found in these major LULC types.Figure 7e shows the percentage of land area by vegetation productivity change (direction and magnitude).By LULC categories as of 2018, vegetation productivity decreased in about 98% and 94% of tree-covered areas (such as forests and woodlands) and croplands, respectively, and 94% of grasslands.In other words, the majority of treecovered, cropland and grassland areas experienced decreasing vegetation productivity as of 2018.Seventeen percent (17%) of wetlands and settlement areas, respectively, and 12% of otherland experienced increasing trends, signifying improved vegetation productivity.

Drought Severity by Land Use-Land Cover Type
Focusing on the drought years identified earlier in the time series analysis of vegetation variability and drought severity (refer to Figures 2 and 5), we examined how drought severity differed between the major LULC types.The 2002-2003 drought-stricken growing season was used as an example because it was the worst drought experienced in the entire time series (Figure 8a-d).For example, land areas most impacted by extreme and severe droughts, respectively, in 2003-2004 are: over otherland (32%, 31%), grassland (19%, 39%) and cropland (12%, 37%).Drought intensities for areas under each LULC type for the identified drought years are shown in percentages alongside drought severity classes (Figure 8e).Between years 2000 and 2018 in the loss areas, the greatest percentage of decreasing vegetation productivity (above 90%) were found in tree-covered, settlement and cropland areas.Areas with increasing trends in loss areas are otherland (34%), wetlands (137%) and grasslands (2%).In areas of persistence, vegetation productivity declined mostly in the same LULC types as in loss areas, whereas it improved in 27% of settlement areas, 17% of wetlands, 11% of otherlands and 4% of grasslands.Minimum and maximum values of NDVI trends in areas of land loss varied between forests (−0.3, 0.21), grassland (−0.28, 0.70), cropland (−0.27, 0.24), wetland (−0.26, 0.54), settlement (−0.24, 0.02) and otherland (−0.23, 0.35).

Drought Severity by Land Use-Land Cover Type
Focusing on the drought years identified earlier in the time series analysis of vegetation variability and drought severity (refer to Figures 2 and 5), we examined how drought severity differed between the major LULC types.The 2002-2003 drought-stricken growing season was used as an example because it was the worst drought experienced in the entire time series (Figure 8a-d).For example, land areas most impacted by extreme and severe droughts, respectively, in 2003-2004 are: over otherland (32%, 31%), grassland (19%, 39%) and cropland (12%, 37%).Drought intensities for areas under each LULC type for the identified drought years are shown in percentages alongside drought severity classes (Figure 8e).

Characterising Drought in the Constituencies Drought Severity in Constituencies in Comparison with Drought Declaration
Drought and household food security vulnerability assessments are conducted annually during the mid-growing season in Botswana [42].Since assessment and interventions are conducted at local levels, we further examined the severity of the drought in the 17 constituencies (Figure 9).The drought frequency and heatmap reveal how the constituencies were affected by droughts of differing magnitudes throughout the entire time series.

Drought Severity in Constituencies in Comparison with Drought Declaration
Drought and household food security vulnerability assessments are conducted annually during the mid-growing season in Botswana [42].Since assessment and interventions are conducted at local levels, we further examined the severity of the drought in the 17 constituencies (Figure 9).The drought frequency and heatmap reveal how the constituencies were affected by droughts of differing magnitudes throughout the entire time series.9), eight were evident in the time series, whereas the other years were favourable for the CDB.

Vegetation Condition Change and Drought Severity
There was high spatial and temporal variability in vegetation productivity during the growing seasons in the 18-year study period.This is typical of dryland ecosystems which are often non-equilibrium and dynamic in response to both climatic and anthropogenic perturbations [43].Recovery of vegetation during the start of the time series (2000)(2001) after the prolonged droughts of the 1990s was evident.This finding is corroborated by [21], which noted above-normal rainfall in the CDB in the year 2000.For example, analysing rainfall amount between 1960-2015 for Palapye, the authors in [32] noted that vegetation condition improved in 1999-2000 after 649 mm rainfall was recorded that year, which was well above the long-term average of 351 mm.
Comparing vegetation productivity in the CDB during drought and non-drought years with the mean for the entire time series revealed the impacts droughts have on vegetation.For example, when compared to the mean, vegetation productivity was very limited in 2002-2003 (the worst drought episode in the time series).Drought occurrence was evident during the growing seasons of eight declared drought years in line with [21,30].FAO special alert for Southern Africa in 2015 noted the retarded growth of early-planted crops as soil moisture was very low at the beginning of the growing season in most parts of the southern African region, including Botswana [44].The drought and household food security outlook report of 2017 [45] attributed the decrease in vegetation productivity in most parts of Botswana to negative drought impacts on vegetation.
The lower than normal vegetation productivity during some of the drought-stricken growing seasons can be attributed to droughts linked to El Niño southern oscillation (ENSO).For example, the prolonged droughts in the growing seasons of 2002-2004 and 2015-2016 coincided with the El Niño years in the recent records [4].Relating the association of ENSO to drought severity during the growing season as utilised in this study, the authors of [4] found the highest statistically significant correlations in January, February and March in Botswana, whereas they found negative non-significant correlations at the start of the season in October.At the regional level in southern Africa, the authors of [46] associated droughts with anomalies of negative Standardised Precipitation Evaporation Index and positive Sea Surface Temperature.At the global level, studies have also documented the effects of ENSO on drought severity, such as [46].
Comparing the ongoing growing season (2020-2021) with the mean for the entire time series, we found above-normal vegetation productivity after mid-October 2020 until February 2021.Thus, this suggests the full recovery of vegetation productivity during this season from the impacts of the prolonged droughts in the last couple of growing seasons.However, this observation is somewhat fraught with uncertainty judging from the below-normal vegetation productivity at the start of the season.Moreover, the growing season has not ended yet.The growing season spanning the first dekad in October to the third dekad in March was chosen to align the cropping and the raining season in Botswana, which enabled the exclusion of the dry season from the drought analyses.

Vegetation Trend and Drought Severity by Land Use-Land Cover and Change
Many studies on droughts have not examined how drought effects differ between LULC types.For those that have land use incorporated, little is known of the influence of drought severity on LULC-in terms of the differences in annual sizes and configurationseither in changing and/or persistent areas.Processes driving decreasing vegetation trends, either climatic or anthropogenic, are better identified when LULC and change are incorporated in the examination of drought severity.For example, in CDB between 2000 and 2018, vegetation productivity declined in most forests, woodlands, croplands and grasslands.In land change areas, the trend of declining vegetation was equally high.In areas of persistence, the greatest percentage of improved vegetation productivity was in wetlands, settlements, and otherlands.Minimal improvement of vegetation trends in forested areas can be attributed to the overall increase in tree-covered areas during the study period.
Previous studies in dryland contexts such as Botswana and elsewhere associate the improvement in vegetation productivity partly to bush-encroachment which remote sensing vegetation indices capture as vegetation greening but bush-encroachment is undesirable in cattle-based systems [32,47,48].
Relating LULC change to land degradation, conversion of tree-covered areas to grassland, otherland and cropland, is degrading.This is because these land transitions drive the removal of vegetation cover and contribute to land degradation processes.Similarly considered as degrading land transitions are those involving the conversion of grasslands into croplands, artificial surface areas and otherlands.For example, in Palapye, the authors of [32,49] found increases in barelands and rock outcrops with limited vegetation growth because of prolonged droughts.This bareland condition due to drought-induced vegetation decline is further exacerbated by human activities, such as overgrazing.In some instances, such as in Bobonong, researchers [30] found increases in bareland patches in communal grazing areas.Despite declining vegetation conditions during a four-year prolonged drought (2002)(2003)(2004)(2005), livestock overgrazed and natural pastures degraded, as pastoralists had no incentive to destock or sell their cattle because of the slump in prices due to the prevalence of Foot-and-Mouth cattle disease during the drought.In other dryland instances, such as in the Sahel, bareland areas with minimal vegetation growth alternate with grasslands in response to rainfall variability and drought [50].
Drought impacts on grasslands, forests and wetlands imply negative impacts on the cattle system and biodiversity, including wildlife in savannas with the associated tourism and hospitality sector, whereas effects on croplands impact food production.For example, the water crisis of 2015-2016 resulting from relatively low, erratic rainfalls reduced water levels and water inflows into dams drastically across the country [51].Regarding drought impacts on agriculture, for example, maize production in Botswana declined from 35,322 tons in 2011 to 13,911 tons in 2017 mainly due to drought constraints.This necessitated an expenditure of about P389 Million (~36 Million American Dollars) on maize import to meet cereal requirements not met through domestic crop production [31].Other changes in LULC, such as the changes in the extent of settlements, as observed between 2000 and 2018 in this study, are not drought-related.Settlement expansion is more a reflection of the increasing land demand for human habitation and infrastructure due to the population growth experienced in the CDB.

Drought Severity in the Constituencies
Relating drought severity to the years declared as drought-stricken by the government, the results reveal that not all constituencies were equally affected by drought, as severity differed from severe to mild drought.Moreover, drought severities in some declared drought years were not as widespread in the CDB as in other parts of the country.For example, the growing season of the year 2009-2010 had improved vegetation productivity in response to above-normal rainfall recorded in previous months, which resulted in flooding events in five sub-districts in the CDB (Serowe/Palapye, Tutume, Boteti, Mahalapye and Bobirwa) [52].In other drought-prone contexts in southern Africa, such as Zimbabwe, researchers [53] found that a drought's distribution and effects differed geographically and from season to season.
Drought severity was further gauged by the frequency at the constituency level.Over the study period, the constituencies experienced between eight to two drought events.Examples are Mahalapye west and east, with eight and seven drought occurrences, respectively, ranging from moderate to mild drought.Boteti west experienced seven drought events with severity ranging from severe to mild.Confirming these drought frequencies in the CDB, the authors of [30] found an average drought frequency of two to four (depending on the index) in the Bobonong region between 2000 and 2015.This is in line with our finding of five drought occurrences in Bobonong as our time series extending up to 2018 captured more drought events.Drought effects on land-based resources and livelihoods will vary depending on the drought severity, the land use as well as the land management practices.For example, water use strategy adopted as part of land management practices was found to have impacted the responses of tree plantations to drought in China [54].Droughts cannot be avoided, since these are an integral part of the climate cycle; however, the impacts on social-ecological systems are better minimised with monitoring as input for the use of both proactive and reactive measures [4,21].

Conclusions
This study proposed a spatial and temporal analysis of drought evolution in the Central District in eastern Botswana from 2000 to 2018.The results highlight the usefulness of incorporating land use-land cover and change in assessing the spatio-temporal variability of drought severity in drylands.Remote Sensing-based vegetation time series metrics were used, as complementary to climatological indices.Indicators characterising changes in vegetation conditions and drought severity during the growing seasons (October to March) from 2000-2001 to 2017-2018 were used.These indicators are NDVI difference-NDVI diff , NDVI anomaly, NDVI trends, Vegetation Condition Index-VCI, Drought intensity and frequency.The use of these different NDVI-based indicators, which might seem redundant, are useful as complementary measures since they differ computation-wise.The NDVI difference as utilised in this study captured the in-season variability in vegetation productivity, whereas the VCI compared each growing season with the long-term minimum and maximum conditions.For example, limited vegetation productivity found during the growing seasons of 2002-2004 and 2015-2016, which was based on both NDVI difference and NDVI anomalies, agreed with the heightened levels of drought severity over the same periods as derived from the VCI.
Results further showed high temporal and spatial variability in vegetation productivity between drought and non-drought conditions in our case studies.The associated negative impact of droughts on vegetation resulting in limited vegetation productivity was further confirmed by results from this study.Drought effects on vegetation productivity during the study period were characterised by decreasing vegetation trends in most parts of the district.Although varying intensities of drought severity (severe, moderate and mild) occurred in the constituencies, the 2002-2003 and 2003-2004 growing seasons were found to be the worst drought periods in the entire series, as most parts of the district were affected.Assessing drought severity and intensities by LULC in selected drought years revealed varying drought effects.We found that drought effects differed between LULC types as well as whether these were areas of land change or persistence.Further examination of drought impacts in areas of no change is required, as our understanding of drought effects in areas with no change is still limited.More empirical studies in this regard will provide useful insights.Using the example of the 2002-2003 drought-stricken growing season, the highest percentage of land impacted by extreme and severe droughts were found in tree-covered areas, croplands and grasslands, whereas improved vegetation trends were found mostly in wetlands and some instances in otherland areas including barelands.Moreover, the results suggest that even in declared drought years, droughts severity varied, and the effects differed between constituencies.A further insight provided is that the magnitude of drought severity in some declared drought years was not as widespread in the CDB.For example, no other severe drought levels were recorded in the CDB after the extended drought which affected the growing seasons of 2002-2004.Differences in spatial resolution of the datasets utilised and the coarse spatial resolution of the 1 km NDVI datasets compared to 300 m annual LULC are limitations identified in this study.With the increasing availability of images of higher spatial resolution, such as from SENTINEL-2, results from RS-based analysis of drought can be improved.However, methodological challenges ensue with the need to incorporate the newer images into the existing NDVI archives.For example, consideration ought to be given on parameterising across sensors and balancing the trade-offs between taking advantage of the superior spatial resolution of the newer satellite missions (e.g., Sentinel-1, 2 and 3) and the temporal resolution of images from the older missions (e.g., SPOT VGT and PROBA-V).As the 1 km, NDVI time series datasets extend way back to the 1990s, their use is indispensable for drought analysis at the current time.Moreover, as drought years were easily detected in the time series analyses, this is proof of the usefulness of the 1 km NDVI time series.For example, the lower than normal vegetation productivity during the prolonged drought periods that negatively impacted the 2002-2004 and 2015-2016 growing seasons coincided with strong El Niño years.With the above-normal vegetation productivity in the ongoing season (2020-2021), results suggest the reversal of the negative vegetation trends observed in the preceding growing seasons.How much these negative trends have been reversed remain uncertain, as the season is still ongoing.For clarity, future studies should examine the usefulness of RS-based indices for understanding the ongoing season's phenology in dryland contexts such as the CDB.
Remote Sensing-based time series enabled us to extend the analysis up to the ongoing season, demonstrating its usefulness for better characterisation of drought events.Remote Sensing-based results such as those obtained in this study, when provided at multiple administrative scales in a timely and cost-effective manner, have the potential to aid decision-makers to better plan and respond to drought situations.Scientific evidence is needed as input into the decision-making process to aid national resource mobilisation for drought management.Botswana requires both proactive and reactive approaches for drought management, for which remote sensing-based assessment and monitoring foster the implementation of drought early warning systems.

Supplementary Materials:
The following are available online at https://www.mdpi.com/2072-4292/13/5/836/s1, Table S1: Geographic and socioeconomic characteristics of the constituencies in the Central District, Table S2: Land cover classification scheme, Table S3: Land use-land cover transition matrix in km 2  Funding: This research received no external funding.

Figure 1 .
Figure 1.Study location: (a) Central District in eastern Botswana; (b) Land use-land cover for 2018; (c) Peak distribution of rainfall in the 17 constituencies examined in the study.

Figure 1 .
Figure 1.Study location: (a) Central District in eastern Botswana; (b) Land use-land cover for 2018; (c) Peak distribution of rainfall in the 17 constituencies examined in the study.
and 2015-2018.In other years, such as 2002-2003 and 2010-2011, vegetation productivity was low mostly in the northern and the eastern parts of the district.Improvement in vegetation performance occurred during 2005-2006, 2009-2010, 2010-2011 and 2014-2015.Figure 3a,b compared maximum, mean and cumulative NDVI in the entire time series (2000-2018) with those of the drought years (2002-2003 and 2003-2004) and non-drought year (2009-2012).Vegetation productivity for 2020-2021 as the ongoing growing season is also depicted.Vegetation productivity during the 2009-2010 non-drought growing season was higher than the mean of the entire time series for most months, except in mid-February to March, whereas, for this current growing season (2020-2021), vegetation productivity is well above the mean of the entire time series from mid-October 2020 to February 2021.NDVI anomalies (Figure 3c) captured lower or higher than normal vegetation productivity in the time series.There were some extended periods of very low vegetation productivity over multiple growing seasons.Examples are the period from the middle of the 2002-2004 growing seasons and 2004-2005 to the start of 2005-2006 growing seasons.The improved vegetation productivity during the growing seasons of 2005-2006, 2009-2010, 2010-2011 and 2014-2015 was further confirmed by the NDVI anomaly (Figure3c).NDVI anomalies (Figure3c) captured lower or higher than normal vegetation productivity in the time series.

Figure 2 .
Figure 2. Variability in vegetation productivity based on the seasonal NDVI difference for the growing seasons (October to March) in the time series.

Figure
Figure 3a,b compared maximum, mean and cumulative NDVI in the entire time series (2000-2018) with those of the drought years (2002-2003 and 2003-2004) and nondrought year (2009-2012).Vegetation productivity for 2020-2021 as the ongoing growing season is also depicted.Vegetation productivity during the 2009-2010 non-drought growing season was higher than the mean of the entire time series for most months, except in mid-February to March, whereas, for this current growing season (2020-2021), vegetation productivity is well above the mean of the entire time series from mid-October 2020 to February 2021.NDVI anomalies (Figure 3c) captured lower or higher than normal vegetation productivity in the time series.There were some extended periods of very low vegetation productivity over multiple growing seasons.Examples are the period from the middle of the 2002-2004 growing seasons and 2004-2005 to the start of 2005-2006 growing seasons.The improved vegetation productivity during the growing seasons of 2005-2006, 2009-2010, 2010-2011 and 2014-2015 was further confirmed by the NDVI anomaly (Figure3c).NDVI anomalies (Figure3c) captured lower or higher than normal vegetation productivity in the time series.

Figure 2 .
Figure 2. Variability in vegetation productivity based on the seasonal NDVI difference for the growing seasons (October to March) in the time series.

Figure 3 .
Figure 3. Profiles of NDVI metrics for the growing seasons: (a) Mean and maximum NDVI in the time series compared to multi-drought seasons (2002-2003 and 2003-2004), non-drought season (2009-2010) and 2020-2021 as the current growing season; (b) Seasonal cumulative NDVI for these metrics as in Figure 3a; (c) NDVI anomaly.Grey circles indicate when the anomaly occurred and sizes interpreted as follows: small circle = a part of the growing season was affected (indicated either as start or mid), medium circle = entire growing season was affected, large circle = multi-year growing seasons were affected.

Figure 3 .
Figure 3. Profiles of NDVI metrics for the growing seasons: (a) Mean and maximum NDVI in the time series compared to multi-drought seasons (2002-2003 and 2003-2004), non-drought season (2009-2010) and 2020-2021 as the current growing season; (b) Seasonal cumulative NDVI for these metrics as in Figure 3a; (c) NDVI anomaly.Grey circles indicate when the anomaly occurred and sizes interpreted as follows: small circle = a part of the growing season was affected (indicated either as start or mid), medium circle = entire growing season was affected, large circle = multi-year growing seasons were affected.

Figure 4 .
Figure 4. Spatial patterns of vegetation productivity change between 2000-2001 and 2020-2021 (Makgadikgadi salt pans in the north is indicated as white using the waterbody mask): (a) Vegetation trends; (b) Significance of trends.

4 . 20 Figure 5 .
Figure 5. Evolution of VCI drought severity and intensity through the time series in the Central District: (a) Growing season VCI (October to March); (b) Percentage of land area affected by varying drought severity and non-drought conditions based on the VCI by years.

Figure 5 .
Figure 5. Evolution of VCI drought severity and intensity through the time series in the Central District: (a) Growing season VCI (October to March); (b) Percentage of land area affected by varying drought severity and non-drought conditions based on the VCI by years.

Figure 6 .
Figure 6.Changing land use-land cover conditions between 2000 and 2018: (a) Areas of land use-land cover loss and persistence; (b) Land use-land cover transitions from the year 2000 to 2018 in percentages (areas of persistence for each land use-land cover class are in bold, that is, areas of no change that remained in the same land class over the 18 years).

Figure 6 .
Figure 6.Changing land use-land cover conditions between 2000 and 2018: (a) Areas of land use-land cover loss and persistence; (b) Land use-land cover transitions from the year 2000 to 2018 in percentages (areas of persistence for each land use-land cover class are in bold, that is, areas of no change that remained in the same land class over the 18 years).

Figure 7 .
Figure 7. Distribution of vegetation trends for major land use-land cover categories for the year 2018: (a) Tree−covered area; (b) Cropland; (c) Otherland; (d) Grassland; (e) Vegetation change as percent area of land in (i) 2018 land categories, (ii) areas of land change between 2000 and 2018 and (iii) areas of land persistence where no change occurred between 2000 and 2018.

Figure 7 .
Figure 7. Distribution of vegetation trends for major land use-land cover categories for the year 2018: (a) Tree−covered area; (b) Cropland; (c) Otherland; (d) Grassland; (e) Vegetation change as percent area of land in (i) 2018 land categories, (ii) areas of land change between 2000 and 2018 and (iii) areas of land persistence where no change occurred between 2000 and 2018.

20 Figure 8 .
Figure 8. Distribution of drought severity for major land use-land cover types: (a) Tree-covered area; (b) Cropland; (c) Otherland; (d) Grassland, during the growing season of 2002-2003 (a drought year); (e) Percentage of land area under varying drought severity and non-drought conditions based on VCI for selected drought years by land use-land cover types (T = Tree-covered areas, G = Grassland, C = Cropland, W = Wetlands/Waterbodies, S = Settlement, O = Otherland).

Figure 8 .
Figure 8. Distribution of drought severity for major land use-land cover types: (a) Tree-covered area; (b) Cropland; (c) Otherland; (d) Grassland, during the growing season of 2002-2003 (a drought year); (e) Percentage of land area under varying drought severity and non-drought conditions based on VCI for selected drought years by land use-land cover types (T = Tree-covered areas, G = Grassland, C = Cropland, W = Wetlands/Waterbodies, S = Settlement, O = Otherland).

Figure 9 .
Figure 9. Heatmap of drought severity at constituency level during the growing seasons of 2000-2001 to 2017-2018.The years with dashed lines across were declared drought years by the Botswana government [31].

Figure 9
Figure 9 depicts drought intensities for each year's growing season and drought frequency per constituency in the entire time series.This heatmap is based on the median VCI value for each constituency (heatmaps of drought severity using the minimum and mean VCI values are shown in Figures S1 and S2, respectively).The prolonged drought during the multiple seasons of 2002-2004 is evident in the heatmap.Based on the count of drought occurrences, irrespective of severity classes between 2000-2018, constituencies with the most frequent droughts in descending order are Mahalapye west (eight), Mahalapye east (seven), Boteti west (seven), Shoshong (six), Bobonong (five), Boteti east (five) and Palapye (five).Other constituencies experienced between two and four drought occurrences with lesser severity.Of the 16 declared drought years by the government (the years with dashed lines in Figure 9), eight were evident in the time series, whereas the other years were favourable for the CDB.

Figure 9 .
Figure 9. Heatmap of drought severity at constituency level during the growing seasons of 2000-2001 to 2017-2018.The years with dashed lines across were declared drought years by the Botswana government [31].

Figure 9
Figure 9 depicts drought intensities for each year's growing season and drought frequency per constituency in the entire time series.This heatmap is based on the median VCI value for each constituency (heatmaps of drought severity using the minimum and mean VCI values are shown in Figures S1 and S2, respectively).The prolonged drought during the multiple seasons of 2002-2004 is evident in the heatmap.Based on the count of drought occurrences, irrespective of severity classes between 2000-2018, constituencies with the most frequent droughts in descending order are Mahalapye west (eight), Mahalapye east (seven), Boteti west (seven), Shoshong (six), Bobonong (five), Boteti east (five) and Palapye (five).Other constituencies experienced between two and four drought occurrences with lesser severity.Of the 16 declared drought years by the government (the years with dashed lines in Figure 9), eight were evident in the time series, whereas the other years were favourable for the CDB.
(2000-2018), Figure S1: VCI minimum value heatmap of drought severity at constituency level during the growing seasons of 2000-2001 to 2017-2018.The years in bold are declared drought years by the Botswana government (Source: [41]), Figure S2: VCI mean value heatmap of drought severity at constituency level during the growing seasons of 2000-2001 to 2017-2018.The years in bold are declared drought years by the Botswana government (Source: [41]).