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

: The recent droughts that have occurred in di ﬀ erent parts of Ethiopia are generally linked to ﬂuctuations 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 di ﬀ erence 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 signiﬁcant 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, signiﬁcant countrywide droughts occurred during the El Niño 2009 and 2015 years. Furthermore, the Pearson correlation coe ﬃ cient 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 ﬁrst attempt to apply a recently developed index, the Normalized Di ﬀ erence 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 simpliﬁed crop monitoring model with satisfactory accuracy and easiness can be developed.


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 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.

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].
Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 19 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.

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 identifiedthe 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 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. 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 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. 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.

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

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 Remote Sens. 2019, 11, 1828 6 of 19 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: 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.

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: where x i is the seasonal mean of variable x, x is the long term seasonal mean and σ 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: where NDVI, NDVI min , and NDVI 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.

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 (X k , k = 1, 2, . . . , n) is given as: 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: When n ≥ 8, the statistic S is approximately normally distributed with mean E[S] = 0, and variance σ 2 given by the following equation: 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: Now, Z follows a standard normal distribution whereby its positive (negative) value indicates an upward (downward) trend. If Z is greater than Z α/2 , where α 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 |z| > 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: 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.

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 , . . . , x k and a response can be written as: and ε 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), β 0 is the intercept, and β 1 , β 2 , · · · , β 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

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].

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. 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. 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. 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.

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 Figures 6 and 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.
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

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 Figures 6 and 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.   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 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.

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. 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.

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 Figures 4, 6, and 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. 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 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. 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 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 × 10 , 2.86 × 10 , 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 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 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 × 10 −16 , 2.86 × 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 −5 + 3.75 × NDLI + 0.518 × NDWI − 1.386 × 10 −3 × MEI) was obtained via the backward elimination procedure in a stepwise regression analysis, which was achieved by dropping the least significant feature.

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.