Spatiotemporal variability of agricultural drought and its association with climatic variables in the Upper Awash Basin, Ethiopia

Drought is a serious threat to agriculture in Ethiopia. This study examined the spatiotemporal variability of agricultural drought and its association with climatic variables in the Upper Awash basin. Mann–Kendall (MK) trend test was employed to examine the drought trend while Sen’s slop estimator and pixel-based linear regression model were used to analyze the magnitude of drought changes. The association between agricultural drought and climatic variables was evaluated by the Pearson correlation coefficient (r). High spatiotemporal variability of drought was observed in Kiremit (June–September) and Belg (February–May) seasons. The Belg season spatial average vegetation condition index (VCI) trends were decreased insignificantly from 2001 to 2019 at a 5% significant level, whereas the spatial average VCI trends of Kiremit season were increased insignificantly. The return period of severe droughts during the Belg season was less frequent than the Kiremt season severe drought. The correlation between spatial average VCI and precipitation was positive for Belg and Kiremit seasons. Likewise, the correlation between average VCI and land surface temperature (LST) was negative in Belg and positive in Kiremit season. Moreover, the correlation between mean VCI and Pacific Ocean Sea Surface Temperature (SST) was positive for Belg and Kiremit seasons. The influencing factor of precipitation and LST on VCI during Belg season was higher than Kiremit season. The findings of this study are vital for decision-making systems and preparing plans to adjust sowing time, select drought-resistant crops, practice in situ water conservation, practice small-scale irrigation and diversify the income of smallholder farmers.


Introduction
Drought is one of the most devastating natural hazards [1], which are affecting ecosystem functions and services [2,3]. Depending on the duration and impact of the water cycle component, droughts can be classified into a meteorological drought (e.g., no precipitation over a long period), hydrological drought (e.g., insufficient surface and groundwater water supply), and agricultural drought (e.g., lack of availability of water for a plant or vegetation growth) [4,5]. Drought can result in soil degradation, desertification, water shortages, plant death, and fires [6][7][8]. Drought can also affect crop development, agricultural and socio-economic activities, and also contributes to social crises and political problems [9][10][11]. Drought-affected areas develop slowly as the signs of plant moisture stress often change gradually [12]. Thus, understanding the extent and frequency of drought and its relationship with climatic variables is imperative to improve agricultural productivity and ensure sustainable socio-economic development [13].
Several methods are available in the literature for monitoring drought [14][15][16]. Of these methods, the traditional drought detection methods are based primarily on precipitation, soil moisture, temperature, evaporation, and surface runoff [17]. In this case, spatial interpolation can be obtained from point data. However, several factors affect the interpolation process, so there may be some uncertainty in the interpolation results. These include questions about sample size, spatial distribution, and missing data [13,15,16]. As a result, satellite remote sensing has become one of the important techniques for earth observation and drought monitoring and can provide comprehensive information on the dynamics and processes of terrestrial systems [18]. Recently, remote sensing is known to be a powerful tool for assessing the spatiotemporal variability of drought [19][20][21].
Globally, numerous remotely sensing vegetative drought indices have been established to monitor droughts, such as Normalized Difference Vegetation Index (NDVI), Vegetation Condition Index (VCI), and Vegetation Temperature Condition Index (VTCI) [22][23][24][25][26]. Of these indices, the NDVI is easy to calculate and can be used to consider long-term series. But it is less reliable when monitoring the drought in the heterogeneous area. This is due to the effect of geographical location, ecological system, and soil conditions [27,28]. To overcome these problems, Kogan [27] developed VCI for a specific region. The VCI reflects the overall effect of rainfall, soil moisture, weather, and agricultural practices. Hence, in areas such as the Upper Awash basin, which has different ecosystems and heterogeneous topography, VCI is vital to compare the impact of weather in areas with different ecological and economic resources, as the index captures rainfall dynamics better than NDVI, particularly in geographically heterogeneous areas. For this reason, it is considered an ideal indicator for large-scale drought monitoring [22,27,29]. In earlier times, VCI was applied in drought monitoring and analysis [13,[30][31][32][33], and its reliability has been verified by many studies [16,[33][34][35][36]. As a result, VCI was used in this study as an indicator of drought [13]. The VCI can describe the spatiotemporal variation in land cover and vegetation. It also helps to distinguish meteorological impacts on vegetation [16].
In East Africa in general and in Ethiopia in particular, droughts have recurred in the last decades and is a major natural disaster that contributes to food insecurity and poverty [37][38][39]. In Ethiopia, the availability of drought has increased due to climate change and will cause a decline in water and agricultural production [40]. The economy of more than 85% of the Ethiopian population mainly depends on rain-fed agriculture, which is vulnerable to climate change [41]. According to Funk et al. [42], 2015 was the driest year in most parts of Ethiopia. As a result, the main rainy season was late and below normal [41]. The El Niño-Southern Oscillation (ENSO) was the ultimate cause of this drought. More specifically, the warmer stage (El Niño) is closely associated with reserve rainfall during the main rainy season in central Ethiopia, including the Upper Awash basin [43]. In such conditions, the moisture requirements of plants cannot be met, leading to a sharp decline in plant production. Drought is one of the most frequent environmental threats in the Upper Awash basin of Ethiopia [44]. Besides, Edossa et al. [45] has also reported the existence of extreme meteorological drought events in the Upper Awash basin. Based on the socio-economic analysis, Desalegn et al. [44] also indicated that the Upper Awash basin experiences drought every two years. However, the spatiotemporal variability of agricultural drought and its association with climatic variables are not well-understood in the Upper Awash basin. Besides, the spatial and temporal variation of agricultural drought across agro-ecological zones in the basin is not known. This type of study has immense importance for developing an understanding of the basics of basin vegetation dynamics and thus helps in policymaking. Evaluation of spatiotemporal variability of agricultural drought and its association with climatic elements is important for policy-makers and planners for establishing effective and comprehensive monitoring and early warning system to reduce the adverse impacts of drought. Therefore, this study was aimed to examine the spatiotemporal variability of agricultural drought and its association with climatic elements (such as precipitation, LST and Pacific Ocean SST) in the Upper Awash basin of Ethiopia. Specifically, the objectives of this study were: i) to assess the spatiotemporal variability of agricultural drought, ii) to examine the frequency of agricultural drought, and iii) to evaluate the association between VCI and climatic variables in the Upper Awash basin.

