Assessment of Climate Change Impacts on the Hydroclimatic Response in Burundi Based on CMIP6 ESMs

: Burundi is susceptible to future water-related disasters, but examining the inﬂuence of climate change on regional hydroclimatic features is challenging due to a lack of local data and adaptation planning. This study investigated the inﬂuence of climate change on hydroclimate-focused changes in the climatology of heavy precipitation (and streamﬂow) means and extremes based on the multi-model ensemble mean of earth system models in the sixth phase of the Coupled Model Intercomparison Project (CMIP). For runoff analysis, hydrologic responses to future climate conditions were simulated using the Soil and Water Assessment Tool (SWAT) model over the Ruvubu River basin, Burundi. Temperature increases by 5.6 ◦ C, with strong robustness, under future climate conditions. The mean annual precipitation (and runoff) undergoes large seasonal variations, with weak robustness. Precipitation (and streamﬂow) changes between the wet and dry seasons differ in signal and magnitude. However, alterations in both the amount and frequency of precipitation reveal the intensiﬁcation of the water cycle due to anthropogenic climate change. Thus, the highest variability in the maximum daily streamﬂow is shown in months of long wet seasons, especially in the far future (2085). Without considering the regional climate characteristics and shared socioeconomic pathway (SSP) scenarios, this behavior is expected to be enhanced in 2085 (compared with 2045) and increase the severity of extreme precipitation and ﬂood risk. Climate change will cause alterations in the magnitude and seasonal distributions of extreme precipitation (and streamﬂow). These ﬁndings could be important for ﬂood planning and mitigation measures to cope with climate change in Burundi.


Introduction
Sub-Saharan Africa is plagued by frequent extreme water-related disasters such as floods and droughts as a result of global-scale atmospheric circulation and regional climate variability. For instance, global-scale tropical atmospheric circulation (i.e., Hadley circulation) in Sub-Saharan Africa, which leads to a twice-a-year seasonal migration, causes a substantial difference in precipitation between the wet and dry seasons. Concurrently, atmospheric circulation causes extreme climate events (i.e., heavy precipitation) resulting from interannual climate variations (i.e., El Niño-Southern Oscillation and Indian Ocean Dipole) in Sub-Saharan Africa [1][2][3][4][5]. Moreover, continuing changes in Earth's climate system will enhance the probability of hard, prevalent, and irreversible impacts on humans and ecosystems across the globe (International Panel on Climate Change, IPCC) [6]. For Sub-Saharan Africa, climate projections show an increase in the variability of rainfall as the temperature rises, which might markedly affect the hydrological cycle [4,7,8]. An intensified hydrological cycle is considered to be the dominant factor influencing extreme climate events (e.g., storm surges, heavy precipitation, and floods) at both global and regional scales [7,[9][10][11]. For this reason, it is imperative to examine potential influences of climate change on climate and hydroclimatic extremes over the Sub-Saharan African region.
To improve the understanding of these concerns in Sub-Saharan Africa, several studies have investigated the impacts of climate change on precipitation and runoff using outputs of the global climate model (GCM) or regional climate model (RCM) projections [12][13][14][15][16]. Most studies have predicted further increases in annual precipitation and enhanced spatial and temporal variations in precipitation under warmer conditions across Sub-Saharan Africa. However, the distribution of the change in annual precipitation is not the same as that of the change in seasonal precipitation. At the seasonal scale, Beyene et al. [12] suggested that the models project an increase in precipitation in the winter (DJF) season, while summer (JJA) precipitation increases or decreases according to the subregion in the Nile River Basin. Thus, changes in precipitation directly contribute to changes in runoff [11]. To understand the changes in runoff under future climate conditions, Elshamy et al. [13] assessed the annual mean changes in precipitation and evapotranspiration impact the streamflows of the upper Blue Nile based on ensemble GCMs. Liersch et al. [14] suggested an increase in the annual mean discharges as well as a shift in seasonal discharge with different change signals in the Upper Blue Nile catchment. These studies indicate that changes in temperature and precipitation will inevitably complicate water resource concerns (i.e., water availability and water-related disasters) on an annual scale as well as on a seasonal scale in Sub-Saharan Africa.
Among the many countries in Sub-Saharan Africa, Burundi is highly prone to be vulnerable to climate change. Burundi has experienced flood-related damage, including loss of life and property, due to heavy rainfall every year under the current climate [17]. Social and economic aspects of Burundi are essentially dependent on rain-fed agriculture, which sustains its economy and fast-growing population. More than 90% of the population lives in rural regions, where people are engaged in agricultural activities (the United Nations Framework Convention on Climate Change; UNFCCC) [18]. Therefore, the climate in Burundi plays the key role in controlling natural and human systems, including not only water resources, but also vegetation growth and food security [6,19]. In Burundi, to reduce catastrophic consequences induced by climate change, examining changes in hydroclimate as well as climate extremes under various climate change scenarios is necessary. However, a primary challenge for this region is that regional responses to future climate remain unknown, especially with respect to hydroclimate changes. Although previous studies on Burundi have suggested a high spatiotemporal variation in extreme events, this suggestion is based only on the observed precipitation data [7,17]. From a macro perspective, climate change affects various socioeconomic sectors of Burundi and has been surveyed solely by a research organization associated with the government agency in the form of a national report [18]. Several reports offer insight into the climate change effects in this area; however, they provide limited projections of hydroclimate changes in detail under the current global warming and likely future conditions.
The major objectives of this study are to assess the potential influence of climate change on hydroclimatic features and examine future extreme hydroclimate events (i.e., heavy precipitation and floods) considering regional seasonal climate features with a focus on Burundi (see Figures 1 and 2). Figure 1 presents the flowchart of the study procedure. To improve the reliability of climate projections, we based the results on the multi-model ensemble mean (MME) of the Earth System Model (ESM) projections involved in the sixth phase of the Coupled Model Intercomparison Project (CMIP6) [20]. The degree of agreement between the multiple results is used to assess the confidence of climate projections, which will reflect the reliability of future projections of hydroclimatic variables [21,22]. Due to the coarse resolution of these ESMs, a statistical downscaling method (i.e., quantile mapping) was used to correct the biases recorded in the ESM outputs. Next, a hydrological model (the Soil and Water Assessment Tool (SWAT) model) was set up to simulate the streamflow of the Ruvubu River in Burundi. Here, we analyzed the long-term mean changes in hydroclimatic variables (e.g., precipitation and streamflow) to explore the regional hydroclimate effects of climate change. In particular, the focus of our study was on the hydroclimatic response to extreme precipitation (and runoff) changes which can cause floods, with consideration of seasonal distinctions and future climate conditions. To the best of our knowledge, no studies have assessed the potential impact of climate change on hydrological responses in Burundi. This study can provide various interconnected sectors with scientific information to alleviate the negative impacts of climate change on Burundi.
Sustainability 2021, 13, x FOR PEER REVIEW 3 of 23 hydrological model (the Soil and Water Assessment Tool (SWAT) model) was set up to simulate the streamflow of the Ruvubu River in Burundi. Here, we analyzed the longterm mean changes in hydroclimatic variables (e.g., precipitation and streamflow) to explore the regional hydroclimate effects of climate change. In particular, the focus of our study was on the hydroclimatic response to extreme precipitation (and runoff) changes which can cause floods, with consideration of seasonal distinctions and future climate conditions. To the best of our knowledge, no studies have assessed the potential impact of climate change on hydrological responses in Burundi. This study can provide various interconnected sectors with scientific information to alleviate the negative impacts of climate change on Burundi.   Here, we analyzed the longterm mean changes in hydroclimatic variables (e.g., precipitation and streamflow) to explore the regional hydroclimate effects of climate change. In particular, the focus of our study was on the hydroclimatic response to extreme precipitation (and runoff) changes which can cause floods, with consideration of seasonal distinctions and future climate conditions. To the best of our knowledge, no studies have assessed the potential impact of climate change on hydrological responses in Burundi. This study can provide various interconnected sectors with scientific information to alleviate the negative impacts of climate change on Burundi.

