*2.5. Wavelet Transforms and Wavelet Coherence*

In this study, we opted to employ the wavelet transform analyses method [55] because of its ability to obtain a time–frequency representation of any continuous signal. Basically, the continuous wavelet transform (CWT) of a given geophysical (in this case) time series is given by transforming the time series into a time and frequency space. While there are several types of wavelets, the choice of the wavelet function is determined by the data series, of which, for geophysical data, the Morlet wavelet function has been shown to perform well [55–57]. Thus, the CWT [*Wn* (s)] for a given time series (*xn, n =* 1, 2, 3, . . . , *N*) with respect to wavelet Ψ<sup>0</sup> (η) is defined as:

$$\mathcal{W}\_n^X(\mathbf{s}) = \sum\_{n'=1}^{n-1} X\_{n'} \Psi^\* \left[ \frac{\left(n' - \ \mathbf{n}\right)}{\mathbf{s}} \delta t \right] \tag{11}$$

where *s* is the wavelet scale, *n* is the translated time index, *n* is the localized time index, and Ψ\* is the complex conjugate of the normalized wavelet*. δt* is the uniform time step (which is months in this case). The wavelet power is calculated from *|Wn (n)|2*. In this study, the CWT statistical significance at a 95% confidence level was estimated against a red noise model [55,57]. Using a continuous wavelet transform analysis, it is also possible to quantify the relationship between two independent time series of the same time step *δt*. In this study, the aim was to quantify the relationship between NDVI averaged for the study area and selected climatic parameters. Following Grinsted et al. [57], for the two time series of *X* and *Y*, with different CWT *W<sup>X</sup> <sup>n</sup>* (*s*) and *W<sup>Y</sup> <sup>n</sup>* (*s*) values, the cross-wavelet transform *Wxy <sup>n</sup>* (*s*) is given by

$$\mathcal{W}\_n^{XY}(\mathbf{s}) = \mathcal{W}\_n^X(\mathbf{s})\mathcal{W}\_n^{Y\*}(\mathbf{s}) \tag{12}$$

where "\*" represents the complex conjugate of the *Y* time series. The output of the above equation can also assist in calculating the wavelet coherence. Basically, wavelet coherence is a measure of the intensity of the covariance of the two time series in a time–frequency domain. This is an important parameter because the cross-wavelet only gives a common power. Another important process is to calculate the phase difference between the two time series. Here, the procedure is to estimate the mean and confidence interval of the phase difference. Following a study by Grinsted et al. [57], we used the circular mean of the phase-over regions with relatively high statistical significance that are inside the cone of influence (COI) to quantify the phase relationship between any two independent time series. As defined in a study by Zar [58] and also later by Grinsted et al. [57], the mean circulation of a set of angles (*ai*, *i* = 1, 2, 3, . . . , *n*) can be defined by the following equation:

*Climate* **2018**, *6*, 95

$$a\_{\mathfrak{m}} = \arg(X, Y)\text{ with }X = \sum\_{i=1}^{n} \cos(a\_i) \text{ and }Y = \sum\_{i=1}^{n} \sin(a\_i) \tag{13}$$

Following the studies by Torrence et al. [55,57], the wavelet coherence between two independent time series can be calculated using the following equation:

$$R\_n^2(s) = \frac{\left| S\left(s^{-1} \mathcal{W}\_n^{XY}(s)\right) \right|^2}{S\left(s^{-1} \left| \mathcal{W}\_n^X(s) \right|^2\right) \times S\left(s^{-1} \left| \mathcal{W}\_n^Y(s) \right|^2\right)}\tag{14}$$

where the parameter *S* is the smoothing operator defined by *S(Wn (s)) = S*scale (*S*time *(Wn(s*))). The parameter *Stime* represents the smoothing in time. For further details about the theory of wavelet analyses, the reader is referred to [55,57,59].

#### **3. Results and Discussion**

To investigate whether the El Niño event of 2014–2016 can be classified as a strong El Niño event, a time series for the period from the beginning of the satellite era (1980) to 2017 was plotted (see Figure 2a). We also considered the DMI index (Figure 2b) as a measure of climatic conditions of the eastern part of southern Africa [43]. A general classification of ENSO events should contain 5 consecutive overlapping 3-month periods with SST anomalies below −0.5 for the La Niña events and above +0.5 for the El Niño events.

In Figure 2a,b, both the El Niño events and positive DMI are shaded in red, whereas La Niña and negative DMI are indicated in blue. To identify the strength of the ENSO events, the threshold is further broken down to weak (0.5–0.9 SST anomaly), moderate (1.0–1.4 anomaly), and strong (≥1.5 anomaly) events. Figure 2a shows that the 2014–2016 El Niño was one of the strongest since the beginning of the record. Other notably strong El Niño occurrences were in 1982/1983, 1997/1998, and 2009/2010. On the other hand, there were many episodes of positive DMI, with one such event in the 2014–2016 period, which seems to be in phase with the recent strong El Niño of 2014–2016.

A composite of the NDVI index averaged for each year from 2002 to 2017 is shown in Figure 3. In this figure, regions where there are greener colors indicate higher NDVI values, whereas the brownish colors indicate low NDVI values. These results show that there seems to be a direct influence of the ENSO in the vegetation of the HiP, especially during strong El Niño years (2014–2016). It is evident that during El Niño years, there was a decline in NDVI values especially in the southern and western parts of the study area. This is presumably because the vegetation of the northern part of the HiP is dominated by a forest which is consist of indigenous trees which are believed to be drought resistant (see Figure 1). Additionally, the contributing factor could be that the eastern part of the HiP is benefiting from orographic lifting as it is situated in a high terrain (see Figure 1). The evidence of the influence of El Niño is more prominent during the strong El Niño years such as 2003 and the recent intense 2014–2016 drought period, as well as the 2008 non-ENSO drought period.

Figure 4a shows the deseasonalized monthly averaged MODIS NDVI time series for HiP from 2002 to 2017 (red line) plotted together with the 12 months running mean smooth trend (black dotted line). The monthly mean NDVI values plotted in Figure 4a were calculated by taking an averaged of MODIS images available in that month. In this study area, the MODIS satellite records four images per month. In general, there is a steady trend of NDVI measured at the HiP beside some anomalies observed in specific parts of the time series. This seems to be the case for southern Africa because other studies also indicated a steady trend for this region [10]. Remarkably, during the 2014–2016 period, a period that coincided with the recent intense El Niño, there was a sudden decrease in the NDVI values which reduced to the lowest minimum value of about 0.3 in November 2015. During this period, EVI values also decreased to minimum values of about 0.11 (results not shown here).

**Figure 3.** The spatiotemporal variability of normalized difference vegetation index (NDVI) at the Hluhluwe-iMfolozi Park for the period from 2002 to 2017. The scale represents the range of NDVI values from 0 to 1.

**Figure 4.** (**a**) The deseasonalized monthly mean NDVI time series for HiP. The continuous red line indicates the trend estimate and the dashed red lines show the 95% confidence interval for the trend based on resampling methods. (**b**,**c**) show the histogram and yearly mean time series, respectively.

A study by Mberego and Gwenzi [60] investigated the temporal patterns of precipitation and vegetation variability over Zimbabwe during extreme dry and wet rainfall seasons using data covering the period 1981–2005. Their NDVI time series indicated a steady trend over this period. However, it seemed to be strongly affected by severe dry conditions, an observation which is consistent with the results presented here. In this study, the deseasonalized monthly mean NDVI time series in Figure 4a

(red line) indicates the possible response that corresponds to both dry and wet years, especially during the most recent strong El Niño events of 2003 and 2014–2016. In relation to the strength of the influence of El Niño in the south-western part of southern Africa, a study by Manatsa et al. [61] analyzed agricultural drought in Zimbabwe using the standardized precipitation index (SPI). They reported the 1991–1992 period as the period which experienced the most extreme drought conditions. A little later, observations by Mberego and Gwenzi [60] reported the year 2002–2003 as the drought period with the most prolonged time of relatively low NDVI values in their time series. While our study does indicate a significant influence of the 2003 El Niño event in the HiP NDVI values, the observations presented here indicate that 2014–2016 was the longest period with low NDVI values. Thus, 2014–2016 could be regarded as the most recent intense El Niño period, with a maximum effect on vegetation in the HiP. The NDVI values dropped from a value ~0.65 in November 2013 to 0.3 in November 2015. This is also verified by a much smoother representation of the NDVI in Figure 4c in which a reduction in NDVI values is observed. This reduction coincides with the most recent strong El Niño. Additionally, a reduction which coincides with the El Niño of 2003 (see Figure 4c). Another significant feature is a strong peak, which reaches ~0.8 just after the Irina tropical storm, which occurred in early March 2012 (see Figure 4).

Figure 5 shows monthly mean time series values plotted together with their corresponding 12-month running-mean smooth trend for NDVI (Figure 5a), EVI (Figure 5b), BAI (Figure 5c), soil temperature (Figure 5d), precipitation (Figure 5e), evapotranspiration (Figure 5f), and NDII (Figure 5g). These monthly mean values are plotted together with their respective smooth trends, which were calculated using the Breaks For Additive Season and Trend (BFAST) method, which is described in details by Verbesselt et al. [62,63]. Basically, the BFAST method integrates the decomposition of time series seasonal, trend, and remainder components of any satellite image time series, and can be applied to any other type of time series in the geosciences that deals with seasonal or non-seasonal time series. The period of the most recent intense drought (2014–2016) is indicated by the grey shaded box in each figure. In general, all the parameters show a seasonal cycle in terms of monthly means.

There is an expected resemblance between the NDVI and EVI observations in both the monthly mean time series and the smooth trend, with a clear indication of the effect of the 2014–2016 drought period. These observations are consistent with a study by Xulu et al. [10], who investigated the response of commercial forestry to the recent strong and broad El Niño event in a region which is 70 km south-east of the HiP. In their study, Xulu et al. [10] reported a significant decline of NDVI values which corresponded to the 2014–2016 El Niño years [10]. Although the influence of the 2014–2016 El Niño in the HiP seems to be the strongest, it follows the same pattern as that reported by Anyamba et al. [40] in their study of the influence of both El Niño and La Niña in the vegetation status over eastern and southern Africa. Considering the level of browning of vegetation demonstrated in Figure 3 for the years 2014–2016, it is necessary to consider the possible fire activity given the relatively dry conditions in the HiP. Figure 5c indicates that during the period of the intense drought of 2014–2016 there was an increase in fire incidences in the HiP. This is revealed by a rise in the BAI values of the smooth trend to its maximum level of approximately 50 in November 2015. During the 2014–2016 period, the HiP experienced an unprecedented decline in the total precipitation per month (see Figure 5d). During the same period, the soil temperature increased to its highest maximum (see Figure 5e). The GLDAS monthly mean ET time series shown in Figure 5f indicates a declining trend during the period 2014–2016, which indicated a possible vegetation stress. In order to investigate the moisture content at root-zone, the NDII index was used. The NDII (Figure 5g) indicates a similar pattern to that of the NDVI and EVI time series. It is observed in Figure 5a that the NDII had a steady trend (0.10) during the period 2002–2013 which was followed by a sudden decrease which reached a minimum value of −0.06 in November 2015.

**Figure 5.** The monthly mean time series values of (**a**) NDVI, (**b**) Enhanced Vegetation Index (EVI), (**c**) Burned Area (BAI), and Modern Retrospective Analysis for the Research Application (MERRA-2) model soil temperature (**d**) and precipitation (**e**), Global Land Data Assimilation System (GILDAS) evapotranspiration (**f**) and Normalized Difference Infrared Index (NDII) (**g**). The dotted lines represent the 12-month smooth trends.

The 12-month running mean smooth trends extracted using the BFAST method for NDVI, EVI, and BAI plotted against anomalies of climatic forcers Niño3.4 and DMI are shown in Figure 6. This plot was constructed to investigate any possible 2-dimensional teleconnection between vegetation condition and the Niño3.4 and DMI climatic forcers, respectively. The panels on the left represent the vegetation indices versus Niño3.4, and the panels on the right show the vegetation indices versus DMI. Both the NDVI (Figure 6a) and EVI (Figure 6b) values show a fairly steady pattern for most parts of the time series, which vary between NDVI values of 0.50 and 0.60, and between EVI values of 0.28 and 0.34. However, both the NDVI and EVI values seem to be enhanced by the extreme amount of rainfall that was brought by the tropical cyclone Irina during early 2012 in the eastern part of southern Africa. In that year, NDVI values increased to a maximum value of approximately 0.62, whereas the more sensitive EVI index reached its maximum of approximately 0.38. The strong peaks that were observed during 2004 for both the NDVI and EVI time series correspond to the greening of vegetation in the HiP which was produced by heavy rainfall that was brought by tropical cyclone Elite in January 2004 [64]. NDVI values were observed to decrease sharply from late 2013 until they reached their minimum of approximately 0.40 in 2015 before beginning to recover to normal average conditions in 2017. This pattern is also depicted in the EVI time series and is directly linked to the stronger and more extensive 2014–2016 El Niño event. Similar results were also presented in a study by Xulu et al. [10], who investigated the influence of recent droughts on forest plantations in Zululand. The notable browning observed in Figure 3 for years 2014, 2015 and 2016, which was also revealed by the NDVI and EVI time series (Figure 5), seems to represent favorable conditions for biomass burning in the

HiP. This is revealed by the unprecedented sudden increase of BAI values to its highest maximum of approximately 50, which coincides with the enhancement of Niño3.4 during 2014–2016 (Figure 6e).

**Figure 6.** The (**a**,**b**) NDVI, (**c**,**d**) EVI and (**e**,**f**) BAI (blue dashed line) 12-month smooth trends versus Niño3.4 (left panels) and DMI (right panels) for the period from 2002 to 2017 for the HiP.

The DMI was highly variable compared to the Niño3.4 climatic forcer throughout the study period, with several distinctive positive DMI values that reached a maximum of just below 1.0. Remarkably, there is a strong peak that extends up to approximately 0.8 during the period band corresponding to 2014–2016 that coincided with the decline in NDVI and EVI and the increase of the ENSO and BAI time series. We note here that the widespread browning observed during the 2014–2016 drought period could have been accelerated by the fact that the climatic forcers, which are known to influence the south-eastern part of southern Africa, may have been in phase during this period. This, of course, needs further investigation; and is discussed below.

#### *3.1. Correlations Statistics and Mann–Kendall Test*

The Pearson correlation between NDVI and climatic variables used in this study for the whole study record was derived. Figure 7 shows the heat map which summarizes the linear relationships between all the parameters monitored in this study. In this figure, it is clear that there is a strong correlation between NDVI and Soil temperature (*r* = 0.35), precipitation (*r* = 0.43), ET (*r* = 0.68), and NDII (*r* = 0.92). On the other hand, there is a significance strong negative correlation between

the NDVI and BAI, which is not surprising because greener vegetation reduces chances of biomass burning, while the possibility of the satellite detecting a charcoal signal from burnt vegetation during dry conditions is high. There is also a noteworthy negative (*r* = −0.27) correlation between NDVI and Niño3.4. The results shown in Figure 7 also reaffirm the strong relationship between soil temperature with a strong correlation coefficient of *r* = 0.77. Considering that Figure 2 indicates some episodes where a strong Niño3.4 peak is in phase with the DMI peaks, the noteworthy correlation of *r* = 0.28 between these two climatic indices seems to reaffirm this.

**Figure 7.** The heat map of Pearson correlation coefficients for NDVI, NDII, precipitation (Prec), soil temperature (Soil.Temo), ET, BAI, Niño3.4, and DMI.

Figure 8 shows the inter-annual variability of the Pearson linear correlation between the HiP NDVI values and parameters such as BAI, Soil Temp, Prec, Niño3.4, DMI, ET, and NDII for the period from 2002 to 2017. The correlation between NDVI and EVI was not analyzed because the two parameters closely resemble each other. In general, NDVI is positively correlated to soil temperature, precipitation, ET and NDII through the study period. The NDVI–NDII correlation was the strongest positive correlation with an average value of *r* = 0.91. This reaffirms the strong relationship between vegetation water stress and soil moisture at the root-zone. The NDVI–ET correlation was observed to be steady at an average correlation coefficient value of *r* = 0.65 during the period from 2002–2013. However, this linear relationship decreased to *r* = 0.40 and *r* = 0.30 in 2015 and 2016, respectively. A study by [65] also used MODIS NDVI and GILDAS evapotranspiration data in order to investigate the relationship between NDVI and evapotranspiration. In their study, they reported a steady positive inter-annual variability of correlation coefficients with an average value of *r =* 0.58. As expected, the NDVI–Niño3.4 correlation is dominated by negative values which are observed to decrease during the periods corresponding to El Niño years. This is consistent with previous studies such as those of Xulu et al. [10,40] who reported a significant influence of ENSO on the vegetation of southern Africa, especially the north-eastern part. Moreover, a salient observation is that the greatest minimum correlation recorded was in 2015, a year with a particularly strong El Niño. The negative correlation between DMI and NDVI also seems to be greater during the recent intense drought period, which could indicate that Niño3.4 and DMI were in phase during this time. The correlation between NDVI and the BAI is expected to be strongly negative as greening is not conducive to biomass burning. However, the results presented in Figure 8 indicate that there was a sudden increase in correlation between NDVI and BAI in 2015 before it returned to its average position in 2016 and 2017. Overall, the inter-annual

variation of almost all the study parameters indicates a noticeable change during El Niño events in the years 2003 and more prominently during the 2014–2016 period.

**Figure 8.** The inter-annual variability of linear correlations between NDVI and BAI, Soil Temp, Prec, Niño3.4, DMI, ET, and NDII for the period 2002 to 2017.

A comprehensive summary of the MLR analysis statistics encompassing NDVI, temperature, precipitation, Niño3.4, and DMI is shown in Table 1. It should be mentioned that the soil temperature, precipitation, ET, NDII, and Niño3.4 were used in this model because of their well-known possible influence on NDVI variability. The DMI climatic parameter was not used as an explanatory variable in the MLR model because of its weak correlation with the NDVI. The results in Table 1 reveal a statistically significant relationship between NDVI and soil temperature and between NDII and ET, with *<sup>p</sup>*-values of 0.000386, <2.00 × <sup>10</sup>−<sup>16</sup> and 0.000173, respectively. Both Precipitation and Niño3.4 indicate a statistically insignificant association with the NDVI because of *p*-values which are far greater than 0.05. A positive significant correlation between NDVI and Soil temperature, NDII, and ET, which is also represented as in Figures 7 and 8, indicates that soil moisture, soil temperature, and evapotranspiration play a significant role in vegetation health in the HiP. The significant but negative correlation between Niño3.4 and NDVI confirms the notion that ENSO variability plays a role in the climatic conditions of southern Africa [35,52].

**Table 1.** The output of the Multiple Linear Regression (MLR) model in which Normalized Difference Vegetation Index (NDVI) is a dependent variable and soil temperature, precipitation (Soil Temp), Niño3.4, Normalized Difference Infrared Index (NDII), Dipole Model Index (DMI), and Evapotranspiration (ET) are independent variables.


Significance codes: 0 '\*\*\*' 0.001 '\*\*' 0.01 '\*' 0.05 '.' 0.1 ' ' 1.

In this study, the Mann–Kendall trend test was used for the analysis of the trend in the HiP NDVI time series. The main advantage of this technique is that it provides a non-parametric test that does not require data to be normally distributed, and it is not dependent on the magnitude of data. Furthermore, this non-parametric test method has a low sensitivity to abrupt breaks in heterogeneous

time series [66]. The Mann–Kendall test model was applied to the NDVI data, and the results are shown in Figure 9. In summary, the z-score and *p*-value for the entire NDVI time series period (2002–2017) were found to be −1.22 and 0.224, respectively. Both the *z*-score and the *p*-value seem to indicate that there was a downward but not significant trend in the NDVI data. The indication of an insignificant downward trend (negative *z*-score) presumably due to the unprecedented sudden reduction of the NDVI values which coincided with the 2014–2016 drought. In order to investigate the influence of drought conditions in the study area using the Mann–Kendall method, it is necessary to calculate the inter-annual variation of Mann–Kendall *z*-scores. These Mann–Kendall *z*-scores were calculated from monthly means for each year starting from 2002–2017.

**Figure 9.** (**a**) The inter-annual variation of Mann–Kendall *z*-scores (α = 0.05, *Z*1 = –1.96, *Z*2 = 1.96) for the HiP from 2002 to 2017. (**b**) Sequential statistics values of progressive (Prog) *u*(*t*) (solid red line) and retrograde *u* (*t*) (black solid line) obtained by Sequential Mann-Kendall (SQ-MK) test for HiP monthly mean NDVI data for the period from 2002 to 2017.

Figure 9a shows the Mann–Kendall *z*-scores based on the 16 years of monthly average NDVI data for the game reserve. In general, it is expected that vegetation will respond to climate fluctuating conditions, and this is clearly depicted by significantly negative *z*-scores (less than *Z*1 = −1.96) during strong El Niño events (e.g., in 2003 and 2014/2015). The significant downward trend observed between 2014 to 2015 is the strongest such downward trend in the history of the MODIS NDVI data used in this study; it demonstrates a clear response of the vegetation of the reserve to the strong El Niño event of 2014–2016. Similar analysis and results comparable with those presented here were reported by Hou et al. [24] in their study on the inter-annual variability in growing-season NDVI and its correlation with climate variables in the south-western Karst region of China.

The sequential version of the Mann–Kendall test was applied to the NDVI monthly mean time series so as to determine the approximates time periods of the beginning of a significant trend. Figure 9b shows the sequential statistic values of forward/progressive (Prog) *u*(*t*) (solid red line) and retrograde (Retr) *u* (*t*) (black solid line) obtained by SQ-MK test for HiP monthly mean NDVI data for the period from 2002 to 2017. There is a noticeable statistically significant downward trend which seems to coincide with the 2003 and strongly the 2014–2016 strong El Niño event. These independent calculations are in agreement with the inter-annual variation of the Mann–Kendall *z*-scores results presented in Figure 9a. In the case that seems to be associated with the 2014–2016 strong El Niño events, there is an apparent downward trend (indicated by the retrograde) that begins in November 2012 and reaches the negative significant trend limit (−1.96) in April 2014. The retrograde statistic values stay in significant negative territories during the period from April 2014 to May 2016 before it starts to revert back to be within the 95% confidence level limits (±1.96). This trend is regarded as significant because the progressive and retrograde curves intersect each other (turning point) within the limits of the 95% confidence level. This significant trend turning point took place during November 2012. Another significant downward trend was observed in late 2003, with the significant trend turning point observed in June 2005.

The intensity of the 2014–2016 drought impact in the HiP has been identified to be identical to that of the early 1980s [11]. Some of the additional factors that reportedly intensified the impact of the 2014–2016 drought include the reduction in the grazing lawns, siltation of rivers, and the increasing number of carnivores [11]. The impact of the 2014–2016 drought did not only affect this natural protected area (HiP), but also the comical plantations which are situated at about 70 km southwest of the HiP [10], [67]. A study by Crous et al. [67] reported a large-scale dieback of *Eucalyptus grandis* × *E. urophylla* (SClone) in the Zululand coastal plain, KwaZulu-Natal, South Africa, during the recent intense drought. This was later supported by Xulu at al. [10], where they reported that the commercial forest of kwaMbonambi, northern Zululand suffered drought stress during 2015.