Description of the study area
The Upper Awash basin is found in the central part of Ethiopia mainly at the western margin of the Main Ethiopian Rift system. Geographically, the basin is located between 8°16′N-9° 18′ N and 37° 57′E-39°17′E, covering 10,640 km 2 . Its physical settings are characterized by the heterogeneity of the large natural systems such as orographic, the high plains, mountains, and plateaus [46]. The topography undulates between 1587 to 3561 m a.s.l. The climate of the study area is influenced by undulating mountain chains and the circulatory systems that interact with orography, cross-equatorial wind system, and the movement of the Inter-Tropical Convergence Zone (ITCZ) [46]. The basin's climate is humid at the highlands and arid to semiarid in the escarpment and rift valley [47]. The land-use types are intensively cultivated (67%), moderately cultivated (25.5%), bushland or shrubland or wooded grassland (4.5%), and urban area and alpine vegetation (3%) [48]. Hence, agriculture is the main economic activity. The mixed crop-livestock system is the most important farming system in the plateau and the highland areas. This farming system is dominated by smallholder farmers in rain-fed agriculture with limited but complementary livestock production. Landholdings are generally small and are fragmented into many plots with different land uses. Besides, the main towns found in the Upper Awash basin are characterized by industry and agro-industry.  [49] and Zaroug [50], the NINO 3.4 SST data has NINO3 and NINO4 characteristics. Therefore, the Pacific Ocean SST of the NINO3.4 area was also used to evaluate the relationship between the SST and the index of vegetative drought.

