Spatio–temporal Assessment of Drought in Ethiopia and the Impact of Recent Intense Droughts

Abstract

The recent droughts that have occurred in different parts of Ethiopia are generally linked to fluctuations in atmospheric and ocean circulations. Understanding these large-scale phenomena that play a crucial role in vegetation productivity in Ethiopia is important. In view of this, several techniques and datasets were analyzed to study the spatio–temporal variability of vegetation in response to a changing climate. In this study, 18 years (2001–2018) of Moderate Resolution Imaging Spectroscopy (MODIS) Terra/Aqua, normalized difference vegetation index (NDVI), land surface temperature (LST), Climate Hazards Group Infrared Precipitation with Stations (CHIRPS) daily precipitation, and the Famine Early Warning Systems Network (FEWS NET) Land Data Assimilation System (FLDAS) soil moisture datasets were processed. Pixel-based Mann–Kendall trend analysis and the Vegetation Condition Index (VCI) were used to assess the drought patterns during the cropping season. Results indicate that the central highlands and northwestern part of Ethiopia, which have land cover dominated by cropland, had experienced decreasing precipitation and NDVI trends. About 52.8% of the pixels showed a decreasing precipitation trend, of which the significant decreasing trends focused on the central and low land areas. Also, 41.67% of the pixels showed a decreasing NDVI trend, especially in major parts of the northwestern region of Ethiopia. Based on the trend test and VCI analysis, significant countrywide droughts occurred during the El Niño 2009 and 2015 years. Furthermore, the Pearson correlation coefficient analysis assures that the low NDVI was mainly attributed to the low precipitation and water availability in the soils. This study provides valuable information in identifying the locations with the potential concern of drought and planning for immediate action of relief measures. Furthermore, this paper presents the results of the first attempt to apply a recently developed index, the Normalized Difference Latent Heat Index (NDLI), to monitor drought conditions. The results show that the NDLI has a high correlation with NDVI (r = 0.96), precipitation (r = 0.81), soil moisture (r = 0.73), and LST (r = −0.67). NDLI successfully captures the historical droughts and shows a notable correlation with the climatic variables. The analysis shows that using the radiances of green, red, and short wave infrared (SWIR), a simplified crop monitoring model with satisfactory accuracy and easiness can be developed.

1. Introduction

In the era of climate change, there is a continuous need to thoroughly assess vulnerabilities caused by complex environmental, ecological, and anthropogenic factors. Drought, as a natural phenomenon, creates numerous multidimensional effects on agriculture, human health, and disease prevalence [1]. Various drought management and vulnerability schemes were thus developed to mitigate the influences of natural and human-made disturbances at regional [2,3] and global scales [4,5]. Vulnerability assessment of natural disasters has become a necessity for policy-makers and practitioners in reducing the impacts associated with them [6,7].

Drought is dryness due to an acute shortage of water, which lasts for several months or years. Drought considerably endangers food and water security. As a complex natural event, it stems from a lack of precipitation over a prolonged period of time, and its effect can be only witnessed slowly over a period of time [8,9]. Besides the shortages of precipitation, droughts are associated with differences between actual and potential evapotranspiration, soil moisture deficits, and reduced groundwater or reservoir levels. These characteristics make the definition of drought complex and, thus, there is no single universally accepted definition. Owing to the lack of comprehensiveness of a single agreed definition, the identification and monitoring of key characteristics of drought is difficult.

Several studies have provided comprehensive reports on indices that are used to monitor the impacts of droughts [10,11,12,13,14,15]. Generally, a variety of drought indices were developed from climatic and satellite data. The most widely used indexes include the Palmer Drought Severity Index (PDSI), Standardized Precipitation Index (SPI) [16], normalized difference vegetation index (NDVI) [17], Normalized Difference Water Index (NDWI) [18], Vegetation Condition Index (VCI), and Temperature Condition Index (TCI) [13]. Remote sensing data-based indices have been widely used and compared with the other approaches for assessing drought, as they are among the best in detecting the onset of drought and measuring the intensity, duration, and impact of drought globally [19]. The remote-sensing based indices for quantifying the state of vegetation, namely the combination of visible and infrared bands, provide unique characterization for the vegetative area, including biomass, growth status, and leaf area coverage, and serve as a basis for the estimation of vegetation condition [20]. Surface temperature may serve as a basis for the estimation of vegetation condition and evapotranspiration [21]. The performance of drought indices generated based on MODIS reflectance and land surface temperature (LST), in association with the standardized precipitation index (SPI), were extensively investigated to assess drought conditions on a global scale to regional scale in the southern Great Plains, USA [22], China [23], in eastern Africa, and in southern and southeastern Africa [24,25,26,27,28].