Study Area
Burundi is a landlocked country in Sub-Saharan Africa confined within latitudes 2.3 • S-4.5 • S and longitudes 29 • E-31 • E ( Figure 2). The area of Burundi is 27,834 km 2 , and its elevation ranges from 739 m to 2686 m. Generally, the west, east, and northeast are lowlands, while the central, northern, and southern parts are highlands. Burundi is bordered to the west by the Democratic Republic of Congo (DRC), to the east and south by Tanzania, and to the north by Rwanda.
This study was conducted over a major part of the Ruvubu River basin located in the center of Burundi, as shown in Figure 2. The Ruvubu basin is delineated into three subbasins (subbasins 1, 2, and 3, including the main tributary of the Ruvubu River), and their characteristics are shown in Table 1. The length of the Ruvubu River is approximately 480 km, and the basin area is 9277.7 km 2 , which is approximately one-third of the total area of the country. The Ruvubu River joins the Akagera River in southern Rwanda and eventually joins the most distant southern part of the Nile River Basin. The climate in Burundi can be broken down into two major seasons within a year, wet from October to May and dry from June to September (below 50 mm/month), with the annual mean precipitation between 800 mm and 1500 mm. However, Burundi typically experiences bimodal rainfall patterns in a year since the general circulation of the atmosphere in the Intertropical Convergence Zone (ITCZ), which is one of the main factors regulating rainfall seasonality, crosses Burundi twice a year [17,18]. There is an additional short dry season in January in the middle of the wet season (from October to May), which lasts for approximately 15 days. Therefore, the dry season is divided into two seasons: the long dry (June-July-August-September, JJAS) and the short dry (January, J) seasons [18]. Additionally, the wet season is divided into two seasons: the long wet (February-March-April-May, FMAM) and the short wet (October-November-December, OND) seasons. For Burundi, four subdivisions of the year, which are more commonly used to separate seasons, are implemented in this study. The precipitation over Burundi presents predominant seasonal variation (see the details in Section 4). The primary economic sector in Burundi is rain-fed agriculture, in which 90% of its inhabitants are employed. These seasonal climate features cause excessive rain in wet seasons and drought in dry seasons, which result in crop failure and food security risks [19]. Therefore, changes in precipitation will directly influence food security and water availability in this region. In particular, heavy rain recently affected Burundi as well as parts of neighboring Rwanda, where the loss of life and damage to structures increased due to heavy rain and flooding. Figure 3 shows the spatial distributions of the long-term climatology (1985-2014) of four seasons (i.e., two wet seasons and two dry seasons) derived from 10 observation precipitation gauge data points across Burundi. One of the remarkable features of precipitation is its large variability both temporally and spatially. Individual seasonal precipitation features show somewhat greater variability across the region, even in the observations. The seasonal climate in Burundi is modified by altitude and is largely dominated by the general circulation of the atmosphere in the ITCZ, which is controlled by the convergence of trade winds [5,18]. However, differences in the mean monthly precipitation between long wet (FNAM), short wet (OND), and short dry (J) seasons are not clear in Figure 3. features show somewhat greater variability across the region, even in the observations. The seasonal climate in Burundi is modified by altitude and is largely dominated by the general circulation of the atmosphere in the ITCZ, which is controlled by the convergence of trade winds [5,18]. However, differences in the mean monthly precipitation between long wet (FNAM), short wet (OND), and short dry (J) seasons are not clear in Figure 3.

Observational Climate Datasets
In this study, observational meteorological data were used to correct the biases in the climate model outputs and simulate streamflows for the Ruvubu basin. Observational meteorological stations were selected considering the availability of climate variables (i.e., minimum and maximum temperature, precipitation, relative humidity, solar radiation, and wind speed) and long-term observation data periods . Table 2 shows the information, including the annual mean climatology (e.g., maximum temperature, TMAX; minimum temperature, TMIN; and precipitation, PRE), from all the selected ten weather stations. A 10 weather gauge stations network for the whole area of Burundi was implemented in this study. In this network, only six gauge stations (i.e., No. 2-5, 9, and 10 suggested in Table 2) are located within the Rubuvu River basin. This six weather gauge stations network is equivalent to approximately one station every 39.3 km at the horizontal spatial resolution considering the catchment area (i.e., 9277.7 km 2 ). Here is an assumption that each station is uniformly distributed in the basin. The observed weather data were obtained from the Geographic Institute of Burundi (IGEBU) at a daily timescale for the historical periods. In addition, geophysical datasets, which include DEM, soil, and land