Remote sensing data processing
The obtained MOD13Q1 and MOD11A2 data are supported with a sinusoidal projection method in the Hierarchical Data Format (HDF). Through format conversion and reprojection using MODIS Reprojection Tool (MRT), these data were prepared for the Geographic Information System (GIS) program. Using metadata attributes stored in the dataset, the filtering of no data values and cloud deletion from the imagery was performed [51]. MODIS NDVI for sixteen days and MODIS LST datasets for eight days were composite of maximum daily values throughout the year. The noise of these composite datasets was removed using the fast Fourier transform (FFT) algorithm [51]. Using the forward transformation process, the Fourier transform transforms the spatial domain image into a frequency domain image. Then, using the frequency domain images filtered through the inverse Fourier transform (IFT), an enhanced noise-free image was created [52]. For the two-dimensional cases, the Fourier transform representation for the discrete f unc tio ns i s e xpr essed as the sum of sine and cosine weight and is given in Eq. 1: where u and v are spatial frequencies, F(u, v) is a function of the frequency domain, f(x, y) is a function of the spatial domain, i is √ −1 , N is the number of x-direction pixels, and M is the number of y-direction pixels. X ranges from zero to N-1 and y ranges from zero to M-1. Fourier se que nces (Fu) have b een mu ltipl ied by a low pass filter (w), giv ing th e freque ncy domain a filtered signal (wFu). The enhanced image in the spatial domain was reconstructed using the inverse transform after filtering to remove noises related to highfrequency components (Eq. 2): The MODIS NDVI was converted into required NDVI values between −1 and 1 using the Arc GIS environment using a raster calculator by multiplying the improved image data with a scale factor (0.0001). The total numbers of NDVI images downloaded (2001-2019) for Belg (February-May) and Kiremit (June-September) seasons were 304. Since the MODIS data are available in every 16 days composite, it was converted into monthly solutions. Two 16 days composite images were added and divided by two in the raster calculator to get the mean monthly NDVI for each season. Similarly, the LST of 2001-2019 periods for Belg and Kiremit seasons were derived from eight days of composite MOD11A2. The digital number (DN) of LST was converted to degree Celsius ( • C) by multiplying the input digital number (image) with a scale factor (0.02) and then subtracting 273.15 °C. A total of 608 LST images were acquired (2001-2019) during Kiremit and Belg seasons. Then, four eight-day composite images were added and divided by four in the raster calculator to get the mean monthly LST for each season.
Most precipitation data from in situ meteorological stations within the Upper Awash basin had an outsized percentage of missing data problems. Moreover, the spatial distributions of stations were not evenly dispersed. During this case, the CHIRPS satellite (https:// data. chc. ucsb. edu/ produ cts/ CHIRPS-2.0/) is an important source of precipitation data [53][54][55][56]. To determine the accuracy of the CHIRPS satellite precipitation data in the Upper Awash Basin, 19 meteorological gauge station data from the National Meteorological Service Agency (NMSA) of Ethiopia were used as references. The details of these stations are given in Table 1. By taking the average of the previous and subsequent months, monthly missing values were filled. However, those years with an entire missing data were omitted from the analyses [53,57]. In this study, a surface map in the form of a grid precipitation map for the study area was constructed using the ordinary kriging geostatistical interpolation method. Ordinary kriging is the best linear unbiased estimator and was found to be the best method because it produces little root mean square error [58]. Besides, the Pearson correlation coefficient (r) was used to assess the efficiency of the areal average CHIRPS precipitation product by using the areal average interpolated meteorological rainfall [54,59].

Methods for identifying drought
The short (Belg) and main (Kiremit) rainy seasons determine rain-fed agricultural production in Ethiopia, including the Upper Awash basin [60]. Almost 85% of the population in the study area practice rain-fed agriculture, which depends on the main and short rainy seasons. The seasonal characteristics of rainfall have a great influence on the production potential of crops in the rain-fed agricultural systems since the availability of water in the soil is essential for the growth of crops and vegetation. MODIS NDVI data were used to derive the Belg and Kiremit seasons VCI. VCI is most useful during the main growing (Kiremit) and short growing (Belg) seasons because it is a measure of vegetation vigor. The VCI was determined in Eq. 3 [61,62] below: where VCI is the vegetation condition index of the pixels, NDVI i is the NDVI value of the pixels, and NDVI max and NDVI min are the maximum and minimum NDVI values during 2001-2019 periods, respectively. The numerator refers to the difference between the actual value of the NDVI and the minimum value of the NDVI for a given time and is representative of the meteorological and vegetation data for a given time. The maximum and minimum denominator values represent the best and the worst vegetation growth conditions, respectively, and the difference of them reflects somewhat the local vegetation condition [22,27,28]. The VCI includes the NDVI with both real-time and historical details. VCI results range between zero and 100, where lower VCI values imply poor growth of vegetation and higher drought levels [27,28,63,64]. On the other hand, higher VCI values are an indicator of good vegetation conditions and characterize lower drought events. Based on VCI, droughts have been classified and the spatial and temporal variations of drought over 2001-2019 time spans were analyzed. In the present analysis, three forms of VCI drought classes that are defined by Kogan [28] were used (Table 2).

Trend detection of the VCI drought index
The pixel-based linear regression model was used to analyze the spatial and temporal trends of VCI changes pixel-wise using the Belg and Kiremit seasons datasets (2001-2019) (Eq. 4). The positive VCI slope values represent an upward trend while the negative VCI slope values indicate a downward trend [33]. An upward trend of VCI indicates enhanced vegetation growth and drought reduction while a downward trend of VCI specifies a reduction of vegetation cover and increasing drought. The trend in the VCI was calculated using Eq. 4 [13,65]: where VCIi is the vegetation condition index for year i, n is the duration of the time series (n = 19), and ti is the index number for the years 2001 to 2019 (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19).
In addition to the pixel-based linear regression model, Sen's slope was used to estimate the magnitude of temporal shifts of the areal average VCI drought index. In contrast to linear regression, this approach is less affected by missing values and outliers [59,66]. If there is a linear trend, the magnitude of the monotonous trend can be determined using the nonparametric slope estimator of Sen using Eq. 5 [67]: where β is the median value of the slope values between the yi and yj data measurements in phases i and j (i < j), respectively. The positive value of β reflects an upward trend, while a downward trend is shown by the negative value of β. The sign of β represents the course of the data trend, while its value shows the trend's steepness.
The Mann-Kendall trend test was used to determine trends in the time series of areal average VCI values for the whole basin. This approach is less influenced by missing values and unequal distribution of data and less vulnerable to outliers because the ranks of the observations are taken into account rather than their actual values [55,68]. The null hypothesis ( H0) of no trend, that is, the Yi observations are randomly ordered in time, against the alternative hypothesis ( H1), according to the Mann-Kendall trend test, where a monotonic (increasing or decreasing) trend was checked in the time series. MK statistics S are computed based on Yue et al. [68] and Mann [69] using Eq. 6: where Yi and Yj are sequential data values for n-length data of the time series and If the dataset is distributed identically and independently, then S's mean is zero, and S's variance is given by Eq. 8: where n is the dataset length, m is the number of tied groups in the time series (a tied group is a collection of sample data with the same value), and ti is the number of data points in the ith group.
The Z statistics are calculated using Eq. 9: To test either an upward or downward monotone trend, a significance level alpha (α) was used. By comparing the computed Z with critical values, the decision for the twotail test was made. The null hypothesis is rejected if the computed Z absolute value is greater than the critical value, or if the p-value is less than the selected significance level (α = 0.05 or 0.1). The direction of trends is upward for positive Z-value and downward for negative Z-value when the null hypothesis is rejected [70]. The result is said to be statistically significant if the null hypothesis is rejected as in Figs. 1 and 2.

Exceedance probability and return periods
Using Weibull's frequency distribution equation, the VCI exceedance probability and return periods were for S > 0 determined [71]. The return period and the probability of exceeding each other are reciprocal (Eq. 10).
where p(xm) is the probability of exceedance, Tr is the return duration that implies an average number of years that will be equalled or exceeded during a given case, n is the total number of study years (19 years), and m is the rank of observations in descending order.

Correlation analysis of the VCI and climate variability
The degree and direction of the relationship between two variables are expressed in the Pearson correlation coefficient (r) [33,72]. The greater association between the two variables is demonstrated by a larger absolute value. To evaluate the relationship between drought and climatic factors (precipitation, LST, and SST), the Pearson correlation coefficient was used. The coefficients of Pearson correlation were determined using Eq. 11 [73]: where r is the coefficient of correlation, n is the time series duration and i is the number of years during the periods studied (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19). Xi and Yi are the VCI and the value of climate variability in the year i, respectively, and X and Y are the mean VCI and the mean climate variability during the study periods, respectively.

Evaluation of CHIRPS precipitation data
The results of the comparison between the areal average CHIRPS and areal average interpolated meteorological  Fig. 3a and b. A very good agreement was observed between areal average observed rainfall and CHIRPS precipitation product ( Fig. 3a and b). The findings of this study agreed with a study carried out by Ayehu et al. [74] in the Upper Blue Nile basin; Belay et al. [55] in the Beles basin; Alemu and Bawoke [59] in the Amhara regions of Ethiopia; and Dinku et al. [54] in Eastern Africa.

Belg and Kiremit seasons yearly variations of drought index (VCI)
The yearly variations of the VCI drought index from 2001 to 2019 periods during the Belg and Kiremit seasons are shown in Fig. 4. The lines parallel to the X-axis in Fig. 4 represent the threshold values of VCI for severe drought (VCI ≤ 35%) and drought (35% ≤ VCI ≤ 50%). During the  (Fig. 4).
The Belg and Kiremit season spatial mean VCI from 2001 to 2019 periods was 45.4% and 42.0%, respectively. The Kiremit season VCI was increased insignificantly at the rate of 0.395%yr −1 over the whole basin while Belg season VCI was decreased insignificantly at the rate of 0.693%yr −1 over the basin at a 5% significant level (Table 3). It indicates that the trend of drought in the Upper Awash basin was decreased and increased insignificantly during Kiremit and Belg season, respectively. Similarly, Shen et al. [16] indicated that the mean annual VCI of China from 1982 to 2010 was slowly increased, indicating that the enhanced vegetation growth and the drought alleviated. Liang et al. [13] also reported the increased trend of VCI in China from 1981 to 2015 during the spring, summer, and autumn seasons, indicating that drought was decreased in China during these periods (1981-2015). On the contrary, Gidey et a1 [75]. reported that the Vegetation Health Index (VHI) value in

Belg and Kiremit season yearly variation of vegetation condition index (VCI) across agro-ecological zones (AEZs)
Agro-ecological zoning can be characterized as a spatial division of the landscape into comparatively similar agricultural and ecological features [76], providing a context for understanding the complexity of agro-ecological systems [77]. Topography plays an important role in the agricultural zone in mountainous countries such as Ethiopia, where Africa's most prominent mountain system is located [76]. The elevation-based agro-ecological zoning developed by Hurni [76] was therefore adopted to classify the basin into three agro-ecological zones (AEZs). These AEZs are locally referred to as Weyna Dega (subtropical) (1500-2300 m.a.s.l), Dega (temperate) (2300-3200 m.a.s.l) and High Dega (Alpine) (3200-3700 m.a.s.l) [76]. Because of its good vertical and horizontal precision, this classification was performed using a 30 m SRTM DEM spatial resolution [78]. To improve agricultural development planning, agronomic zoning is generally significant, as agronomic regions are heavily dependent on climatic parameters such as rainfall quantity and variability, temperature, and vegetation characteristics [76]. It is important to understand the occurrence of drought in the key agronomic regions of the basin. This then helps to recognize the most drought-affected agro-ecological areas, which in turn will enable decision-makers to establish environmentally sound drought-friendly strategies.
Understanding the yearly variation of VCI distributions across agro-ecological zones (AEZs) enables us to identify areas experiencing severe droughts. The annual variation of vegetative drought index in different AEZs of Upper Awash basin during Belg season in the study periods (2001-2019) is shown in Fig. 5. In the Belg season, VCI was decreased insignificantly at the rate of 0.817%yr −1 and 0.552%yr −1 in Woyna Dega and Dega AEZ, respectively (Table 4). But VCI was decreased significantly at the rate of 2.174%yr −1 in High Dega AEZ (Table 4) (Fig. 6).  (Fig. 7). The drought years in Ethiopia include 1984Ethiopia include , 1987Ethiopia include , 1991Ethiopia include /1992Ethiopia include , 1993Ethiopia include /1994Ethiopia include , 2002Ethiopia include , 2008Ethiopia include /2009Ethiopia include , 2011Ethiopia include /2012Ethiopia include , 2015. These drought years either coincide or follow El Nino events [53]. Likewise, severe droughts were identified in the year 2002,2003,2009,2010,2012,2013, and 2014 during the Kiremit season (Fig. 8).  Figs. 9 and 10, respectively. In the Belg season, severe droughts were observed in the central and southern parts of the basin (Fig. 9). Likewise, severe droughts in the Kiremit season were experienced in the central, eastern, and northeastern parts of the study area (Fig. 10). Conversely, none drought areas were experienced in the northern and western parts during the Belg season and in the northern, western, and southeastern parts in the Kiremit season. In Belg season, maximum severe drought areas were observed in Woyna Dega AEZ while most none drought areas were experienced in Dega and High Dega AEZs. Similarly, most severe drought vulnerable areas in the Kiremit season were observed in Woyna Dega AEZ whereas maximum none drought areas were identified in Dega AEZ. Figure 11 shows

Spatial and temporal drought trends
The spatial VCI trend (2001-2019) of the Upper Awash basin in the Belg season was varied spatially from −6.4 to 5.54 (Fig. 12). Negative VCI slope values were observed in the mid-northern (in and around Addis Ababa), southern, southwest, and eastern parts of the study area due to the rapid expansion of anthropogenic influences such Ababa. Urban expansion due to higher population invasion and rapid socio-economic development in Addis Ababa and neighboring cities inevitably reduces natural   (Fig. 13). In the Kiremit season, the negative VCI slope values were observed in the mid-northern (in and around Addis Ababa) and northwestern parts because of urban and industrial expansion. In contrast, positive VCI slope values were experienced in the eastern and western parts of the area (Fig. 13). The highest negative VCI trends (−6.7-−2.5) were experienced in the northern parts of the basin, and it indicates the potential drought and droughtvulnerable area due to reductions in forests and shrubland areas. In the Kiremit season, negative VCI slope values have been identified in the built-up area and grassland area while negative VCI slope values have been identified in the forest and built-up areas during the Belg season. On the contrary, positive VCI slope values were identified in most of the cropland areas during the Belg and Kiremit seasons. Similarly, Qian et al. [33] find that the VCI trend was increased in most agricultural areas of China from 1982 to 2010. The negative VCI trend area coverage was higher in the Belg season than the Kiremit season in the Upper Awash basin. Therefore, most drought vulnerable areas were observed in the Belg season than Kiremit season during the study periods (2001-2019) of the Upper Awash basin. The Belg season VCI spatial average temporal trend was decreased from 2001 to 2019 (Table 3) and indicating the decline of vegetation growth and the rise of drought. On the contrary, the Kiremit season VCI spatial average temporal trend was increased through 2001 to 2019 (Table 3) and indicating the enhancement of vegetation growth and the decline of drought in the study area.

Exceedance probability and return periods
The probability of exceedance (P(xm)) of the Upper Awash basin areal average VCI equal to or greater than 35% and 50% during the analyzed years (2001-2019) in Belg and Kiremit season is shown in Fig. 14. The Belg season VCI corresponded to P(xm) of 0.65 and 0.5, respectively, equivalent to or greater than 35% and 50%. These findings indicate that in the Upper Awash basin there was a 35% chance of severe drought occurring (VCI < 35%) and a 50% chance of normal drought occurring (35% < VCI < 50%). Similarly, the average VCI equal to or greater than 35% and 50% corresponded to P(xm) of 0.7 and 0.2, respectively, in the Kiremit season, which forecasts the likelihood of severe drought occurring and normal drought to be 30% and 80%, respectively.
The return period of the spatial average Belg VCI is equivalent to or greater than 35% and 50%, respectively, Similarly, the return periods during the Kiremit season of severe drought (VCI < 35%) and normal drought (35% < VCI < 50%) were 1.43 and 5 years, respectively. Similarly, Kogan et al. [82]] reported that droughts affected the Horn of Africa (including the Upper Awash basin) annually. Besides, Gidey et al. [75] find that, from 2001 to 2015, the incidence of agricultural drought events in the low land areas of Raya was 10-11 times higher. Moreover, Edossa et al. [45] reported that meteorological droughts occurred once in two years in the Upper Awash basin (Nazareth Area), which is consistent with the findings of this study. The normal drought return cycle (35% < VCI < 50%) during the Kiremit season was less frequent than the Belg season. On the other hand, the return period of severe drought (VCI < 35%) during the Belg season was less frequent than the Kiremt season.

