Integrating a GIS-based approach and a SWAT model to identify potential suitable sites for rainwater harvesting in Rwanda

According to the UNESCO World Water Assessment Programme (2019), since the 1980s, global water consumption has increased about 1% per year, due to a combination of socio-economic development, population growth, and the change in the consumption patterns. Consequently, freshwater scarcity is becoming a serious issue in several parts of the world, especially in developing countries (Gasirabo et al. 2019). Climate change and its related rainfall variability led to variations in water resources as well as extreme events such as droughts and floods, which exacerbate water scarcity and threaten sustainable development (Kaini et al. 2020b; Ayugi et al. 2021). Similarly, changes in the hydrological regime are anticipated owing to the reported climate change impacts on the spatial and temporal variations in precipitation patterns (Oo et al. 2020; Kaini et al. 2021), which have forced mankind to find alternative sources of water supply. This calls for a springboard for the design of future adaptation strategies to climate change. Among different alternatives, rainwater harvesting (RWH) is the most promising and has been listed by the Intergovernmental Panel on Climate Change (IPCC 2014) as a structural approach for intended adaption.

Africa is facing serious water constraints (Loucks & Jia 2012). The rapid population growth is pushing up water demand and hastening water resource degradation. By 2025, around 230 million people will live in African countries where water is scarce (UNESCO 2009). Therefore, concentrating more efforts on RWH at the catchments scale is one of the viable approaches for improving water resource management. Aside from its significant contribution to boosting groundwater supplies, the RWH provides several other benefits, notably the absence of excessive runoff, flood control in the watershed, greater availability of water from the soil, and soil conservation (Karamage et al. 2016).

Rwanda has one of Africa's lowest per capita water supply and storage capabilities (670 m³ per capita per year) (Mahreen & Shivali 2019). The safeguard of the available water and its sustainable management have been designated as a priority for the 2030 vision as goal 6 of the Sustainable Development Goals (SDGs). This states that optimizing water demand management activities by solving the challenges of water resources is crucial (Mahreen & Shivali 2019). Although Rwanda receives a relatively high annual average rainfall (Uwizeyimana et al. 2019), the full potential of rainwater is not exploited. Features such as steep slopes and uncontrolled terrain speed up surface runoff and so much of the rainwater go unevenly. Consequently, the country is often perceived to be susceptible to increased erosion, flooding, and drought especially in the Nyabugogo catchment area (Mind'je et al. 2019b). Owing to the above, the RWH system is the best technique, which can be effectively used to trap unused surface runoff and thereby increase groundwater recharge. However, RWH structures have to be located in places where water is available in large quantities with conditions favorable for increased infiltration (Kumar & Jhariya 2017).

Various techniques and approaches have been established in selecting suitable sites for RWH. Many of these produced thematic maps derived from the geographic information system (GIS) integrated with remote sensing (RS) data in the identification of suitable RWH sites using multi-criteria decision-making (MCDM) (Karimi & Zeinivand 2019). The identification of suitable sites for RWH structures relies on numerous factors as suggested by the Food and Agriculture Organization (FAO), such as climate (rainfall), hydrology (rainfall–runoff), topography (land slope), land cover land use, soil properties, and different socio-economic conditions (distance to settlement/streams/roads) (Adham et al. 2018). In the same perspective, several hydrological models are available to estimate rainfall–runoff, of which the Soil and Water Assessment Tool (SWAT) was applied in this research owing to its availability, convenience, friendly interface, and simple operation. The SWAT offers components and features that are useful to simulate water balance, sediment depletion, climatic change, agricultural development, and land management methods (Guiamel & Lee 2020). High surface runoff often causes soil erosion, which results in sludge deposits at the bottom of reservoirs, lowering water quality (Liu et al. 2020). Therefore, soil erosion becomes an important factor to consider in the improvement of the rationality and accuracy of results.

