*2.2. Data*

In order to investigate the variability of vegetation in the HiP in response to climatic conditions as well as the recent intense drought of 2014–2016, we opted to use the monthly averaged MODIS Terra/Aqua 16-day datasets measured for the period from 2002 to 2017 (16 years). With its considerable time resolution (about for images per month) compared to other satellites, MODIS images were the most appropriate for this study because of the size of the geographic area. The MODIS data used here are archived in the GEE as image collection. This data product is generated from a MODIS/MCD43A4 version 6 surface reflectance composite. More details about the MCD43A4 MODIS/Terra and Aqua nadir BRDF-adjusted reflectance daily level 3 global 500 m SIN grid V006 data can be found in a study by Schaaf et al. [36]. The data were extracted and processed using the JavaScript code editor in the GEE platform (https://earthengine.google.com/, Mountain View, CA, USA) (see Appendix A), which provides possibilities of parallel computing and large data processing for even very large study areas. For the purpose of this investigation, our main parameter is the NDVI, but we also considered other vegetation indices such as the Enhanced Vegetation Index (EVI), the Burned Area Index (BAI), and Normalized Difference Infrared Index (NDII). The BAI was also included in order to determine the possible vegetation burning activity, which may have been triggered by drier conditions associated with an intense drought period. NDII has been recently proven to be a robust indicator for monitoring the moisture content in the root-zone from the observed moisture state of vegetation [19,21]. These spectral indices were calculated using the formulas:

$$\text{NDVI} = \frac{NIR - R}{NIR + R} \tag{1}$$

$$\text{EVI} = 2.5 \, \frac{NIR - R}{NIR - 6 \, R - 7.5 \, B + 1} \tag{2}$$

$$\begin{array}{l}\text{BAI} = \frac{1}{\frac{(0.1+R)^2 + (0.06+NIR)}{NDI1}}\\\text{NDII} = \frac{NIR - SINIR1}{NIR + SWIR1}\end{array} \tag{3}$$

where *R*, *NIR*, and *SWIR*1 are spectral bands in the blue (450–500 nm), red (600–700 nm), near-infrared (700–1300 nm), and shortwave infrared (1550–1750 nm) regions.

In this study, we derived the precipitation values averaged for the study area for the period from 2002 to 2017 using the Climate Engine Application (CEA, http://climateengine.org/, Moscow, ID, USA), while soil temperature monthly mean data was derived from the National Aeronautics and Space Administration (NASA, Washington, DC, USA): http://giovanni.gsfc.nasa.gov. Both the soil temperature and precipitation data are an output of the Modern Retrospective Analysis for the Research Application (MERRA-2) model [37]. The MERRA model is an American global reanalysis tool operating from 1979 onwards that is based on the NASA Goddard Earth Observation serving Data Assimilation System version 5 (GEOS-5). The MERRA-2 model data are given at a spatial resolution of 0.67◦ × 0.50◦ at 1-hourly to 6-hourly intervals.

There is always an expected variability of surface water content due to changes in both weather and climatic conditions. Therefore, in a study such as this one, it is essential to always investigate the water lost to the atmosphere through both evaporation and transpiration. This can be an important process as it could explain details about vegetation water stress. Given that the study area is a remote area which does not have evaporation and/transpiration measurements records, we opted to use the Global Land Data Assimilation System (GLDAS) evapotranspiration (ET) data. The GLDAS system was designed to generate optimal fields of land surface and fluxes, and it is also capable of generating quality controlled, spatially and temporally consistent, terrestrial hydrological data including ET and other related parameters [38].

The ENSO phenomenon influences rainfall and temperature conditions largely over southern Africa [39,40]. Previous studies have demonstrated how vegetation responds significantly to ENSO [40] and the DMI [41] index as a measure of climatic conditions [42–44] in some parts of southern Africa. Thus, in order to investigate changes in vegetation in the HiP due to variability in climatic conditions, it is important to consider these climate indices. In this study, we used the Niño3.4 monthly mean time series retrieved from the National Oceanic and Atmospheric Administration (NOAA) website (https://www.esrl.noaa.gov/psd/gcos\_wgsp/Timeseries/, Washington, DC, USA). The Niño3.4 index is calculated by taking the area-averaged sea-surface temperature (SST) within the Niño3.4 region, which is at 5◦ N–5◦ S longitude and 120◦ W–170◦ W latitude in the Pacific Ocean. On the other hand, the DMI is calculated by taking the difference between the SST anomalies in the western (50◦ E–70◦ E; 10◦ S–10◦ N) and eastern (90◦ E–110◦ E), (10◦ S–0◦ N) sectors of the equatorial Indian Ocean [41]. The DMI data were downloaded from the website: http://www.jamstec.go.jp/frcgc/research/d1/ iod/iod/dipole\_mode\_index.html. The relevant time series of Niño3.4 and DMI are shown in Figure 2.

**Figure 2.** The standardized monthly Niño3.4 (**a**) and dipole mode index (DMI) (**b**) time series for the period from 1980 to 2017.