Ethiopia faces drought conditions every eight–ten years [29]. The country has been facing drought at a growing incidence throughout the past many decades [30]. Among these, the 1984–1985 drought affected the lives of more than two hundred thousand people and millions of livestock [31]. The climate in Ethiopia is changing, even though significant trends are not clear [32]. OXFAM reports that according to the survey made questioning local people in Ethiopia, the climate is experiencing an increase in the rate of drought [33]. The farmers report that good harvests are less common due to an extended extreme dry season and strong rain in the wet season, followed by a prolonged absence of precipitation, which is likely due to a manifestation of global warming. Both the rise in temperature and the long absence of precipitation are major factors for causing droughts. The projected increase of weather events such as droughts due to climate change derails the availability of water and will lead to a cut in agricultural production.

Ethiopia’s economy is essentially dependent on rain-fed agriculture, which is vulnerable to climate change [34]. 2015 was one of the driest years in large parts of Ethiopia [35]. The main rain season, locally called ‘kiremt’, was late and below normal conditions [36]. Consequently, the government called for emergency assistance for 10.2 million people [37]. The ultimate causes of this drought event originated from great distances, through atmospheric and oceanic circulations. The El Niño–Southern Oscillation (ENSO) phenomenon hugely impacts Ethiopian rainfall [38]. In particular, the warm-phase El Niño is closely linked with reserved rains during kiremt, over northern and central Ethiopia [39]. Under these circumstances, the evapotranspiration needs of plants were not met, leading to an intense reduction in vegetative production. Thus, the need to assess long-term vegetation trends and investigate the relationship between these changes and the variability in climatic conditions is increasingly important in Ethiopia.

The specific objectives of this research are: (i) to detect any long-term hydro-meteorological trends using the Mann–Kendall statistical test; (ii) to assess the drought patterns using the vegetation condition index; and (iii) to identify the main causes of NDVI change in relation to rainfall, soil moisture, LST, and ENSO. Additionally, this paper will be the first to attempt to incorporate the Normalized Difference Latent Heat Index (NDLI) as a proxy to evapotranspiration needs of the plant. NDLI, a combination of the green, red, and SWIR channels of the electromagnetic spectrum, has been found to be useful for the detection of plant water content [40]. It is highlighted that a better analysis of drought allows for the development and implementation of successful policies to better understand disruptive climate change in the region, to improve food security and strengthen climate resilience.

2. Study Area and Data

2.1. Study Area

The study area, Ethiopia, is located between 3°00′ to 15°00′N and 32°00′ to 48°00′E in the inner part of the Horn of Africa, as shown in Figure 1. The country has a total area of 1.1 million square kilometers, is landlocked, and has the second largest population in Africa, second to Nigeria. The elevation ranges from 194 to 4539 m above mean sea level. The highland, with an altitude of 1500 m or above, is located at the central and northern parts of the country and constitutes roughly 35% of the country [41]. In a traditional way, based on elevation, at least three climatic zones are identified—the tropical (lowland zone), which is below 1830 m in elevation and has mean annual temperatures of 20–28 °C; the subtropical zone, which includes the highland areas of 1830–2440 m in elevation and with mean annual temperatures of 16–20 °C; and the cool zone, which is above 2440 m in elevation and with mean annual temperatures of 6–16 °C [42].

Due to its complex topographical and geographical features, the climate of Ethiopia exhibits strong spatial variation and different rainfall regimes [43]. Thus, rainfall shows considerable spatial heterogeneity in Ethiopia [44]. Much of the region is generally bimodal, with long rains in JJAS (June–September) and short rains during OND (October–December). The meridional translation of the Intertropical Convergence Zone (ITCZ) across the equator is the main factor of the MAM (March–May) and OND seasons [45]. Topography also plays a role in affecting the annual cycle of precipitation. The highland areas receive an annual rainfall of about 1200 mm, with the least temperature variation, whereas the lowland areas (Afar and Somali regions) receive an annual rainfall of less than 500 mm with larger temperature variations [41]. The spatial distribution of Ethiopian drought indicates that most of the drought and food crises events are concentrated in the central and northern highlands, extending from North Shewa through Wollo to Tigray [46].