Association between VCI and climate variability
The Pearson correlation coefficient (r) between spatial average VCI and precipitation was 0.64 and 0.10 for Belg and Kiremit seasons, respectively ( Table 5). The positive correlation between spatial mean VCI and precipitation implies that enhanced precipitation supports vegetation growth and drought reduction [83]. Similarly, a study conducted by Tiruneh et al. [84] reported a strong positive correlation (r = 0.62) between NDVI and precipitation during the Belg season from 2001 to 2016 in the Upper Awash basin. The result of this study also reported a weak positive correlation (r = 0.10) between VCI and precipitation during the Kiremit season in the Upper Awash basin. This weak correlation between precipitation and VCI during the kiremit season may be due to the signal saturation of certain biomass values, the lack of solar radiation used for photosynthesis due to cloud. Liou and Mulualem [40] reported a strong positive correlation (r = 0.83) between NDVI and precipitation during the Kiremit season in Ethiopia from 2001 to 2018. Likewise, Measho et al. [79] reported that a strong positive correlation   [86] also reported that the linear correlation between vegetation temperature index and monthly precipitation in the southern Great Plains of the USA. The correlation between spatial average VCI and land surface temperature (LST) was negative (r =−0.77) and positive (r = 0.06) during Belg and Kiremit seasons, respectively (Table 5). These results indicate that land surface temperature (LST) was the main influencing factor on spatial average VCI during the Belg season. Similar to the finding of this research, Tiruneh et al. [84] find that the negative (r =−0.67) and positive (r = 0.41) correlations between spatial average NDVI and LST during Belg and Kiremit season, respectively, in the Upper Awash basin of Ethiopia from 2001 to 2016. The positive correlation between spatial mean VCI and LST during the Kiremit season indicated that an increase in land surface temperature caused an upward trend in the VCI, which implies a decline in drought. The result of this study was also supported by Baniya et al. [83] and reported that a positive correlation between the spatial average VCI and temperature in Nepal from 1982 to 2015 during annual and seasonal monsoon time scales. A study undertaken by Qian et al. [33] also reported the strong positive correlation between VCI and mean annual temperature in the agricultural area of China from1982 to 2010. The increased land surface temperature will enhance vegetation growth by an accelerated release of nutrients and improved availability from the soil until the optimum temperature for photosynthesis [83,87]. A similar result was also reported in this study during the Kiremit season. On the contrary, Liou and Mulualem [40] reported a strong negative correlation (r =−0.76) between NDVI and LST during the Kiremit season in Ethiopia from 2001 to 2018.
The Pearson correlation coefficient (r) between spatial mean VCI and SST was 0.32 and 0.106 for Belg and Kiremit seasons, respectively (Table 5). Likewise, Tiruneh et al. [84] find the positive correlation, r = 0.42 and r = 0.22, between the spatial mean NDVI and SST anomaly during Belg and Kiremit season, respectively, in the Upper Awash basin. Similar to the finding of this research, Philippon et al. [88] reported that a positive correlation between NDVI and NINO 3.4 SSTA in the autumn season over northwestern Africa.