From an extensive literature review, it was clear that there is a distinct lack of research information on the technical and managerial aspects of RWH technologies, particularly in developing countries including Rwanda. Furthermore, the paucity of linkages between research and learning institutions, the corporate sector, and key organizations are also significant barriers to the establishment of new RWH structures. Moreover, due to several challenges, including the spatial and temporal scarcities of meteorological data, biophysical data, and infrastructure, this topic has received less consideration in Rwanda. However, integrating geospatial analysis and modeling approaches have become a promising tool to identify suitable sites for RWH. Previous studies in the study area have been conducted with full attention on assessing flood mitigation, soil erosion, and landslide (Hakizimana et al. 2021; Icyimpaye et al. 2021). A very recent study assessed the hydrological response to rainfall events through discharged flow and volume simulation (Mind'je et al. 2021), while none has considered the assessment of rainwater retention mechanisms which is one of the best solutions for soil erosion and flood control (Karamage et al. 2016). Therefore, this study is conducted to bridge the gap identified in the literature related to the provision of scientific-based information on the potential sites for RWH structures in Rwanda. This will positively influence the scientifically based decision-making on the effective management and conservation of rainwater. A GIS-based conceptual framework is used in line with MCDM using the Analytic Hierarchy Process (AHP) and the SWAT. Thus, this study aims at (1) estimating the surface runoff potential and sediment yield; (2) assessing the potential sites for RWH; and (3) proposing the appropriate sites for the RWH structures in the study area.

Nyabugogo is a level 2 catchment within the Nile basin and a tributary of the lower Nyabarongo River spreading over the central, eastern, and northern parts of Rwanda. The catchment (Figure 1) has a total surface area of 1,661 km², representing 6.31% of Rwanda's total area. It has rural areas but also a mostly urbanized and densely populated catchment in Rwanda. The western part of the study area has more shale contents, whereas the middle and east had mostly schist and quartzite alteration with significant granite and pegmatite. Throughout the catchment, alluvial material is found at the valley's bottom. Nitosol, acricol, alisol, and lixisol are the most common soil types in the catchment, with ferralsols in the eastern part near Lake Muhazi and acricol in the western part. The camisole is found in abundance in the western part. Clay soils with low infiltration and a flat topography define the catchment's central part and valley bottom (Manyifika 2015).

Similar to the entire country, the Nyabugogo catchment enjoys a tropical temperate climate with mean annual precipitation varying between 992 and 1,128 mm and temperatures varying between 19 and 21 °C. It has two rainy seasons – one in middle September to December and a peak rainy season during March up to early May. The Nyabugogo River traverses the City of Kigali and has main tributaries such as Mwange, Muyanza, Rusine, Kajevuba, and Yanze rivers. It is 45.97 km long, measured from the outflow of Lake Muhazi to its confluence with the lower Nyabarongo River near Kigali at an altitude of about 1,360 m.a.s.l. The highest point of the catchment is 2,281 m.a.s.l in the northern part. The central feature of the catchment is Lake Muhazi, which is about 80 km long in the east–west (Munyaneza et al. 2013). Agriculture, fishing, and forestry (24%) are the most important employment sectors in the area, followed by trade (20%) and other services (21%).

In this study, all necessary thematic layers were considered using RS and GIS to identify suitable sites for RWH. Using ArcGIS 10.7, basic thematic layers such as the soil properties, land cover/land use (LCLU), elevation, slope, and distance to rivers/roads/settlements (Table 1) were prepared. ArcSWAT, an ArcGIS extension, and interface for the SWAT were applied to estimate the surface runoff potential and sediment yield.

The LCLU map (Figure 3(a)) was produced from the Landsat-8 Operational Land Imager (OLI) image (path/row: 172/61) obtained from the United States Geological Survey (USGS) through a global visualization tool. LCLU is the main component of runoff and evapotranspiration control, as it is a significant factor for the detection of RWH structure (Kadam et al. 2012). The selection of Landsat imagery was influenced by its quality particularly for those with limited or low cloud cover. The acquired image was then used to perform classification using the maximum likelihood classification (MLC) technique in ENVI software version 5.5.2 (Mind'je et al. 2019a). The main emphasis for accuracy assessment pixel samples was on areas that could be visibly identified on Google Earth imagery. Thus, a total of 80 random points in all land-use types were sampled and 72 total ground true value data were found. The classification accuracy of the LULC map exhibited an overall accuracy of 90%, which was very satisfactory for the purpose of this study.

The topographical features of an area play a vital role in the generation of surface runoff and the implementation of the RWH structure as they influence the amount of sedimentation and the speed of water flow (Adham et al. 2018). For this, the Shuttle Radar Topographic Mission (SRTM) provided the digital elevation model (DEM), 30 m resolution data of the USGS Earth Explorer community used to produce the slope (Figure 3(b)), the distance to roads (Figure 3(c)) using road networks attained from the Rwandan Transport Development Agency (RTDA), distance to rivers (Figure 3(f)) through river networks collected from the DIVA-GIS, and distance to settlement (Figure 3(d)) using the Euclidean distance tool of the Spatial Analyst tool. The latter representing the socio-economic factors have been used as they play a significant role in the prevention of any future conflict between the constructed RWH structures and the development of the infrastructure (Al-Adamat et al. 2010). The soil properties (Figure 3(e)) control the infiltration characteristic of the soil, making it one of the most significant factors in RWH planning (Tumbo et al. 2013; Harka et al. 2020). Data for soil texture were derived from the Digital Soil Map of the World (DSMW) (http://www.fao.org/geonetwork/srv/en/metadata.show%3Fid=14116) assembled by the United Nations FAO.

For the runoff volumes prediction, the SWAT uses a revised computational efficient version of the Soil Conservation Service Curve Number (SCS-CN) method (USDA 1972), relating runoff to land use, management practices, and soil type. In addition, the SWAT model uses the Modified Universal Soil Loss Equations (MUSLE) to calculate HRU-level soil erosion.

The input data needed by the ArcSWAT include the DEM, soil map, LCLU, and weather data. The Nyabugogo catchment divided into 35 sub-basins was delineated using the DEM and provided topographic parameters such as overland slope and slope length for each sub-basin. The final LCLU map (Figure 3(a)) was divided into six classes with cropland (64%) and forest (17.13%) as dominant classes. According to the slope map (Figure 3(b)), the study area was classified into five classes by the predominance of slope classes between 15 and 30%, which extend over approximately 32.29% of the total area of the basin. In the part where the grades of the slopes are strong to an extreme, they extend on 32.48%, located mainly in the center and the north-west, whereas the low slopes represent only 3.34% of the basin. The soil map (Figure 3(e)) was classified into Fo42-3b-497 (17.51%), Fo96-3b-558 (0.40%), Fo-97-3b-559 (32.08%), Nh7-2-3c-853 (46.59%), and WATER-1972 (3.42%). To discretize the basin into 429 HRUs, we used the option of various HRUs, each representing a distinct combination of land cover and soil type. Data on the weather (daily rainfall and minimum and maximum temperatures) were used for the prediction of streamflow and sediment yield.

The model calibration was accomplished with multiple runs using both manual and auto-calibration. Initially to the calibration and validation steps, the identification of the most sensitive parameters for the study area was carried out and then was used to calibrate the model. Table 2 summarizes the sensitive parameters used in calibration, whereby numbers 1–10 were used to calibrate the simulated flow and numbers 11–14 were used in sediment yield calibration. To accomplish the uncertainty analysis, the flow output was calibrated and validated using the Sequential Uncertainty Fitting, ver. 2 (SUFI-2) algorithm (Zettam et al. 2017). This program is tied to the SWAT and is included in the SWAT Calibration Uncertainty Procedures (SWAT-CUP) calibration software application (Abbaspour et al. 2004). A total of 2,500 runs of calibration, divided into five iterations of 500 runs each, were performed for each simulation.

In multi-criteria modeling studies, multicollinearity analysis is crucial for optimizing and carefully selecting the relevant influencing factors that will suit and match the modeling process. This assists in the selection of the suitable influencing factors to use by excluding those that are not critical for optimal modeling (Dormann et al. 2013; Isner 2014). A multicollinearity test was performed to exclude the most highly correlated factors that could lead to errors and reduce the accuracy of the modeling results. Multicollinearity between continuous variables was assessed using the most preferred method, the variance inflation factor (VIF) calculation (Zheng et al. 2021). The VIF value greater than 10 signifies serious multicollinearity (Rahimian Boogar et al. 2019). The maximum VIF value was evaluated at 1.33. Therefore, in this research, there was no multicollinearity found among the effective independent factors.

In considering sites for RWH, not all factors are regarded as equal. In this study, the MCDM (Rana & Suryanarayana 2020) process was used while the weighted overlay analysis approach was applied based on the criteria following the FAO guidelines. The suitability levels were assigned a value on a scale of 1 (not suitable) to 5 (highly suitable). This study proposed suitable sites for constructing check dams, bench terraces, and farm ponds. Numerous factors (Table 3) were considered to determine the spatial site of RWH structures within the study area.

To estimate the amount of flow from the catchment, sensitivity analysis was done to identify the parameters most relevant in affecting the streamflow simulation (Mind'je et al. 2021). Table 2 shows the detected parameters to be important in regulating the streamflow generation for all the simulations. In the study area, number 1 was ‘very sensitive’, from numbers 2–5 ‘sensitive,’ while numbers 6–10 were slightly sensitive. From the SWAT initial estimations, these parameters were related to runoff and adjusted to suit the model simulations with the observed flow data. In the SWAT, these parameters are usually used to calibrate the base flow, which was confirmed by Harka et al. (2020). The entire simulation was performed daily from January 2005 to December 2017 (excluding a 3-year warm-up from 2005 to 2008). Generally, as confirmed by the previous literature related to model validations (Chaibou Begou et al. 2016), there are no rules of thumb for selecting the proportion of calibration and validation datasets depending on the available amount of data (years). Due to limited observed discharge data similar to this study, Mourad et al. (2005) suggested using the larger portion of the data for calibration and the smaller portion of the data for validation. Chaibou Begou et al. (2016) also recommended using the split-sample approach, by splitting the discharge measurement data into two datasets: two-thirds for calibration and the other one for validation. Therefore, the available observed data in the study area were split from January 2011 to December 2012 for calibration and January 2013 to December 2013 for validation (Figure 4). The statistical performance of the model was considered satisfactory (Moriasi et al. 2007) with 0.9 of R² and 0.73 of the NSE for calibration, and 0.86 of R² and 0.7 of the NSE for validation.

The balance between rainfall and water loss due to evapotranspiration, groundwater recharge, and runoff is the water budget for a watershed (Xu et al. 2015). The runoff potential of the area (Figure 5) affects the recharge and movement of surface water and is an important parameter for RWH. The runoff potential was classified into five classes from very high to very low runoff potential: very low (<30 mm), low (30–60 mm), moderate (60–90 mm), high (90–120 mm), and very high (>120 mm). A large portion of the study area (55.3%) fell under the high runoff potential class with an area extent of 903.27 km², while a very small part (1.8%) fell under the very low runoff potential class. The settlement areas increased the impervious layers, thus reducing the infiltration rate, and this leads to an increase in runoff (Maina & Raude 2016; Harka et al. 2020). Settlement areas are the most vulnerable in the catchment, where the retention capacity is low. The average annual runoff potential for cropland, settlement, and grassland was 98.73, 139.19, and 56.63 mm, respectively. Very high to moderate runoff potential areas were suitable for harvesting rainwater. Previous studies (Ibrahim et al. 2019; Sayl et al. 2020) have used runoff potential to identify suitable RWH sites, with the conclusion that runoff can be harvested using RWH structures such as check dams and ponds (Chowdhury & Paul 2021).

As argued by White & Arnold (2009) and Ouyang et al. (2010), there can be many impacts on the movement of soil from one location to another, such as the openness of the ground and the erosion of plant nutrients. Deposition of eroded soils can cause waterway disruptions (e.g. water deposition in rivers and reservoirs) and may reduce the quality of water. Soil erosion is nearly always associated with rainfall, soil properties, topography, and land use. The kinetic energy of soil erosion and the probability of erosion increase as rain and runoff intensity increase. Previous studies (Prinz & Singh 2000; Zheng et al. 2018) have suggested that when identifying RWH sites, considering sediment yield can improve the rationality and precision of the results. The annual average rate of sediment delivered by the Nyabugogo catchment is approximately 40.8 tons/ha/year, which resulted in an annual average soil loss of 294,465.9 tons. This rate is comparable to that estimated using the RUSLE (51.3 tons/ha/year). To determine the potential severity of soil loss, the estimated average annual loss of soil from the basin has been classified into five erosion classes. The results obtained with the RUSLE (Figure 6(a)) application show that 74.6% of the study area was exposed to a moderate to very low risk of erosion and 22.2% had a high to very high risk of erosion. The soil erosion results obtained using the SWAT (Figure 6(b)) indicate that 70.6% of the catchment area was exposed to a moderate to very low risk of erosion, and 27.4% had a high to very high risk of erosion. Earlier studies (Salem et al. 2014; Ibrahim et al. 2019) echoed severe soil erosion happening in agricultural areas, mostly under higher slopes. This is in accordance with the results of this current study where a higher soil erosion level is attributed to the high percentage of cropland (about 88.7%) with a mountainous slope (>30%) in sub-basin 1. Low to moderate erosion risk areas were considered to be suitable for RWH as revealed by an earlier study (Zheng et al. 2018). Previous researches (Polyakov et al. 2014; Mekonnen et al. 2015) disclosed that RWH techniques such as stone bunds, check dams, ponds, and terraces can lessen catchment sediment yield. Hence, similar techniques can be used in the study area for the effectiveness of sediment yield reduction and to meet the amount of water demanded by the population to some extent. However, the combination of some of these techniques can decrease sediment yield from the catchment, despite the reported benefits of individual RWH techniques by earlier studies (Mekonnen et al. 2015; Grum et al. 2016).

The potential RWH area (Figure 7) was divided into classes ranging from 1 to 5, where 1 represents the least suitable areas and 5 represents the most suitable areas (Table 4). The regions that were considered not suitable covered 1.3% of the total area, while the most suitable and suitable regions covered 4.8 and 16.35% of the entire area, respectively. The moderately suitable areas covered 36.5% and the marginally suitable areas covered a huge portion of the study area at 39.9%. Results showed that some of the southern and eastern parts of the catchment were most suitable for RWH. These areas were deemed suitable as they are located in areas with milder slopes and several tributaries, which is consistent with previous studies (Maina & Raude 2016; Adham et al. 2018; Harka et al. 2020). Most of the areas with very high to high suitability ranges texturally are clay loam, surface runoff potential (more than 120 mm), and soil erosion within low to very low risk. Among the key parameters, surface runoff, sediment yield, and slope are essential factors in identifying the suitability of RWH areas. Zheng et al. (2018) showed that areas with high runoff and low erosion risk were highly recommended for RWH sites.

To optimize water resources, previous literature (Kadam et al. 2012) unveiled important criteria (Table 3) in selecting zones suitable for bench terraces, farm ponds, and check dams (Figure 8). Farm pond structures are suggested for agriculture areas since supplemental irrigation for crops is needed. During the construction of new irrigation structures and the repair of existing irrigation canal systems, climate change influences on irrigation water requirements should be considered (Kaini et al. 2011). Areas suitable for farm ponds were mostly located in the eastern part of the catchment with 6% (98.5 km²) of its total area. Farm ponds are located in places with gently undulated slopes (<8%) and with clay loam soils in agreement with an earlier study (Srivalli & Singh 2019). The possibility of diverse uses of water should be investigated for the long-term sustainability of water supplies and irrigation systems (Kaini 2016). The soil texture is a key criterion for determining a suitable area for bench terracing. Since soils with a high degree of aggregation and stability resist soil erosion better than those with a low degree of aggregation and stability, the higher the clay content of the soils, the better the chances of finding bench terraces (Mugo & Odera 2019). Bench terraces cover 25.9% (423 km²) wherein the most suitable areas for bench terraces are in the northern part of the watershed, the bench terraces are mostly located on slopes ranging between 15 and 30%, which is in agreement with findings obtained by Mbilinyi et al. (2007). Bench terraces were found at distances ranging from 0 to 350 m along streams, but no bench terraces were found at distances greater than 350 m. These results agree with the findings by Tumbo et al. (2013). Moreover, check dams are more important as they assist in water conservation and soil erosion. For lower-order streams, check dams are generally recommended, and the terrain slope should be flat and gentle to retain the maximum amount of water. Check dams were found to be 26.1 km² (1.6% of the total study area) mostly located on steep slopes below 15% in line with previous studies (Kadam et al. 2012). The study highlights that suitable areas for RWH structures account for approximately 33.6% of the total area.

Subsequently, the drainage and road network maps were combined with the RWH zone map to identify 50 potential recharge sites for farm ponds, bench terraces, and check dams. In general, these RWH structures in the study area can help with surface storage and aid in soil erosion and flood control (Karamage et al. 2016). Additionally, a settlement buffer map was overlaid to ensure that no RWH structures were located within 500 m of the settlement (Figure 9). Of the 50 recharge sites, only three sites are located in the southern parts of the study area. About 29 sites are located in the most suitable area for RWH. The proposed RWH structures provide insights into future development, which is beneficial for policymakers and development planners. However, socio-economic factors at the local level influence the extraction and use of water resources (Kaini et al. 2020a; Suhardiman et al. 2020).

This study attempted to identify possible sites for RWH structures across the Nyabugogo catchment. Here, the biophysical factors, including the slope, LCLU, soil, runoff potential, sediment yield, and socio-economic factors (distance to rivers/roads/settlements), were considered for site identification. Therefore, the study has proved the ability and effectiveness of relying on RS and GIS integrated with both SWAT and MCDM for identifying potential sites for RWH. This could be substantial for the development and management of RWH programs in Rwanda. The results revealed that 4.8% of the study was highly suitable for RWH, suitable (16.3%), moderately suitable (36.5%), marginally suitable (39.9%), and not suitable (0.3%). The zones suitable for constructing bench terraces, farm ponds, and check dams covered 423 km² (25.9%), 98.5 km² (6%), and 26.1 km² (1.6%), respectively. A total of 50 sites for the RWH structures have been identified.

Regardless of the scarcity of data, the current research gave an insight into how integrating MCDM in RS and GIS integrated with the SWAT model is an effective solution for assessing runoff potential, sediment yield, and RWH potential site identification. Although this study identified and proposed RWH structures to combat the increasing water scarcity and aid in flood control in the study area, further research is suggested to consider more factors such as heavy metals and organic pollution in the water to meet different water use purposes, including irrigation and domestic use, among others. Finally, appropriate and extensive fieldwork should be conducted on the chosen sites to ensure that they do not conflict with other land uses in the region that are not spatially identified by the available data.

The authors appreciate the United States Geological Survey and FAO for making data accessible at no cost. We also express sincere gratitude to the Rwanda Natural Resources Authority and Rwanda Meteorological Agency for the provision of data. As well, the authors greatly acknowledge the supports received from the Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences (CAS).

This research has been supported by the Strategic Priority Research Program of Chinese Academy of Sciences (no. XDA20060303), the CAS International partnership Project (no. 902020E01010), the CAS Interdisciplinary Innovation Team (no. JCTD-2019-20), the Regional Collaborative Innovation Project of Xinjiang Uygur Autonomous Regions (no. 2020E01010), the National Natural Science Foundation of China (no. 41761144079), and the project of the Research Center of Ecology and Environment in Central Asia (no. Y934031). A special acknowledgment is expressed to K.C. Wong Education Foundation and China-Pakistan Joint Research Center on Earth Sciences that supported the implementation of this study.

All authors declare no conflict of interest

All relevant data are included in the paper or its Supplementary Information.