The climate hazards group database provides free daily climate data with a 5 × 5 km   2 spatial resolution with quasi-global coverage (50°S–50°N, ftp://ftp.chg.ucsb.edu/pub/org/chg/products/). Daily climate data is available for the study area (Figure 1A) from CHIRPS (Climate Hazards Group InfraRed Precipitation with Stations) for precipitation (Funk et al., 2015) and from CHIRTS (Climate Hazards Group InfraRed Temperature with Stations) for temperature (Verdin et al., 2020) (Figures 1B–D). The data were generated in several stages by blending satellite records and in situ station data and are available since 1981 for precipitation and from 1983 to 2016 for temperature. The data has been evaluated for its ability to reliably reproduce the climatology and major meteorological systems of Ethiopia (Dinku et al., 2018; Belete et al., 2020; Taye et al., 2020; Kabite Wedajo et al., 2021; Malede et al., 2022) and other regions (Zambrano-Bigiarini et al., 2017; Saeidizand et al., 2018; Muthoni et al., 2019; Muthoni, 2020).

This work was conducted based on previously downscaled climate projections covering the period of 2020–2100 from a list of CMIP6 GCMs under three emission scenarios known as shared socioeconomic pathways (SSPs) (Rettie et al., 2023). The data was produced by applying a geospatial statistical downscaling technique to downscale the original coarse spatial resolution (i.e., > 70 × 70 km²) GCM outputs to 10 × 10 km² spatial resolution covering entire Ethiopia (3^oN—15^oN and 32^oE—48^oE). The data includes projections of temperature and precipitation from 13 to 9 CMIP6 GCMs, respectively (Table 1). For a detailed description of the downscaling procedure and evaluation see there. In this study, we considered the downscaled climate projections datasets under three SSPs. The selected SSPs are SSP2-4.5, SSP3-7.0, and SSP5-8.5, which represent medium, medium-high, and high-forcing scenarios based on middle-of-the-road, regional rivalry, and fossil-fueled socioeconomic development scenarios, respectively (O’Neill et al., 2017; Meinshausen et al., 2020). The emission scenarios span a broad range of CO   2 concentration, with radiative forcing level of 4.5 Wm   − 2 and 8.5 Wm   − 2 by 2100 for SSP2-4.5 and SSP5-8.5, respectively (IPCC, 2021).

Table 2 lists the 23 climate extreme indices investigated in this study. The indices were among the 27 climate change indicators which were developed by the Expert Team on Climate Change Detection, Monitoring Indices (ETCCDMI) and have been recommended by the World Meteorological Organization (WMO) (Zhang et al., 2011). The R package climdex.pcic (Pacific Climate Impacts Consortium (2020), http://cran.r-project.org/web/packages/climdex.pcic/index.html) was used to derive the climate indices at 10 × 10 km   2 grid resolution covering entire Ethiopia (i.e., 3^oN - 15^oN and 32^oE - 48^oE). To assess the skill of the GCMs in reproducing the observed climate extremes, the indices were calculated for a common 30-year period for which CHIRPS and CHIRTS data were both available. Hence, for each grid cell, the 23 climate extreme indices (Table 2) were calculated for the observed climate data (i.e., for the CHIRPS and CHIRTS data) and the historic GCMs climate outputs for the common 30-year period (i.e., 1983–2012). Similarly, the same climate extreme indices were calculated for the projected climate period (i.e., 2020–2100) for each GCM under each SSP. The indices were computed on High-Performance Cluster bwUniCluster 2.0.

The performance of the GCMs in reproducing the observed climate extremes was evaluated based on the Taylor diagram (Taylor, 2001). Taylor diagram is a widely used tool to evaluate how well a model matches observed climate states (Guo et al., 2018; Rao et al., 2019; Li et al., 2021; Yang et al., 2021; Liu et al., 2022). A Taylor diagram simultaneously visualizes three summary statistics: the standard deviation ( σ ) of simulated ( Y ) data normalized to that of the observed ( X ) data, the correlation coefficient (r) between simulated and observed data, and the centered root mean squared error (RMSE c ) between simulated and observed data (Table 3). For the Taylor diagrams, we computed the model skills of climate extreme indices over the area-averaged data across the entire grid cells for 30 years (i.e., 1983–2012).

We used the Mann-Kendall (MK) method (Mann, 1945; Kendall, 1962; Pohlert, 2020) to test the trends in the respective climate extreme indices for both the baseline climate and future projections. The MK is a non-parametric test (i.e., the data does not have to meet the normality assumption) and widely used method because of its simplicity (Cattani et al., 2018; Esayas et al., 2018; Afuecheta and Omar, 2021; Li et al., 2021; Simanjuntak et al., 2022). The MK test determines the presence of monotonic (i.e., consistent) increasing or decreasing tendency of data in a given time. The magnitude of the trend is determined by using Sen’s slope estimator (Sen, 1968) which is a non-parametric approach to estimate the overall slope in a data series (Beyene et al., 2022; Malaekeh et al., 2022; Pervin and Khan, 2022). All the MK tests and slope estimates were computed using the “trend” R software package (Pohlert, 2020).

Climate change projections usually involve three main sources of uncertainty, namely: uncertainty due to GCMs (M   t ), SSPs (S   t ), and internal climate variability (V). Using the method proposed by Hawkins and Sutton (2009); Hawkins and Sutton (2011), we evaluated the projected climate extreme indices from the GCMs (Table 1) under the three emission scenarios (SSPs) for the period from 2020 to 2100. For the temperature-related indices, this gives a total of 39 projections from the 13 GCMs ( N m = 13) and 3 SSPs ( N s = 3). For the precipitation-related indices, the total was 27 with N m = 9 and N s = 3, because three GCMs simulations did not include precipitation.

The decomposition of the uncertainty was computed on the changes (denoted as Y s , m , t ) of climate extreme indices between a future and a baseline climate period for each grid cell following the studies of Hawkins and Sutton. (2011); Zhang and Chen (2021). Here, the period 1983–2012 was considered as a baseline period to calculate the changes ( Y s , m , t ) due to climate change. The changes in temperature derived indices are expressed in terms of absolute changes as follows: Y s , m , t = p s , m , t − 1 30 × ∑ t = 1983 2012 p s , m , t (1) and the changes in precipitation-related extreme indices are defined as percentage ratio: Y s , m , t = 100 × p s , m , t 1 30 × ∑ t = 1983 2012 p s , m , t − 1 % (2) where s = 1, … , N s , m = 1, … , N m , and t = 1, … , N t refer to the number of SSPs, GCMs, and years, respectively, and p s , m , t refer to the projected extreme indices for the s t h SSPs, m t h GCM, and t t h year.

Subsequently, the resulting changes ( Y s , m , t ) were then smoothed into rolling decadal (i.e., 10 years) means and subjected to uncertainty decomposition. The uncertainty decomposition procedures are summarized as follows. To quantify and decompose the uncertainty, the smoothed mean change ( Y s , m , t ) for all GCMs and SSPs was partitioned into a climate change signal (the smooth fit, i s , m , t ) and a residual ( ε s , m , t ) by fitting fourth-order and second-order polynomial models to precipitation and temperature indices, respectively (Hawkins and Sutton, 2011; Zhang and Chen, 2021). Y s , m , t = i s , m , t + ε s , m , t (3)

The respective means of i s , m , t and ε s , m , t were calculated as follows (Eqs. 4–6); (Eq. 7): i ¯ s , t M = 1 N m ∑ m = 1 N m i s , m , t (4) i ¯ S m , t = 1 N s ∑ s = 1 N s i s , m , t (5) i ¯ S , M t = 1 N m × N s ∑ s = 1 N s ∑ m = 1 N m i s , m , t (6) ε ¯ = 1 N m × N s × N t ∑ s = 1 N s ∑ m = 1 N m ∑ t = 1 N t ε s , m , t (7)

The component of the uncertainty due to the GCMs was then estimated as the variance of the multi-scenario averages: M t = 1 N m ∑ m = 1 N m i ¯ S m , t − i ¯ S , M t 2 (8)

Likewise, the component of the uncertainty due to the SSPs was then estimated as the variance of the multi-model averages: S t = 1 N s ∑ s = 1 N s i ¯ M s , t − i ¯ S , M t 2 (9)

On the other hand, the uncertainty due to the internal climate variability corresponds to the variance of the residuals from the fits over all GCMs, SSPs, and projection period (i.e., 2020–2095, N t = 75). V = 1 N m × N s × N t ∑ s = 1 N s ∑ m = 1 N m ∑ t = 1 N t ε ¯ − ε s , m , t 2 (10)

In our present study, we assumed no interaction effects between GCMs and SSPs, and equal weights were given for all GCMs despite the different performances (Zhang and Chen, 2021). The total uncertainty ( T t ) is, therefore, the sum of the uncertainty due to GCMs, SSPs, and internal climate variability. T t = S t + M t + V (11)

The signal-to-noise (S/N) ratio was calculated to understand the influence of uncertainties on projected climate extreme indices over time (Hawkins and Sutton, 2009). S / N t = i ¯ t 1.65 × T t (12)

A larger S/N ratio implies that the projected uncertainties are smaller relative to the average climate change signal (Hawkins and Sutton, 2011; Zhang et al., 2011). The uncertainty of projected climate extreme indices was analyzed for the nine sub-regions (denoted as R1-R9) of Ethiopia that were identified based on homogeneous rainfall zones by Rettie et al. (2023).

Measuring the performance of the climate models in capturing the observed climate extreme indices is an important part of climate studies. Figure 2 shows the Taylor diagrams comparing the performance of the GCMs and their ensemble average in reproducing the different temperature extreme indices in the observation data (i.e., CHIRPS and CHIRTS) during 1983–2012. For most of the temperature-related extreme indices, the correlation coefficients range from 0.2 to 0.60. A larger range of correlation coefficients is found for extremes such as TN10p, TX90p, TN10p, TN90p, SU, and TR (refer Table 2). For absolute temperature indices like TNn, TXn, WSDI, and CSDI, the correlation coefficients (r) are very low. However, the interannual variation expressed by the standard deviations (SD) of the observation was well reproduced by most of the GCMs. Most of the indices simulated were well below 1.5 SD relative to the SD of the observations. For TX10p, WSDI and CSDI indices, many of the GCMs produced lower variability compared to the observations. The diagrams also show that the centered RMSE is well below 1.5 units relative to the SD of observations. The smallest error was found for TX90p and TX10p. The centered RMSE quantified the differences in two fields, in our case the indices in the observation data and the indices simulated by the GCMs. The diagram further shows that the GCMs ensemble average reduced both the errors (centered RMSE) and the interannual variability (SD) while the correlation coefficients (r) were increased for majority of the temperature-related indices indicating that the GCMs ensemble average performance was better than individual GCMs. The individual model as well as their ensemble average performance also varied at the regional level (Supplementary Figures S1, S2). In their multi-model study in China, Wei et al. (2022) also reported that GCMs ensemble averages of climate extreme outperform those by individual models. The models’ skill was much better for the drier and hotter sub-regions (i.e., R1, R2, and R4) compared to the cooler highland regions. More importantly, the improvement from the GCMs ensemble average was more pronounced at the regional scale (Supplementary Figure S2).

Figures 3–5 show the spatial distribution of annual trends in temperature-related extreme indices in the projected climate (2020–2100) under the three SSPs in comparison to the corresponding trends in the observed climate (1982–2012). The figures clearly show significant trends in both observations and projections for all indices except for projected cold spell duration indices (CSDI). Looking at the percentile indices, higher warm extreme indices (TX90p and TN90p) are expected in future climate ranging from 4% (SSP2-4.5) to 10% (SSP5-8.5) per decade compared to approximately 3% per decade increase in the observation period. Based on the CHIRPS data (i.e., the same observed data as in our study), Gebrechorkos et al. (2019b) reported similar patterns of extreme temperature trends for the pasty climate (1979–2010). Changes of higher magnitude of temperature extremes were reported at local level studies (Birhan et al., 2022). The increasing trend in warm extreme indices is confirmed by increasing trends in other warming indicators such as summary days (SU), tropical nights (TR), and warm spell duration (WSDI). The number of summer days (SU) and tropical nights (TR) are expected to increase by up to 35 days and 50 days, respectively, particularly in the highland regions under the SSP5-8.5 scenario. Similar increasing trends are expected for the warm spell duration index (WSDI), which could increase by up to 45 days per decade under the SSP5-8.5 scenario. The results also indicate that the magnitude the warm percentile indices expected to increase would be higher than the magnitude the cold percentile indices are expected to decrease. The relatively stronger downward trend of cold extreme indices (i.e., TX10p and TN10p) over the observation period (∼4% per decade) is expected to decrease in future climate by ∼1% per decade under the SSP5-8.5 scenario. On the other hand, the absolute extreme maximum and minimum temperature indices (TXx, TXn, TNn, and TNx) show spatially heterogeneous patterns, not solely, but particularly in the observations. However, largely significant increasing trends of the absolute extreme maximum and minimum temperature indices are expected in the future climate, especially under the higher emission scenarios (SSP3-7.0 and SSP5-8.5). For instance, under the SSP5-8.5 emission scenario, the predicted absolute maximum of maximum temperature (TXx) and minimum of minimum temperature (TNn) trends could reach up to 6 °  C (8 °  C) per decade in the warmest (northeastern) region of the country (see Figure1). This is a very relevant finding since a few incidences of such magnitude of extreme temperatures have harmful effects on crop growth (He and Chen, 2022) and would have a devastating impact on crop yields (Vogel et al., 2019). As to be expected, the projected trends consistently increase from lower (SSP2-4.5) to higher (SSP5-8.5) emission scenarios.

The Taylor diagrams in Figure 6 indicate the performance of the individual GCMs with respect to the precipitation-related extreme indices in the observation data (i.e., CHIRPS) during 1983–2012. The GCMs had difficulties reproducing the observed data even at the sub-regional level (Supplementary Figures S3, S4). The correlation coefficients (r) of most precipitation-related extreme indices were below 0.4. The simulated variability in precipitation was considerably higher than the observed one where the standard deviations normalized to the observations are mostly around 2, with R20mm, R95pTOT, and R99pTOT showing values up to 3. The errors between the GCMs and the observation (centered RMSE) were within 2 units. The diagrams reveal that the GCMs ensemble average reduced both the error (centered RMSE) and the interannual variability (SD). However, for the majority of the indices, the correlation coefficients did not improve. Overall, regarding the precipitation-related indices, the GCMs performed much less well than for the temperature-related indices. The spatial distribution of the trend in precipitation-related extreme indices in the observed climate (1982–2012) and the projected climate (2020–2100) under the three SPSs are presented in Figures 7–9. The spatial patterns show that the observed trends were largely statistically significant (i.e., shown as areas with patches) across most of the indices except for R1mm, R10mm, R20mm, and Rx5day for some pocket areas in the western and southeastern parts of the country. The southeastern part of the country which belongs to the driest regions of the country [see Figure 1A; also refer to Rettie et al. (2023)] exhibited a significant increase in maximum 5-day precipitation (Rx5day), which is equivalent to an increase of 10 days per decade. Beyene et al. (2022) also found a significantly increasing trend in Rx5day the southern region, mainly Omo-Gibe and Rift Valley lake basins. Whereas the western part bordering Sudan showed a significantly increasing number of wet days (R1mm) by roughly 5 days per decade. The trends in number of wet days were however not significant in all scenarios. Largely increasing trends were also reported in Beyene et al. (2022) for R10 mm indices for Ethiopia of which 20% of the grids were statistically significant. Cattani et al. (2018) analyzed seasonal rainfall variability and trends over East Africa for 1983–2015. They found that R1mm and R20 mm show decreasing trend during October-December and an increasing trend during March-May seasons for the larger part of East Africa. On the other hand, as in the observed climate, significantly increasing trends are projected for maximum 5-day precipitation (Rx5day) where the trends were significant across the large parts of the country. Likewise, generally increasing trends were projected for number of days with more than 10 mm (R10 mm) and 20 mm (R20 mm) indices. They were significant across large parts of the country. Regarding the precipitation totals and wet days, exceptionally higher increases were projected under SSP3-7.0 in a small pocket region in the northwestern part of the country, which was statistically significant as well (Figure 9). The projected change for the indicated pocket region was more than 25 days per year increase both for very wet (R95pTOT) and extremely very wet (R99pTOT) days. Our results are in line with Gebrechorkos et al. (2019a), who found a significantly increased trend in R99pTOT and R99pTOT in southern Ethiopia for the period of 1981–2010. Unlike the temperature indices, the trends in precipitation are generally heterogeneous in terms of spatial distribution, magnitude, and statistical significance as well as across the different emission scenarios (SSPs). Despite the large spatial inconsistency, the overall results suggested that there are indications for an increase in the frequency of the most intense precipitation events.

We assessed the three components of climate change uncertainty (i.e., model (GCM), scenario (SSP), and internal climate variability) for all climate extreme indices for the 9 sub-regions defined by Rettie et al. (2023); Figures 10–13 present the evolution of the three components over time for the temperature-related indices (Figure 10; Figure 11) and precipitation-related indices (Figure 12; Figure 13). The figures show that, for temperature-related indices, a general decreasing contribution of uncertainties from the GCMs and internal variability to the total uncertainty while the reverse was true for the uncertainty from SSPs. These findings were largely consistent across the different regions of the country. Accordingly, the uncertainty from the GCMs accounts for about 64%–88% of the total uncertainty at the beginning of the projection period (i.e., 2020). This proportion decreased to about 45%–67% by the end of the projection period (i.e., 2100). Meanwhile, the contribution from the internal variability decreased from about 11%–32% at the beginning of the projected climate to less than 15% by the end of the century. Summer days (SU) are an exception here, with a slight projected increase. On the other hand, the fraction of uncertainty from the SSPs increased from less than 1% in 2020 to 18%–54% by the end of the century. Our results are consistent with previous studies that in the beginning projections are usually dominated by uncertainties from GCMs and internal variability (Zhang and Chen, 2021). The results were consistent across the temperature indices except for the cold spell duration index (CSDI) where the SSPs’ contribution to the total uncertainty remained negligible.

In contrast to the temperature indices, the total projection uncertainty of precipitation-related indices was dominated by the contribution from the GCMs and the internal climate variability, with marginal contribution from the SSPs. However, the fractional contribution from the GCMs and internal variability varies among the different indices. For CDD, CWD, R1mm, and Rx1day, a large proportion (∼57–87%) of the total uncertainty was due to the internal climate variability whereas, for the rest of the precipitation indices, the contribution from the GCMs is considerable. Mendoza Paz and Willems. (2022) also found that the larger proportions of the uncertainty in the projected precipitation extremes were related to the GCMs. In addition, the contribution of model uncertainty increases with a lead time for most of the precipitation indices except for CDD, CWD, R1mm, and Rx1day. Results were consistent across the different sub-regions. Compared to the temperature indices, internal climate variability was an important source of uncertainty for precipitation-related indices. Previous studies also showed the relative importance of the uncertainty from the internal climate variability for precipitation-related indices (e.g., Hawkins and Sutton, 2011; Fatichi et al., 2016; Gu et al., 2018).

To reduce the high uncertainty associated with the precipitation indices, we suggest further investigation based on multiple reference data sets (Grusson and Barron, 2021; Madakumbura et al., 2021) and as far as possible, station-based observation data (Cattani et al., 2018; Kim et al., 2019). Fatichi et al. (2016) claimed that the rigorous assessment of historic climate variability may give sufficient information about future changes in precipitation extremes. In addition, seasonal level analysis (Cattani et al., 2018; Ademe et al., 2020; Gemeda et al., 2021; Ali Mohammed et al., 2022; Beyene et al., 2022; Teshome et al., 2022) could also help to reduce the uncertainty compared to annual level analysis (this study) given the high spatial variability in the country.

Finally, we quantified the signal-to-noise ratio (S/N) to demonstrate the influence of uncertainties on projected climate extreme indices over time (Hawkins and Sutton, 2009) and hence to evaluate the robustness of the projected changes in climate extremes (Hawkins and Sutton, 2009; 2011; Zhang et al., 2011). Figure 14 and Figure 15 present the S/N ratio for the temperature and precipitation-related indices, respectively, for the different sub-regions.

Largely consistent across the temperature-related indices, the signal-to-noise values increased in general with time, but with considerable regional variation (Figure 14). This implies that the magnitude of projected changes was greater than the magnitude of the associated uncertainty and hence the projected changes are reliable. However, the S/N ratio reaches peak values by the mid of the century (i.e., between 2050–2060) for temperature intensity indices (TXx, TXn, TNx, and TNn) with slightly decreasing values with projection time. The peaks around the mid of the century could be related to the shift in the contribution of uncertainties from the different sources (i.e., GCMs, SSPs, and internal climate variability). Zhang et al. (2011) reported a similar period where the shift in the contribution of uncertainty among the various sources occur. Regional comparisons show that sub-region R8 (which is the wettest region in the country, cf. Rettie et al. (2023)) shows a higher S/N ratio for most of the temperature-related indices while sub-region R1 (which belongs to the driest regions in the country) shows relatively lower S/N, particularly for coldness indices (i.e., CSDI, and TX10p).

Precipitation-related extremes are key for climate change adaptation as the country’s economy is heavily dependent on rainfed agriculture. The S/N ratio was below unity despite a slight increment with projection time (Figure 15) for most of the precipitation indices indicating that the magnitude of projected changes was smaller than the associated uncertainty. This implies that the projected changes are associated with high uncertainty and make the projection less reliable and not well-suited as a basis for decision-making. Achieving reliable projections for precipitation has been a challenge due to its associated higher uncertainties compared to temperature (Madakumbura et al., 2021; Zhang and Chen, 2021; Birhan et al., 2022). On the other hand, located near the equator and the Indian Ocean, the effects of the bi-annual migration of the Inter-Tropical Convergence Zone (ITCZ) and the El Niño–Southern Oscillation (ENSO) are the most important climate systems governing precipitation across Ethiopia (Korecha and Barnston, 2007). In this regard, one of the limitations of the current study is attributed to the structural deficiency of the climate models used in the study in reasonably simulating these major climate systems. The majority of state-of-the-art GCMs fail to simulate realistic ENSO characteristics (Beobide–Arsuaga et al., 2021) and the double- ITCZ bias remains one of the most outstanding errors in the models (Tian and Dong, 2020). Given these model deficiencies, their long-term prediction of climate extremes might be affected as well. Therefore, the results of our study particularly those of precipitation-derived extreme indices should be taken with caution. Despite these limitations, the regional comparisons suggest that fewer homogeneous clusters could be sufficient for the kind of studies treated in this paper (Ware et al., 2022).