Based on the Climate Hazards Group Infrared Precipitation with Stations (CHIRPS) daily precipitation data obtained from the Climate Hazards Group at the University of California, Santa Barbara (UCSB) [35], the main rainfall season from June to September, locally called kiremt, accounts for 60–80% of the annual rainfall, with the remaining falling in the dry season, from October to May, Figure 2.

The rainfall significantly varies between the northeastern and the western highlands of Ethiopia, where orographic rainfall is substantial. Figure 2 additionally depicts that an average seasonal LST of the land derived from the solar radiation (MOD11A2 Terra v.006 product) of Ethiopia is between 10 °C and 54 °C, with maximum temperatures concentrated on the lowland areas. Similarly, the soil moisture (derived from the FLDAS Noah Land Surface Model L4) distribution for the top 0–10 cm layer increases in the western and northern parts of the country. Moreover, NDVI distributions derived from the MOD13Q1 Terra v.006 product confirm healthy vegetation and forests, mainly located in the western parts of Ethiopia, which match with the rainfall, LST, and soil moisture patterns.

The land cover types for Ethiopia extracted from the European Space Agency (2016) Global Land Cover map are shown in Figure 3. The most dominant cover types are grassland, cropland, and shrubs, covering 29%, 26%, and 21% of the whole study area, respectively. The land cover in the highlands continually changes because of the persistent agricultural activities and higher population density as compared to the lowlands [47]. Large areas of agricultural farms, where people largely depend on the rain-fed farms, are of major concern due to recurrent drought incidents.

2.2. Datasets