Observational Climate Datasets
In this study, observational meteorological data were used to correct the biases in the climate model outputs and simulate streamflows for the Ruvubu basin. Observational meteorological stations were selected considering the availability of climate variables (i.e., minimum and maximum temperature, precipitation, relative humidity, solar radiation, and wind speed) and long-term observation data periods . Table 2 shows the information, including the annual mean climatology (e.g., maximum temperature, TMAX; minimum temperature, TMIN; and precipitation, PRE), from all the selected ten weather stations. A 10 weather gauge stations network for the whole area of Burundi was implemented in this study. In this network, only six gauge stations (i.e., No. 2-5, 9, and 10 suggested in Table 2) are located within the Rubuvu River basin. This six weather gauge stations network is equivalent to approximately one station every 39.3 km at the horizontal spatial resolution considering the catchment area (i.e., 9277.7 km 2 ). Here is an assumption that each station is uniformly distributed in the basin. The observed weather data were obtained from the Geographic Institute of Burundi (IGEBU) at a daily timescale for the historical periods. In addition, geophysical datasets, which include DEM, soil, and land use/vegetation cover data, were collected from various sources and used to set up the hydrological model (i.e., the SWAT model) and simulate streamflows for the Ruvubu basin (see details in Section 2.4). The DEM data from the Shuttle Radar Topography Mission (SRTM) with a 1 arcsecond spatial resolution were obtained from the United States Geological Survey (USGS). Soil data with a 5 arcmin spatial resolution were collected from the Food Agriculture Organization (FAO) [23], and land use/vegetation cover data were collected from the Regional Centre for Mapping of Resources for Development (RCMRD, http://geoportal.rcmrd.org/, accessed on 15 September 2021). The observed discharge data were obtained from the IGEBU for the calibration and validation of the hydrological model parameters. The observed discharge data from the Kibaya station on the Ruvyironza River cover a period of 6 years (2012-2017). Data quality control was carried out before the application of these data (i.e., observational meteorological data and discharge data) in the analysis. Here, the QC process consisted of identifying problems in the data (e.g., outliers and missing values) and adjusting the data value under consideration of the various data features (e.g., regional climatology and spatiotemporal consistency). Table 2. Information on the gauge stations across Burundi used in the study.

No.
Station Name

Location
Annual Climatology

Climate Change Scenarios
To obtain reliable climate projections, we used MME projections based on the CMIP6. CMIP6 climate model simulations show improved regional characteristics when based on more complicated physical processes compared with the fifth phase of the CMIP climate model [24]. In this regard, analyses of the latest CMIP6 outputs based on an ESM can provide updated results and a better understanding of climate impacts on the hydroclimate in this region.
In this study, two ESMs (i.e., KACE-1-0-G and UKESM1-0-LL) in CMIP6 were implemented to examine the impacts of climate change on the hydrology of the Ruvubu basin. The two ESMs were selected in consideration of the availability of data (e.g., climate variables, various model runs, and different shared socioeconomic pathways) at the time of analysis and model performances. These ESMs reasonably simulate observations in a wide variety of contexts [25,26]. The KMA's Advanced Community Earth system (KACE-1-0-G; K-ACE) model was developed by the National Institute of Meteorological Sciences/Korea Meteorological Administration (NIMS/KMA) in collaboration with the KMA-UK Met Office Hadley Centre. The K-ACE model is a coupled climate model (atmosphere-oceansea ice-land; AOIL), and detailed model descriptions were provided by Lee et al. [27]. Additionally, the UK Earth System Model (UKESM1-0-LL; UKESM1) implements a coupled climate model (i.e., HadGEM3-GC3.1) as a new core physical model with additional Earth system components (e.g., terrestrial carbon and nitrogen cycles) interactively coupled to each other. This model was jointly developed by the National Environmental Research Council (NERC) and the UK Met Office. Model runs were performed by the UK Met Office, the NIMS-KMA, and New Zealand's National Institute of Weather and Atmospheric Research (see Sellar et al. [25] for details), and we used ensemble model runs performed by the NIMS-KMA among various model runs. To consider differences in initial conditions for individual model runs, climate change scenarios were derived from the three initial con-dition ensemble members (i.e., r1i1p1f1, r2i1p1f1, and r3i1p1f1 for KACE-1-0-G; r13i1p1f2, r14i1p1f2, and r15i1p1f2 for UKESM1) for every two ESM outputs. Climate model outputs under the historical scenario and two shared socioeconomic pathway (SSP) scenarios (i.e., SSP1-2.6 and SSP5-8.5) were obtained from the NIMS-KMA. SSP1-2.6 is at the low end of the range of future forcing pathways and is based on sustainability-focused growth and equality, whereas SSP5-8.5 is at the high end of the range of future forcing pathways and is based on the rapid and unconstrained growth of economic output and energy use.
Since the uncertainties in input forcing (e.g., maximum/minimum temperature and precipitation data) of the hydrological model are directly transferred to the uncertainties in runoff simulations, it is important to use reliable inputs to drive a hydrological model [28]. Due to the coarse resolutions of the climate projections from the ESMs for the local watershed, the downscaling technique is generally needed to adjust the global scale to the local scale and reduce the uncertainties recorded in the ESM outputs [28,29]. To remove biases in the input forcing data at a daily timescale, statistical bias correction was carried out by applying the quantile mapping (QM) method. The QM method is a simple and effective method to adjust biases by fitting the cumulative distribution function (CDF) of raw ESM outputs to that of observations. Therefore, it is commonly used in climate change impact studies [11,[30][31][32]. The individual raw ESM outputs were downscaled to 10 observational weather stations located in the Ruvubu basin over Burundi (presented in Table 2).
After downscaling and bias correction of the ESMs, these models were validated. The validation consisted of calculating the bias for precipitation and temperature between the downscaled and bias-corrected MME mean and the observed long-term mean climatology. The MME was obtained by calculating the average mean monthly data for six ensembles of the two ESMs. In this study, to examine the future hydroclimatic response, future periods were set as two time slices of thirty years each, which are indicated as the near-term future period (2030-2059; S1) and the long-term future period (2070-2099; S2). The average variations in seasonal mean precipitation and temperature were compared to the reference period (1985-2014; S0) for two separate SSPs (SSP1-2.6 and SSP5-8.5).

SWAT Hydrological Model
The SWAT model was developed by the Agricultural Research Service of the United States Department of Agriculture (USDA). The SWAT model is a continuous, semi-distributed, and process-based model that is operated on a daily basis at the basin scale [33,34]. It is widely used to simulate hydrological processes on various scales, including small catchments and large catchments, under different climate scenarios over a long period of time and is especially suitable for African river basins [35,36].
The SWAT model for the Ruvubu basin was set up by applying five key procedures using geophysical datasets (i.e., DEM and topographic data): watershed delineation, hydrologic response unit (HRU) analysis, write input tables, SWAT simulation, and parameter sensitivity analysis (see Appendix A). The DEM was used to delineate the watershed and analyze the drainage patterns of the land surface terrain. Topographic data (i.e., soil and land use/cover maps) were used for defining the distributions of HRUs in Burundi. Each subbasin was divided into HRUs considering topographic features on the basis of the land use/cover, soil, and slope feature. In the SWAT model, the hydrological cycle is simulated on the basis of the water balance equation at a daily timescale (Equation (1)): where SW t (mm) is the final soil water volume, SW 0 (mm) is the initial soil water volume on day i, t (days) is the time, P i (mm) is the precipitation amount on day i, Q i (mm) is the amount of surface runoff on day i, E i (mm) is the evapotranspiration amount on day i, W i (mm) is the amount of water entering the unsaturated zone from the soil profile on day i, and R i (mm) is the amount of return flow on day i. Figure A1 shows a methodology flowchart for the SWAT setup in the Ruvubu basin.
Hydrological balance is simulated for each HRU based on the basic hydrological components such as evapotranspiration, soil water volume, surface runoff, lateral runoff, and groundwater runoff. In the SWAT, surface runoff and lateral runoff are simulated by the Soil Conservation Service (SCS) curve number (CN) method and kinematic storage model on a daily timescale, respectively. Individual hydrological components derived from the HRUs are summed up in each subbasin and then the components flow into the main channel of the subbasin. Channel runoff routing is estimated by the Muskingum method [33,34].
In the processes, individual hydrological components include the various model parameters. To simulate a realistic hydrological cycle, several model parameters are determined from basin features derived from geophysical data and river networks. Others are estimated for the certain HRUs in data-scarce regions from the transferred parameters estimated in the region with similar basin characteristics. After the model setup, the final step, which consisted of the calibration and validation processes, was conducted. Here, calibration and validation of the SWAT model were performed manually using observed discharge data at the Kibaya station on the Ruvyironza River, which is one of the tributaries of the Ruvubu River. The simulation results were highly influenced by the initial conditions (e.g., the initial soil conditions), especially in the data-scarce regions. Therefore, model runs were conducted during the period from 2010 to 2017 on a daily basis using observation meteorological forcing for calibration and validation processes, including two-year warmup periods (2010-2011). Since the observation streamflow data are only available for four years (i.e., 2012 and 2015-2017), we separated the data period more weighted for the calibration rather than validation. Due to limited data availability, the recent data from 2015 to 2017 corresponding to the above-normal years (2015-2016) and the normal year (2017) were used for calibration, while the data for 2012, which was a normal year, were used for calibration. The SWAT parameters were calibrated considering the sensitivity of individual parameters (see details in [28,37]), which are displayed in Table 3.

SWAT Hydrological Model Performance
To validate the simulated streamflows against the measured streamflows, the SWAT model was run based on observational meteorological input at a daily timescale (see Section 2.2). The simulated flow hydrographs from SWAT during the historical period (2012-2017) were compared with the measured streamflows during the same period except for 2 years (2013-2014) due to poor data quality. The streamflows simulated using observational meteorological input forcing are both quantitatively and qualitatively consistent with the measured streamflows of the Ruvyironza River basin (Figure 4a). The simulated streamflow captures the peaks and base flows of the measured streamflows for the calibration period as well as the validation period. Figure 4b presents box-whisker plots of daily streamflows in the Ruvyironza River basin on the basis of the measured streamflows (referred to as OBS) and simulated streamflows (referred to as SIM) forced by the observational meteorological input. SIM slightly overestimates the maximum OBS value during both the calibration and validation periods. The interquartile ranges represent the ability of SWAT to reasonably simulate the phase of the variation (i.e., lower quartile, median, and upper quartile) and high (and low) flows. To statistically evaluate the accuracy of the simulated streamflows, we considered the coefficient of determination (R 2 ), Nash-Sutcliffe efficiency (NSE), RMSE observation standard deviation ratio (RSR), and percent bias (PBIAS) as quantitative statistics. For daily data, the model performance was satisfactory: R 2 > 0.50; NSE > 0.50; RSR < 0.70; and −25% < PBIAS < +25% [38]. The values of the four statistical indicators show that the model performance is reasonable in the calibration and validation periods (see Table 4). According to the positive PBIAS value, the SWAT model slightly overestimated the measured flows due to the uncertainties in various sources, such as the hydrological model, the uncertainty in the input data, and the lack of measured flow data. Overall, the results indicate that the SWAT model is applicable to simulating the water balance components over the Ruvubu basin.

Evaluation of the Climate Simulation
Prior to examining the hydroclimate response to potential climate change based on ESM ensemble projections, historical simulations were validated to confirm whether biascorrected ESMs can suitably reproduce the fundamental climatological characteristics represented in the patterns observed in the Ruvubu basin. For the comparison during the historical period (1985-2014), the precipitation and temperature (i.e., maximum temperature and minimum temperature) data from individual gauge stations used in the SWAT model ( Figure 2) were extracted from the six ESM ensembles and observation data. Hereafter, the results obtained from the MME of the six ESM ensembles and observations are referred to as MME and OBS, respectively. Figure 5 presents the annual cycle of monthly temperature and precipitation from the MME and OBS over the Ruvubu basin. The monthly maximum and minimum temperatures in the MME are consistent with those in the OBS. The MME captures the typical shape of the annual cycle of OBS precipitation but tends to overestimate monthly precipitation. The positive bias in precipitation influences the error in the streamflow simulation. However, the MME adequately represents the seasonality of precipitation over Burundi. Regarding the magnitude of precipitation, the MME reproduces the precipitation seasonality in the OBS well, with a greater than 30-fold disparity between the wettest and driest months. The variation in precipitation between the six individual ESM outputs was determined. The inter-model spread varied in proportion to the magnitude of monthly precipitation. Therefore, the MME seems to overestimate precipitation. For the approach to seasonal precipitation, the magnitude of overestimation was higher in the months corresponding to the long dry season (JJAS) and lower in the months corresponding to the long wet season (FMAM). However, the range of precipitation in the ESM simulations was less than the interannual variability of the OBS in individual months, which was measured using the standard deviation during the historical period . Although the MME shows very good performance in capturing the monthly climatology of temperature rather than precipitation, these results demonstrate that the MME is capable of simulating the observed patterns in both temperature and precipitation.
Sustainability 2021, 13, x FOR PEER REVIEW 11 of 23 model ( Figure 2) were extracted from the six ESM ensembles and observation data. Hereafter, the results obtained from the MME of the six ESM ensembles and observations are referred to as MME and OBS, respectively. Figure 5 presents the annual cycle of monthly temperature and precipitation from the MME and OBS over the Ruvubu basin. The monthly maximum and minimum temperatures in the MME are consistent with those in the OBS. The MME captures the typical shape of the annual cycle of OBS precipitation but tends to overestimate monthly precipitation. The positive bias in precipitation influences the error in the streamflow simulation. However, the MME adequately represents the seasonality of precipitation over Burundi. Regarding the magnitude of precipitation, the MME reproduces the precipitation seasonality in the OBS well, with a greater than 30-fold disparity between the wettest and driest months. The variation in precipitation between the six individual ESM outputs was determined. The inter-model spread varied in proportion to the magnitude of monthly precipitation. Therefore, the MME seems to overestimate precipitation. For the approach to seasonal precipitation, the magnitude of overestimation was higher in the months corresponding to the long dry season (JJAS) and lower in the months corresponding to the long wet season (FMAM). However, the range of precipitation in the ESM simulations was less than the interannual variability of the OBS in individual months, which was measured using the standard deviation during the historical period . Although the MME shows very good performance in capturing the monthly climatology of temperature rather than precipitation, these results demonstrate that the MME is capable of simulating the observed patterns in both temperature and precipitation.

Future Projections of the Hydrologic Response
The different behaviors of temperature, precipitation, evapotranspiration, and runoff under warmer climates impact the annual (and seasonal) mean precipitation and runoff as well as their extremes. Table 5 shows the relative annual changes in the maximum temperature (TMAX), minimum temperature (TMIN), precipitation (PRE), actual evapotran-

Future Projections of the Hydrologic Response
The different behaviors of temperature, precipitation, evapotranspiration, and runoff under warmer climates impact the annual (and seasonal) mean precipitation and runoff as well as their extremes. Table 5 shows the relative annual changes in the maximum temperature (TMAX), minimum temperature (TMIN), precipitation (PRE), actual evapotranspiration (AET), and total runoff (ROF) during the future periods (i.e., the near-future period, S1, and the far-future period, S2) in relation to the reference period based on the MME. The temperature increase is evident across the Ruvubu basin. The TMAX (and TMIN) is projected to increase in S1 by 1.95-2.29 • C and in S2 by 2.31-5.53 • C. The temperature response is very comparable to the variation in atmospheric greenhouse gas (GHG) concentration, suggesting that the levels of warming in SSP5-8.5 are much higher than those in SSP1-2.6. PRE is projected to increase by 9.15% (26.78%) in S1 and by 7.98% (38.01%) in S2 under SSP1-2.6 (SSP5-8.5). The ROF is projected to increase by 13.65% (62.16%) in the near future and by 4.94% (65.36%) in S2 under SSP1-2.6 (SSP5-8.5). The slower rate of increase in runoff during the S2 period than during the S1 period is caused by a significant increase in the AET compared with the increase in PRE. The MME-based relative changes in both PRE and ROF are slightly higher in the S1 period than in the S2 period, especially in SSP1-2.6.

Scenarios Period Subbasin
Absolute Annual Mean Climatology  Figure 6 illustrates the seasonal changes in the TMAX, TMIN, PRE, AET, and ROF. To evaluate the consistency among climate projections, we selected the regions where 100% and 67% of the projections agreed in terms of the sign of the change. This analysis highlights the confidence of projections. For temperature changes (i.e., TMAX and TMIN), the consistent patterns that are shown in the projections under the two SSPs are an increase in both dry seasons and in both wet seasons. The signals of the changes represented in most of the areas are likely robust, with four out of the six models in agreement (i.e., over 67%). Although the patterns and magnitudes of the changes differ seasonally, significant changes in the individual variables are shown under the higher SSP scenario (i.e., SSP5-8.5). In the majority of the seasons during individual future periods under both SSPs, the TMAX, TMIN, PRE, AET, and ROF in the future are higher than those in the reference period. It is probable that the increased temperatures (i.e., TMAX and TMIN) in response to enhanced GHG emissions are unequivocal during all seasons. However, even under the same global warming conditions, the modeled seasonal mean changes in PRE, AET, and ROF are very different from those in the TMAX and TMIN. Relatively, changes in PRE, AET, and ROF seem to be less consistent than the changes in temperature. However, the seasonal variations in PRE, AET, and ROF differ according to the SSPs (i.e., SSP1-2.6 and SSP5-8.5). Similar to the precipitation and runoff projections, the SSP5-8.5 projections predicted the largest change based on the MME. Moreover, the largest increases in the results under the two SSPs mostly occur in the far-future period (S2), and their statistical significance is limited for almost all variables except the TMAX and TMIN under SSP1-2.6.  Although the annual and seasonal behavior of precipitation shows a changed pattern in future climate conditions, more fundamentally, these features result from the enhancement of the hydrological cycle. Therefore, we examined how the enhanced hydrological cycle impacts precipitation on a daily scale in Burundi based on the frequency distribution of precipitation. To this end, we analyzed changes in the frequency distribution of daily precipitation during future periods compared with the reference period (referred to as the Although the annual and seasonal behavior of precipitation shows a changed pattern in future climate conditions, more fundamentally, these features result from the enhancement of the hydrological cycle. Therefore, we examined how the enhanced hydrological cycle impacts precipitation on a daily scale in Burundi based on the frequency distribution of precipitation. To this end, we analyzed changes in the frequency distribution of daily precipitation during future periods compared with the reference period (referred to as the REF). Figure 7 shows the frequency distribution of daily precipitation at 10 stations based on the six ESM outputs for the two future periods and two different SSPs. Compared with the REF, the frequency of high-intensity precipitation increases for the future period. The relative change rate of both the frequency and amount of precipitation is individually calculated by the difference between the two future periods (S1 and S2) and the REF according to the individual intensity span. We counted all precipitation events and added up their amounts contributing to the individual intensity bin at an interval of 2 mm/day. These results imply the linkage between the change in the intensity of precipitation and the changes in total precipitation under future climate conditions. However, the ESM outputs show conflicting characteristics of precipitation according to the various intensity spans. An increase in high-intensity precipitation events occurs more frequently with the loss of a decrease in low-intensity precipitation events under warming conditions regardless of the future period and SSP scenario, as shown in Figure 8. This feature is one of the general results of intensification of the water cycle due to the anthropogenic climate. Although changes in both the amount and frequency of precipitation show relatively large variations under SSP1-2.6, they are clear and robust under SSP5-8.5. Moreover, this changing pattern becomes stronger in the S2 period under the high GHG emission scenario (i.e., SSP5-8.5) than under the low GHG emission scenario (i.e., SSP1-2.6). Thus, the precipitation change under SSP5-8.5 represents a drastic change in terms of both its amount and frequency across Burundi. These results reveal that changes in low-intensity and high-intensity precipitation may lead to increased drought and flood risk, respectively. Figure 8 shows the box-whisker plot of anomalies in the monthly maximum daily streamflow between the REF and individual future periods (i.e., S1 and S2) in the Ruvubu basin. The monthly streamflow anomalies derived from the individual periods and the six ESM outputs during a future period were calculated based on the long-term mean of the maximum daily streamflow during the REF (1985-2014). Each box-whisker represents the temporal variability of the monthly maximum daily streamflow during the period. Similar change patterns of the monthly maximum daily streamflow are shown in the months (i.e., JJAS) corresponding to the long dry season regardless of the type of SSPs. However, in general, individual 30-year anomalies in the maximum daily streamflow respond differently to the GHG emission scenarios. The box-whisker plot of the monthly maximum daily streamflow under SSP1-2.6 shows relatively little difference between the S1 and S2 periods in terms of its variability except for some months (e.g., FMAM and D) compared to SSP5-8.5. Under SSP1-2.6, anomalies of the maximum daily streamflow in the S2 period are expected to decrease in the majority of months. However, under SSP5-8.5, anomalies of the maximum daily streamflow in the S2 period are expected to increase in all the months. Overall, the mean (and above median) value of each box-whisker tends to increase under SSP5-8.5 compared to that in the REF. The highest variability in the maximum daily streamflow during both the S1 and S2 periods is shown in several months (i.e., FMAM) corresponding to the long wet season, especially in the S2 period. These results reveal that flood risk will be exacerbated in future periods due to changes in hydrological elements according to climate change, particularly under high GHG emission conditions (SSP5-8.5). Sustainability 2021, 13, x FOR PEER REVIEW 15 of 23 Figure 7. The changes in the amount and frequency of precipitation are based on the MME of six ESM outputs for the REF and both the S1 and S2 periods. Each precipitation event is accumulated over the total precipitation events corresponding to individual intensity bins. The gray error bars represent the inter-model spread of the six ESM outputs. Figure 8 shows the box-whisker plot of anomalies in the monthly maximum daily streamflow between the REF and individual future periods (i.e., S1 and S2) in the Ruvubu basin. The monthly streamflow anomalies derived from the individual periods and the six ESM outputs during a future period were calculated based on the long-term mean of the maximum daily streamflow during the REF (1985-2014). Each box-whisker represents the temporal variability of the monthly maximum daily streamflow during the period. Similar change patterns of the monthly maximum daily streamflow are shown in the months (i.e., JJAS) corresponding to the long dry season regardless of the type of SSPs. However, in general, individual 30-year anomalies in the maximum daily streamflow respond differently to the GHG emission scenarios. The box-whisker plot of the monthly maximum daily streamflow under SSP1-2.6 shows relatively little difference between the S1 and S2 periods in terms of its variability except for some months (e.g., FMAM and D) compared to SSP5-8.5. Under SSP1-2.6, anomalies of the maximum daily streamflow in the S2 period are expected to decrease in the majority of months. However, under SSP5- Figure 7. The changes in the amount and frequency of precipitation are based on the MME of six ESM outputs for the REF and both the S1 and S2 periods. Each precipitation event is accumulated over the total precipitation events corresponding to individual intensity bins. The gray error bars represent the inter-model spread of the six ESM outputs.

Future Projections of Extreme Flooding Events
We examined the future projections of extreme precipitation and extreme runoff derived from the SWAT simulations fed with the outputs of the six ensembles. Future changes were calculated based on the differences between the reference period (i.e., S0) and the two future periods (i.e., S1 and S2). Figure 9 presents the changes in extreme events derived from the annual maximum precipitation and the annual maximum streamflow for the 30 years of the individual reference period and two future periods. To investigate the extent to which the SSP1-2.6 and SSP5-8.5 scenarios affect the severity of extreme precipitation and streamflow, we compared both the annual maximum daily precipitation (MDX) and the annual maximum daily streamflow (MDF) over the Ruvubu basin. Figure 9a (and Figure 9c) shows the 30 MDXs (and MDFs) in ascending order during the individual 30-year periods according to the historical and SSP (SSP1-2.6 and SSP5-8.5) scenarios. The MDX shows considerable interannual variability, with values ranging from below 30 mm to above 105 mm. The increases in the magnitude of the MDX (and the MDF) in the S1 and S2 periods are strongly connected with the relatively high MDX (and MDF) values among sequences of the 30-year values. On the other hand, the future period seems to have slightly influenced the lower MDF values, especially in SSP1-2.6. However, the MDX (and the MDF) is projected to increase in all sequences under SSP5-8.5, and the largest increases are the maximum values. In addition, the increasing trend in the MDX (and the MDF) is more apparent under SSP5-8.5 than under SSP1-2.6 regardless of ESM outputs (Figure 9b,d). The results imply that future climate is highly likely to increase flood extremes and accompanying risks in the Ruvubu basin.
Sustainability 2021, 13, x FOR PEER REVIEW 16 of 23 8.5, anomalies of the maximum daily streamflow in the S2 period are expected to increase in all the months. Overall, the mean (and above median) value of each box-whisker tends to increase under SSP5-8.5 compared to that in the REF. The highest variability in the maximum daily streamflow during both the S1 and S2 periods is shown in several months (i.e., FMAM) corresponding to the long wet season, especially in the S2 period. These results reveal that flood risk will be exacerbated in future periods due to changes in hydrological elements according to climate change, particularly under high GHG emission conditions (SSP5-8.5).

Future Projections of Extreme Flooding Events
We examined the future projections of extreme precipitation and extreme runoff derived from the SWAT simulations fed with the outputs of the six ensembles. Future changes were calculated based on the differences between the reference period (i.e., S0) and the two future periods (i.e., S1 and S2). Figure 9 presents the changes in extreme events derived from the annual maximum precipitation and the annual maximum streamflow for the 30 years of the individual reference period and two future periods. To investigate the extent to which the SSP1-2.6 and SSP5-8.5 scenarios affect the severity of extreme precipitation and streamflow, we compared both the annual maximum daily precipitation (MDX) and the annual maximum daily streamflow (MDF) over the Ruvubu basin.  To analyze the disparity in the seasonal occurrences of the MDX and the MDF, the absolute seasonal frequency of the individual MDX and MDF values for each period (i.e., S0, S1, and S2) is represented in Figure 10. Figure 10 indicates that changes in the MDX and MDF frequency are related to their seasonal distribution. The increase in the frequency of the MDX during L.WET and the decrease in the frequency of the MDX during L.DRY and S.WET contribute to the higher occurrence of the MDX during L.WET for the S1 period. Moreover, the increase in the frequency of the MDX during S.WET and the decrease in the frequency of the MDX during L.WET and L.DRY contribute to the higher occurrence of the MDX during S.WET in the S2 period. In general, larger seasonal changes in the MDX and the MDF are shown under SSP5-8.5 than under SSP1-2.6. This result indicates that future climates cause a redistribution of the seasonal frequencies of the MDX and the MDF.
jected to increase in all sequences under SSP5-8.5, and the largest increases are the maximum values. In addition, the increasing trend in the MDX (and the MDF) is more apparent under SSP5-8.5 than under SSP1-2.6 regardless of ESM outputs (Figure 9b,d). The results imply that future climate is highly likely to increase flood extremes and accompanying risks in the Ruvubu basin. To analyze the disparity in the seasonal occurrences of the MDX and the MDF, the absolute seasonal frequency of the individual MDX and MDF values for each period (i.e., S0, S1, and S2) is represented in Figure 10. Figure 10 indicates that changes in the MDX and MDF frequency are related to their seasonal distribution. The increase in the frequency of the MDX during L.WET and the decrease in the frequency of the MDX during L.DRY and S.WET contribute to the higher occurrence of the MDX during L.WET for the S1 period. Moreover, the increase in the frequency of the MDX during S.WET and the decrease in the frequency of the MDX during L.WET and L.DRY contribute to the higher occurrence of the MDX during S.WET in the S2 period. In general, larger seasonal changes in the MDX and the MDF are shown under SSP5-8.5 than under SSP1-2.6. This result indicates that future climates cause a redistribution of the seasonal frequencies of the MDX and the MDF.

Discussion
The influence of climate change on regional features of temperature is straightforward, whereas the impacts on the features of hydrological variables are ambiguous due to strong natural variability [39]. Therefore, temporal and spatial variability in precipitation is a critical climatic factor for regulating society, water resources, agriculture, the environment, and water-related disasters in a region [40,41]. Moreover, the precipitation in Sub-Saharan Africa is particularly heterogeneous, resulting from the complexity of topog-

Discussion
The influence of climate change on regional features of temperature is straightforward, whereas the impacts on the features of hydrological variables are ambiguous due to strong natural variability [39]. Therefore, temporal and spatial variability in precipitation is a critical climatic factor for regulating society, water resources, agriculture, the environment, and water-related disasters in a region [40,41]. Moreover, the precipitation in Sub-Saharan Africa is particularly heterogeneous, resulting from the complexity of topography, water/land distribution, maritime effects, and seasonal dynamics [2]. Notably, the temporal and spatial characteristics of precipitation are heterogeneous among the countries in Sub-Saharan Africa, especially Burundi (see Figure 3). However, seasonal precipitation patterns show bimodal regimes (see Figure 5). The rainfall peaks (e.g., larger than 165 mm/month based on the MME) are generally concentrated in the months of the long wet (FMAM) and short wet seasons (OND). On the other hand, a small amount of precipitation is found in the months of the long dry season (JJAS), whereas it is not clear in those of the short dry season (J). The short dry season (J) varies in length and is related to the dry spell lengths rather than the amount of monthly precipitation. Therefore, the precipitation variability in Burundi makes it difficult to manage water resources and deal with water-related disasters in an effective way under the current climate conditions. However, robust responses of an intensified hydroclimatic regime with increased high-intensity precipitation have been suggested at the global and regional scales [39], which is in line with these results in Burundi (e.g., Figures 7 and 9a,b). The increase in extreme climate events will have adverse effects on flooding rather than drought in Burundi (e.g., Figures 8-10). In general, water-related disasters result in regional extreme precipitation. The annual maxima of precipitation, which are recognized as the criterion for determining the magnitudes of extreme precipitation events, are characterized by distinct seasonal regional climate variability at the local scale [42]. In Burundi, extreme precipitation events are impacted by different meteorological factors, i.e., a combination of local factors and enhanced air moisture from the Indian Ocean and the presence of the mesoscale convective vortex [5]. Burundi usually experiences extreme precipitation events in two wet seasons (i.e., S.WET and L.WET) during the REF period (Figure 10), which is in line with previous studies [5,17]. Climate change will increase the amount and frequency of high-intensity precipitation for future periods (Figure 7). In addition, for the far-future period (i.e., S2; 2085), the extreme streamflow will be likely to increase in most of the months (i.e., corresponding to S.WET, L.WET, and S.DRY) (Figure 8). The results indicate that further climate change is likely to increase the variety of extreme events (i.e., the MDX) at the temporal scale and the risk of floods with more intense extreme precipitation in Burundi.
Changes in the climate system originating from emissions of GHGs are expected to transform runoff and extreme runoff patterns. In particular, climate adaptation related to water-related disasters, which disrupt the economies of developing countries, is more important than other aspects to minimize damage. Therefore, the hydrological model has been generally used for simulating runoff and investigating the potential changes in runoff under possible future climate scenarios [11,43]. However, despite this foreseeable relevance and the widespread applicability of hydrological modeling, hydrological impact assessments under climate change have received little attention in Burundi. Hydrological modeling is still a major challenge in some regions, including Burundi, due to data availability (i.e., accessibility of data, lack of reliable data, and missing data) and sparsity of gauge data across time and space for observing runoff features [44][45][46]. The biggest and fundamental problem is that the observation period of the discharge data is too short for estimating the SWAT model parameters. It may contribute to the main uncertainty source of streamflow simulation. Due to these reasons related to the local data, the limitation of this study is that uncertainties are involved in simulating the streamflow in the Ruvubu River basin. Although various uncertainties are associated with simulating runoff under climate change variations (i.e., systematic biases between the ESMs and observations, variations among the ESMs and emission scenarios), Burundi is likely to face challenges in coping with extreme runoff and concurrent streamflow under global warming. Therefore, our results indicate that further research on flood planning techniques and the implementation of adaptation and mitigation measures to address climate change could be necessary for this region.

Conclusions
To our knowledge, except for some national reports, few studies have conducted a climate change impact assessment on the hydroclimate over the Burundi region. In Burundi's water resource management, one of the primary challenges is hydroclimate variability. Therefore, this study investigated the hydroclimatic response (i.e., the mean climatology and extreme events) to climate change based on the MME forced by two SSP (SSP1-2.6 and SSP5-8.5) scenarios and CMIP6 ESM simulations in the Ruvubu basin in Burundi from a seasonal perspective. After removing systematic biases, six ESM climate projections were used as meteorological forcing for the SWAT model to simulate the streamflow across the Burundi river basin. To obtain reliable streamflow simulations, the SWAT model was calibrated and validated using the observational datasets in Burundi. Future changes in hydroclimate variables (e.g., TMAX, TMIN, PRE, and ROF) were calculated by comparing the differences between individual future climate conditions (S1 and S2 periods under the two SSPs) and the reference period. The focus of this study was a comparison of hydroclimatic responses to potential future climate in terms of their seasonal means, characteristics, and extremes.
The annual mean changes in all of the hydroclimate variables (e.g., TMAX, TMIN, PRE, AET, and ROF) except for TMIN during the S1 period are much larger for the two future periods (S1 and S2) under SSP5-8.5 than under SSP1-2.6. The annual mean changes in both PRE and ROF during the S2 period show lower increases than during the S1 period under SSP1-2.6. Meanwhile, changes in other hydroclimate variables show intensified hydrological response in a future period (i.e., S2 period; 2085). From a seasonal point of view, the seasonal means of the TMAX (and TMIN) will similarly increase by 1.7-2.4 • C (1.7-2.2 • C) in all the four seasons (i.e., S.DRY, L.WET, L.DRY, and S.WET) during the S1 period regardless of the SSP scenario (see Table 5). However, a larger difference (i.e., approximately 2.9-3.6 • C) between the two SSPs is shown for both TMAX and TMIN during the S2 period due to the impact of high cumulative GHGs by the end of the twentyfirst century. The potential changes in the TMAX and TMIN are more obvious and robust (with six out of six ensemble models in agreement) than those in other hydroclimate variables. In line with the expected increases in the TMAX and TMIN, the seasonal mean changes in the AET during the S1 period (i.e., approximately -2.8-6.8%) compared with the S2 period (i.e., approximately 7.7-54.9%) slightly differ between SSP1-2.6 and SSP5-8.5. However, individual seasonal changes in PRE and ROF are quite different in the two future periods and two SSP scenarios. The changes in PRE (and ROF) are dominated by various factors, including not only the global temperature rise but also the localized natural variability. Regarding the precipitation extremes, adaptation to enhanced heavy precipitation (measured by the MDX) in terms of both magnitude and seasonal occurrence will be necessary for Burundi. The magnitudes of the MDX and the MDF are projected to increase more in future periods, especially during the S2 period under SSP5-8.5, than during the reference period.
The results of this investigation reveal that future climate change will impact the mean climatology (i.e., the annual mean and the seasonal mean) as well as the extreme values of hydroclimate variables (e.g., TMAX, TMIN, PRE, AET, and ROF) with a certain level of reliability in the Burundi region. In addition, changes in precipitation extremes increase the risks related to flood extremes and damage socioeconomic sectors. Moreover, the daily maximum runoff, which is likely to increase the flood risk, increases in the future period irrespective of the type of climate projection. As suggested in this study, the climate of Burundi is highly influenced by climate change, similar to other countries in Sub-Saharan Africa. Therefore, understanding the long-lasting impacts of global warming on the mean changes in hydroclimate is critical to the planning of water resource management and increasing agricultural production in Burundi. Structural and nonstructural measures for water-related disasters are needed to adapt to climate change, especially to minimize the potential flood risk in Burundi.