Conclusion
This study investigated the spatiotemporal variability of agricultural drought and its association with climatic variables in the Upper Awash basin of Ethiopia. High The Belg season spatial average VCI trends were decreased, whereas the Kiremit season spatial mean VCI trends were increased during the studied periods (2001-2019). However, the decreasing and increasing trends of VCI were not statistically significant at a 5% significant level. The decreasing trend of VCI in the Belg season indicates the increasing trend of agricultural drought while the increasing trend of VCI in the Kiremit season shows the decreasing trend of drought in the basin. In the Belg season, the most severe drought years were identified in Woyna Dega AEZ whereas most none drought years were found in High Dega AEZ. Similarly, during the Kiremit season, the most severe drought years were observed in Woyna Dega and High Dega AEZs, while most none drought years were experienced in Dega AEZ. In the Kiremit season, areas identified with severe droughts were mostly cropland areas, while in the Belg season, the most severe drought areas were identified in most cropland areas and some builtup areas in the basin. The return period of severe drought (VCI < 35%) during the Belg season was less frequent than the Kiremt season. The correlation between spatial mean VCI and precipitation was 0.64 and 0.10 for Belg and Kiremit seasons, respectively. Likewise, the correlation between spatial average VCI and land surface temperature (LST) was negative (r =−0.77) for Belg and positive (r = 0.06) for the Kiremit season in the basin. Moreover, the correlation between spatial mean VCI and Pacific Ocean sea surface temperature (SST) was 0.32 and 0.106 for Belg and Kiremit seasons, respectively. Generally, precipitation and LST were the main influencing factors on VCI during the Belg season in the basin. Therefore, the findings of this study can be used as a useful information source on spatiotemporal variability of drought and its association with climatic variables in the drought-prone areas of the Upper Awash basin for policy-makers and planners for establishing effective and comprehensive monitoring and early warning system to reduce the adverse impacts of drought.

Funding
The research was fully self-funded.

Conflicts of interest The authors declare that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.