The data for this study were acquired from four sources. With extensively high temporal and spatial resolution as compared to the other satellites, the products of MODIS onboard NASA Terra and Aqua satellites were suited for this study because of their large geographic coverage. We used the monthly averaged MODIS Terra 16 day datasets for the period from 2001 to 2018 (18 years) that are archived in the Google Earth Engine (GEE) image collection. Time series NDVI and LST covering the whole study area at 250 m and 1 km spatial resolution were generated from MODIS/006/MOD13Q1 and MODIS/006/MOD11A1 version 6 surface reflectance composite, respectively. Similarly, surface reflectance products of MODIS/006/MOD09GA were generated for computing the NDWI and NDLI [40]. The data were extracted and processed using the JavaScript code editor in the GEE platform (https://earthengine.google.com/ Mountain View, CA, USA), which offers a parallel computing environment for processing large datasets. For monitoring the spatial and temporal conditions of drought, we chose NDVI, but we also included the other parameters that trigger dry conditions. NDVI, the most common index for remote sensing of vegetation, is known to be saturated over areas with high leaf area indexes. Numerous vegetation indexes using the same set of near-infrared and red channels have been developed, even though these indices do not enjoy the same popularity as NDVI, which is known for its capability to distinguish vegetation from other types of land cover, but is not really designed to sense the water content in the vegetation canopy. Nevertheless, remote sensing of the water content has important implications in agriculture and forestry. For the detection of plant water content, the near-infrared region (NIR) and shortwave infrared regions (SWIR) have been found to be useful. Thus, NDWI is defined in a similar way to NDVI but uses the near-infrared channel to monitor the water content of the vegetation canopy. Fluctuations in the vegetation canopy are indicators of drought stress [18]. Besides NIR and SWIR, two spectral regions of the electromagnetic spectrum have been found to be useful for the detection of plant water content: the green and red channels. Liou et al. [40] recently developed the NDLI, which uses the green, red, and SWIR channels. NDLI is sensitive to water availability for different land covers at the land–air interface and outperforms the different versions of NDWI indexes. The spectral indices are calculated using the following formulas:

$$NDVI= \frac{NIR-R}{NIR+R}$$

$$NDWI= \frac{NIR-SWIR}{NIR+SWIR}$$

$$NDLI= \frac{G-R}{G+ R+SWIR}$$

where $G, R, NIR, SWIR$ are the spectral reflectance for MODIS band 4 (545–565 nm), band 1 (620–670 nm), band 2 (841–876 nm), and band 6 (1628–1652 nm), respectively.

The other data source used to generate time series rainfall data for the period from 2001 to 2018 was CHIRPS. CHIRPS is a 30+ year quasi-global rainfall dataset combining satellite observations from the Climate Prediction Center (CPC) and the National Climate Forecast System version 2 (CFSv2) and in situ precipitation observations [35,37]. It is widely used in Ethiopia for drought monitoring [38]. It is well demonstrated that CHIRPS can complement the sparse rain gauge network and provide high spatial and temporal resolution for trend analysis [48,49]. The CHIRPS data were accessed from the ftp server (ftp://ftp.chg.ucsb.edu/pub/org/chg/products/CHIRPS-2.0).

The monthly soil moisture (0–10 cm) was generated from the Famine Early Warning Systems Network (FEWS NET) Land Data Assimilation System (FLDAS) dataset, developed to assist food security assessments in data-sparse developing countries [50]. This is a natural tool to monitor drought conditions and was accessed from https://earlywarning.usgs.gov/fews/product/308.

In this study, we used the multivariate ENSO index (MEI) and the dipole mode index (DMI) to observe how vegetation responds to climatic conditions. The monthly mean MEI time series were retrieved from the National Oceanic and Atmospheric Administration (NOAA) website (https://www.esrl.noaa.gov/psd/enso/mei/, Washington, DC, USA). The MEI time series was calculated by taking the leading principal component time series of the empirical orthogonal function of the five variables, namely the sea level pressure, sea surface temperature, surface zonal winds, surface meridional winds, and Outgoing Longwave Radiation within the 30°S–30°N and 100°E–70°W region. Besides this, the DMI was calculated by taking the differences between the sea surface temperature anomalies in the western (50°E–70°E, 10°S–10°N) and eastern (90°E–110°E, 10°S–0°N) portions of the Indian Ocean [51]. The monthly DMI data was accessed from the Japan Agency for Marine-Earth Science and Technology (JAMSTEC) website (http://www.jamstec.go.jp/frcgc/research/d1/iod/iod/dipole_mode_index.html). The data used in these study are summarized on

Table 1. Datasets characteristics and source.

Data	Source	Characteristics
Precipitation	CHIRPS	Monthly precipitation at 0.05 × 0.05 from Jan. 2001 to Dec. 2018
NDVI	MODIS	Monthly NDVI at 250 m from Jan. 2001 to Dec. 2018
LST	MODIS	Monthly LST at 1 km from Jan. 2001 to Dec. 2018
NDWI	MODIS	Estimated from surface reflectance at 500 m from Jan. 2001 to Dec. 2018
NDLI	MODIS	Estimated from surface reflectance at 500 m from Jan. 2001 to Dec. 2018
Soil Moisture	FLDAS Noah	Monthly soil moisture (0–10cm) at 0.10 × 0.10 from Jan. 2001 to Dec. 2018
MEI	NOAA	Monthly MEI time series from Jan. 2001 to Dec. 2018
DMI	JAMSTEC	Monthly DMI time series from Jan. 2001 to Dec. 2018

.

3. Methodology

3.1. Identification of Drought

A common way to calculate anomalies is to apply the Standardized Anomaly Index (SAI). SAI is a standardized departure from the long term mean and is calculated as:

$$SA{I}_i=\frac{x_i-\overline{x}}{\sigma }$$

where $x_i$ is the seasonal mean of variable $x$, $\overline{x}$ is the long term seasonal mean and $\sigma$ is the standard deviation of the seasonal mean of all data. The anomaly maps were created by subtracting the seasonal climatology mean from the seasonal values and then dividing this by the standard deviation. The resulting maps depict the intensity of how good or bad the current season is compared with the average situation. Seasonal anomaly maps of precipitation, NDVI, LST, and soil moisture for the years 2015–2018 from the 2001–2014 climatology were computed to identify and quantitatively measure which part of Ethiopia was severely affected in the year 2015 and its recovery to normal conditions.

Besides the common SAI, another method to compare the current NDVI with historical values is the Vegetation Condition Index [52]. The VCI has been extensively used to monitor vegetation conditions [53]. It normalizes NDVI on a pixel-by-pixel basis, scaling between the minimum and maximum values of NDVI:

$$VCI = 100*\left(\frac{NDVI-NDV{I}_{min}}{NDV{I}_{max}-NDV{I}_{min}}\right)$$

where NDVI, $NDV{I}_{min}$, and $NDV{I}_{max}$ are the mean seasonal NDVI, and its absolute long-term minimum and maximum NDVI values, respectively, for each pixel. VCI varies in the range of 0 to 100 percent, reflecting relative changes in the vegetation condition from extremely low to high VCI [52]. As proposed by [52] and recently applied by [54], a threshold value of below 35% is used to indicate drought conditions as shown in

Table 2. Drought categories derived from the Vegetation Condition Index (VCI).

VCI Percentage	Drought Severity Level
>35	No drought
20–35	Moderate drought
10–20	Severe drought
<10	Extreme drought

.

3.2. Mann–Kendall Trend Analysis

The Mann–Kendall method is a non-parametric rank-based test method, which is commonly used to identify a monotonic trend in climate, by using remote sensing and hydro-metrological data [55]. The usefulness of a non-parametric test relies on its resilience to outliers, non-normality, missing values, and seasonality and, hence, it is necessary for this study [56,57,58]. The univariate Mann–Kendall statistic $S$ for time series data $\left(X_k,k=1,2,\dots ,n\right)$ is given as:

$$S= {\displaystyle \sum }_{j< 1}^{n}sgn\left(X_i-X_j\right)$$

where $X_i$ and $X_j$ are the seasonal mean values in years $i$ and $j$, respectively, $i>j$, and $n$ is the length of the time series. The sign of all possible differences $X_i-X_j$ is computed as:

$$sgn\left(X_i-X_j\right)= \{\begin{array}{c}+ 1, if X_i-X_j>0 \\ 0, if X_i-X_j=0\\ -1, if X_i-X_j<0\end{array}$$

When $n\ge 8$, the statistic $S$ is approximately normally distributed with mean E[S] = 0, and variance ${\sigma }^2$ given by the following equation:

$${\sigma }^2=\left\{n\left(n-1\right)\left(2n+5\right)- {\displaystyle \sum }_{j=1}^p{t}_j\left(t_j-1\right)\left(2t_j+5\right)\right\}/18$$

where $t_j$ is the number of data points in the jth tie group, and $p$ is the number of tie group in the time series. The test statistics $z$ is computed as:

$$z= \{\begin{array}{c}\frac{S-1}{\sigma }, if S>1 \\ 0, if S=0\\ \frac{S+1}{\sigma }, if S<1\end{array}$$

Now, $Z$ follows a standard normal distribution whereby its positive (negative) value indicates an upward (downward) trend. If $Z$ is greater than $Z_{\alpha /2}$, where $\alpha$ represents the significance level, the trend is considered as significant. In this regard to the $z$-transformation, this study is considered a 9.5% confidence level, where the null hypothesis was no trend was rejected if $\left|z\right|>1.96,$ and the alternative hypothesis that increasing or decreasing monotonic trend exists in the series was accepted. The magnitude of the linear trend was then predicted by the Sen’s slope estimator [59], i.e., the change per unit time of a trend was computed as:

$$\mathrm{Sen}’\mathrm{s} \mathrm{slope} = \mathrm{Median} \left\{\left(x_i-x_j\right)/\left(i-j\right)\right\},i >j,$$

where $x_i$ and $x_j$ are the changing values of the variable at time steps $i$ and $j$, respectively. A value close to zero means there is not much variation through time. A negative value of the slope depicts a negative trend, whereas a positive value indicates a positive trend. This method is recommended for remote sensing time series analysis and has been used for vegetation trend analysis [60]. The trend analysis described above was applied to the seasonal rainfall, NDVI, LST, and soil moisture values using the “spatialEco” package in R-project.

3.3. Multiple Linear Regression

Multiple linear regression is an extension of simple linear regression. It is used when to predict the value of a variable based on the value of two or more other variables. For instance, for analyzing a dependent variable (in this case NDVI) in light of related independent variables (precipitation, soil moisture, LST, NDWI, NDLI, MEI, DMI). It allows us to determine the overall fit of the model and the relative contribution of each of the predictors to the total variance explained. In this paper, we tried to quantify the susceptibility of NDVI to changes in climatic and hydro-metrological variables. Mathematically, a multiple linear regression model with $k$ predictor variables $x_1, x_2,\dots , x_k$ and a response can be written as:

$$y = {\beta }_0 + {\beta }_1{x}_1+ {\beta }_2{x}_2 + \cdots + {\beta }_k{x}_k + \epsilon \mathrm{where} i=1,2,\dots ,k ,$$

and $\epsilon$ is the residual terms of the model, which tries to minimize, $y$ is the dependent variable in this case NDVI, $x_i$ represents the independent variables (precipitation, soil moisture, LST, NDWI, NDLI, MEI, DMI), ${\beta }_0$ is the intercept, and ${\beta }_1, {\beta }_2, \cdots , {\beta }_k$ are the coefficients of $x_i.$ Before we chose to analyze our data using multiple regression, we made sure that assumptions required for multiple regression were met. We checked the existence of a linear relationship by inspecting the scatter and partial regression plots between NDVI and each of the independent variables. By using the variance inflation factor (VIF) values, we further checked whether the explanatory variables were highly correlated with each other or not. A VIF measures the extent to which multicollinearity has increased the variance of an estimated coefficient. It looks at the extent to which an explanatory variable can be explained by all the other explanatory variables in the equation.

4. Results and Discussion

The most recent ENSO, which was developed in 2014 and strengthened in the summer, has caused global impacts [61]. In Ethiopia, the dry kiremt seasons are closely linked to the significantly warmer Pacific sea surface temperatures [39]. Figure 4 depicts that the strongest kiremt precipitation anomalies derived from the CHIRPS datasets are located in the central and northwestern parts of Ethiopia, with maximum −4.6 standardized deviations anomalies around −460 mm/year in 2015. Vegetation in Ethiopia is sensitive to water availability and severely affected by low precipitation. Correspondingly, large area negative NDVI deviations are a result of water stress concentrated in the western, northern, and central parts of Ethiopia, with maximum NDVI departures by approximately −2.5 standardized anomalies below average. In the same way, the 10 cm soil moisture and LST follow the same patterns as those of the precipitation anomalies, by approximately −3 and 3.5 standardized deviations from their corresponding normal conditions, respectively. During 2016, due to the dry conditions linked with La Niña, the negative precipitations of the southern and eastern parts of Ethiopia persisted, with maximum −4.0 standardized deviations anomalies around −305 mm/year. The dry conditions evolved from the north and central regions to the south and east parts. Across the region, however, NDVI did not follow the same pattern and the vegetation productivity did not quickly decline. This may have been due to the extended availability of water stored in soils for growing crops [62]. Following the return of ENSO to neutral conditions in 2017, the central and northern regions of Ethiopia become more favorable for crop development. During this period, the cropland areas experienced enhanced precipitation and vegetation, which was also closely linked to the increase in soil moisture. The agricultural data obtained from the annual agricultural sample survey of the Central Statistics Agency indicated increments from 7.32 to 28.93 quintals per hectare for maize, from 5.05 to 26.76 quintals per hectare for Teff, and from 2.28 to 29.67 quintals per hectare for wheat [63].

While in 2018 the precipitation showed negative anomalies, the maximum soil moisture and NDVI anomalies were about two standardized deviations above the average conditions. Similarly, the minimum LST departure was about −3 standardized deviations above the average conditions. It is worth mentioning that in 2018, compared to 2017, a higher precipitation in the southeast part of Ethiopia was observed, which was well-matched with increased NDVI.

4.1. Drought Patterns Based on VCI

Figure 5 depicts the spatio–temporal persistence of drought detected by VCI during the growing season in Ethiopia over the past two decades. It is shown that the growing season signifies the maximum vegetation growth, and demonstrates the suitability of VCI to detect drought and assist the measures of vegetation health.

In this figure, regions which are greener indicate vegetation levels higher than the average conditions, whereas the red colors indicate poor conditions. Severe to extreme droughts were identified in the years 2002, 2003, 2004, 2009, 2010, 2012, and 2015 for the north, central, west, and southwest parts of the country, where the land is mainly covered by rain-fed agriculture. The results show a direct influence of ENSO on the vegetation of Ethiopia, especially during the El Niño years 2009–2010 and 2014–2015. During El Niño years, the NDVI values gradually declined and remained marginally below average. On the other hand, the years 2001, 2005, 2006, 2007, 2013, 2016, and 2018 reflect the near-normal NDVI throughout most of the rain-fed agriculture regions. In Ethiopia, an El Niño event would cause suppressed rainfall during the kiremt season, causing serious reductions in cereal yields and output [64]. On the other hand, when a La Niña event followed on from an El Niño, favorable and above average vegetation conditions were observed, for instance 2010–2011 and 2016–2017 La Niña events, which followed on from the 2009–2010 and 2014–2015 El Niño events, respectively.

4.2. Spatial and Temporal Trends

The spatial and temporal variability of the trends, together with the significance of the trends in precipitation, NDVI, soil moisture, and LST, are presented in Figure 6 and Figure 7. The Mann–Kendall test was carried out to observe whether the mentioned variables changed over space during the 18 years period in the country. The areas in green (positive slope value) indicate an increasing monotonic trend in precipitation, NDVI, soil moisture, and LST, whereas areas in red (negative slope value) indicate a decreasing monotonic trend in precipitation, NDVI, soil moisture, and LST.

The pixel-based trend analysis shows the growing season trend values of precipitation range from −26 to 11 mm, with significant changes occurring in the central and northern parts of the country. Specifically, the northern, central, and rift valley regions of Ethiopia experienced a decreasing rainfall trend, whereas western Benshangul and the highlands of the central Amhara region show an increasing trend. On the other hand, the lowland pastoral regions of Somali region did not show a significant trend. Generally speaking, 52.8% of all pixels in the country show a decreasing trend and significant trends concentrate on the central and lowlands regions of the country.

With respect to the NDVI trend, the northern and northwestern areas of the Tigrai and Amhara region, as well as the southern region, showed a decreasing trend during the study period. The growing season NDVI values ranged from −0.0142 to 0.0213, and overall 41.67% of the country indicated a decreasing trend. The significant decreasing trends were located in the northwestern part. Similar pixel-based trend analysis for LST depicted in Figure 7 showed that LST increased for the northwestern, central highland, and southern parts of the country, whereas there was an estimated 11% significance decrease concentrated on the western parts of the Gambella region. These results are in agreement with the recent findings of Workie et al. [65], who used a linear regression approach to detect trends. Similar procedures performed for soil moisture convey that decreasing significant trends can be observed in the central and lowland areas of Afar and Somali regions, whereas the southwestern part of Benshangul and western part of the Gambela region are experiencing a greening trend.

4.3. Multi Linear Regression and Correlation Statistics

To facilitate relationships between NDVI and other parameters, a small box region (38E–39E, 9N–10N) which experienced significant decreasing trends, presented in Figure 4, Figure 6, and Figure 7 was extracted. Figure 8 shows the monthly anomalies time series plots for NDVI and soil moisture (Figure 8a), precipitation and LST (Figure 8b), and NDLI and NDWI (Figure 8c). Basically, the anomalies calculated by subtracting monthly climatology values from each month provide additional information about the variations present. The periods of severe droughts that resulted in countrywide drought conditions during the growing seasons are shaded with a box in Figure 8.

NDVI anomalies in this region were near normal for several years. In contrast, it showed slight green up in 2010 and late 2016, as conditions translate to weak La Niña. Maximum NDVI departures were observed in 2009 and 2015, where NDVI gradually decreased and remained slightly below average. In particular, the 2015 events were accompanied by higher precipitation anomalies of about −100 mm. There is an exact resemblance between the other parameters, with a clear identification of the drought and normal years. Considering the spatial drought patterns derived from VCI (Figure 5), the intense drought years certainly resulted in a decline in soil moisture and water availability. The water stress situations in the root zone were well captured by soil moisture values. The NDWI and NDLI indicate a similar pattern to that of NDVI, where they reached peaks in 2010 and 2016.

The Pearson correlation coefficients between NDVI and other factors (precipitation, soil moisture, LST, NDWI, NDLI, MEI, and DMI) on a seasonal time scale for the whole study record were computed to assess the relationship between them. The Pearson correlation coefficient was conducted using the statistics package in R. Figure 9 shows the heatmap, which summarizes the linear relationships between the parameters. There was a strong correlation between NDVI and precipitation (r = 0.83) soil moisture (r = 0.83), NDLI (r = 0.96), and NDWI (r = 0.63). The positive correlation between precipitation and NDVI implies that an enhanced precipitation supports vegetation growth and vice versa [66]. On the contrary, a significant negative correlation between NDVI and LST (r = −0.76) was observed. Furthermore, a less notable negative correlation of (r = −0.43, r = −0.39) was observed between NDVI and the two climatic indices MEI and DMI, respectively.

Since there are substantial correlations among NDVI, Precipitation, LST, soil moisture, NDLI, NDWI, MEI, and DMI (Figure 9), the detection of multicollinearity is crucial before plugging data into a regression model. Multicollinearity denotes predictors that are correlated with other predictors. The most widely-used diagnostic for multicollinearity is the VIF. We can see from

Table 3. The output of the multiple linear regression (MLR) model, in which NDVI was the dependent variable and precipitation, LST, soil moisture, NDWI, NDLI, MEI, and DMI were independent variables.

Variable	Estimate	Std. Error	t-Value	p-Value	Sig	VIF
Precip.	−2.215 × 10−5	3.788 × 10−5	−0.585	0.5594		2.387
LST	4.577 × 10−4	6.140 × 10−4	0.745	0.456		3.433
SM	5.505 × 10−2	6.104 × 10−2	0.902	0.368		4.622
NDLI	3.697	1.363 × 10−1	27.125	<2.00 × 10−16	***	3.407
NDWI	0.528	1.097 × 10−1	4.813	2.86 × 10−6	***	2.804
MEI	−1.470 × 10−3	7.704 × 10−4	−1.909	0.0576	.	1.089
DMI	−4.552 × 10−4	5.439 × 10−3	−0.187	0.852		1.099

Significance. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1.

that the VIFs are all down to satisfactory values; they are all less than 5. Even though there is some multicollinearity in our data, it is not severe enough to warrant further corrective measures.

The results in Table 3 reveal the statistically significant relationship between NDVI, NDLI, and NDWI and MEI, with p-values of $<2.00 \times {10}^{-16},2.86 \times {10}^{-6}$, and 0.0576, respectively. The significant relationships between NDVI, and NDLI and NDWI make it clear that an increase in water availability causes an upward trend in NDVI, which implies a decline in drought [67]. The results indicate that water availability in the soil was the main influencing factor on the spatially averaged NDVI. The significantly negative correlation between MEI and NDVI reaffirms the claim that ENSO variability plays a major role in the climatic conditions and control vegetation growth conditions of central and northern parts of Ethiopia [68]. The overall multiple linear regression is significant, with a multiple R-squared value of 0.978 and adjusted R-squared value of 0.962. However, precipitation, soil moisture, and LST have insignificant regression coefficients due to their p-values, which are far greater than 0.05. This is due to the interaction (correlation) between the independent variables, and often since p-value is a function of sample size, as well as variance, there is no single rule for setting the “significance” threshold [69]. The insignificant association observed between precipitation and NDVI could also be due to the delayed response of vegetation to precipitation [70], where a time lag effect was not considered in this study. For future prediction, an optimal regression equation (NDVI = −7.01 × 10⁻⁵ + 3.75 × NDLI + 0.518 × NDWI − 1.386 × 10⁻³ × MEI) was obtained via the backward elimination procedure in a stepwise regression analysis, which was achieved by dropping the least significant feature.

5. Conclusions

This study assessed the spatio–temporal variability of drought during the growing season in Ethiopia through VCI, anomaly maps, and trend analysis for the past two decades, from 2001 to 2018. The VCI results identified that severe to extreme countrywide droughts were identified in 2002, 2003, 2004, 2009, 2012, and 2015. On the other hand, the years 2001, 2005, 2006, 2007, 2013, 2016, and 2018 reflected near-normal NDVI throughout most of the rain-fed agriculture regions. These results are coherent with the findings of previous studies in indicating the onset, spatial, and temporal dynamics of agricultural drought in Ethiopia [18]. Pixel-based trend analysis showed that a significant precipitation decrease in the central areas is accompanied by a significant increase in LST. The increase in temperature in the growing season is of major concern, as it implies an increase in evapotranspiration and, thus, affects crop yields. Also, the browning in northwestern parts as estimated from NDVI trends was due to low rainfall and an increase in soil temperature. Furthermore, the anomaly maps for precipitation, soil moisture, and LST help us identify the locations and areas of potential concern regarding reduced crop harvest. We found that large areas of the central highland agricultural farms where people largely depend on rain-fed farms are of major concern due to recurrent drought incidents. Moreover, NDLI has a high correlation with NDVI, precipitation, LST, and soil moisture and successfully captured historical droughts (Figure 8). Additionally, the results of multilinear regression indicate that NDLI, NDWI, and MEI play a significant role in the variability of vegetation health. The analysis shows that using the radiances of green, red, and SWIR, a simplified crop monitoring model with satisfactory accuracy and easiness can be developed. Thus, NDLI can be a tool to help us better understand the vegetation vigor and moisture availability, and subsequently effectively assess large-scale temporal and spatial characteristics of drought.

This analysis can serve as an important input for food security studies and the planning of potential relief measures. However, this approach suffers from the low spatial and temporal resolution satellite images utilized, as this hugely impacts the quality of the trend analysis. Further research on detecting and assessing temporal and spatial trends is needed to offer essential information for planning agencies and government policies to monitor factors that trigger drought and to minimize their impact.