The Nzoia basin is situated in the tropical zone at the northeastern edge of Lake Victoria in western Kenya, as shown in Figure 1. It covers an area of about 12,700 km² with altitudes ranging from 1,140 m.a.s.l in the lower plains near Lake Victoria to 4,300 m.a.s.l in the northern mountainous areas. The precipitation regime in the region is trimodal, with precipitation peaks in April-May, in August, and in October, and with annual precipitations varying from around 1,000 mm in the lowlands to above 2,000 mm in the highlands (Githui et al., 2010). More than half of the Nzoia basin is covered by agricultural areas, which leads to over-exploitation of resources and land degradation to fulfill the growing food demand of increasing population (Githui et al., 2010).

The climatic characteristics in the Nzoia basin are evaluated with gridded observational products of precipitation, temperature, energy fluxes, and discharge data. The Integrated Multi-satellitE Retrievals for Global precipitation measurement dataset (IMERG; Huffman et al., 2014) is a near-global and daily product, provided on a regular grid with a 0.1° horizontal resolution and available for the period 2001–2020. The Climate Research Unit temperature dataset (CRU; Harris et al., 2014) is a global and monthly product provided on a regular grid with a 0.5° horizontal resolution and available for the period 2001–2020. The energy fluxes dataset from the machine-learning based initiative to upscale biosphere-atmosphere fluxes from FLUXNET sites to continental and global scales (FLUXCOM; Jung et al., 2019) is a global and monthly product provided on a regular grid with a 0.5° horizontal resolution and available for the period 2001–2013 for FLUXCOM. The discharge data comes from the gauge station 1EEF01 located at the basin outlet (0.124°N, 34.090°E), as displayed in Figure 1B, and was provided by the Water Resource Authority of Kenya as a daily timeseries for the period 2009–2017.

The regional climate model WRF-Hydro selected for this study is based on version 4.4 of the Weather Research and Forecasting model (WRF; Skamarock et al., 2019) coupled with the version 5.2 of the WRF-Hydro hydrological module (Gochis et al., 2020). WRF-Hydro resolves the equations of atmospheric motion on a three-dimensional grid over a limited area, with a set of parameterization options to represent subgrid scale physical processes, such as radiation, turbulence, cumulus convection, cloud microphysics, and terrestrial hydrology. The two model domains shown in Figure 1 are selected for this study. The outer domain covers most of East Africa with an area of 4,000 km × 6,000 km at 50 km resolution and 50 vertical levels up to 10 hPa. The initial condition and lateral boundaries of the outer domain are constrained by the air pressure, geopotential, zonal and meridional winds, temperature and water vapor provided by the ERA5 reanalysis (Hersbach et al., 2020) on pressure levels at a six-hourly interval, with a 0.25° resolution. The inner domain, which encompasses the study region with an area of 800 km × 800 km at 10 km resolution and the 50 vertical levels, is driven by the outer domain using a one-way nesting method. The equations of atmospheric motions in the outer and inner domains are resolved at 180 s and 60 s timesteps, respectively, to ensure numerical stability.

The outer and inner domains share the following physics parameterization setup: the long-wave and short-wave radiation schemes of Dudhia (1989) and Mlawer et al. (1997), the Mellor-Yamada-Nakanishi-Niino Level-2.5 version turbulence scheme of Nakanishi and Niino (2004), the cumulus convection scheme of Grell and Freitas (2014), the six-class microphysics scheme of Hong and Lim (2006), and the community Noah land surface model with multi-parameterization options of Niu et al. (2011) referred to as Noah-MP. This choice of physics parameterization is motivated by its relatively good performance in terms of basin-average daily precipitation for the Nzoia basin during the period chosen for discharge calibration, with a correlation coefficient of 0.74 and a percentage bias of 7 %, as displayed in Figure 2A. Additionally, the inner domain includes lateral terrestrial flow with the WRF-Hydro hydrological module (Gochis et al., 2020).

Noah-MP is the land surface module in WRF-Hydro and describes the fate of snow cover, vegetation canopy, and soil moisture within a soil column of 2-m depth divided into four layers, to update the upward energy and water vapor fluxes at the lower boundary of the modeled atmosphere. The upward water vapor flux provided by Noah-MP is divided into five components that are the snow sublimation, the intercepted canopy water evaporation, the plant transpiration, the soil water evaporation below the canopy, and the soil water evaporation outside the canopy. Moreover, the plant stomatal resistance is calculated according to the Ball-Berry option (Eq. B1 in Niu et al., 2011) and the leaf area index is derived from a leaf dynamics model (Eq. 10 in Niu et al., 2011). The parameters used in these equations, such as, among others, root zone depth, carboxylation rate, leaf reflectance, leaf transmittance, leaf turnover, leaf stress death coefficient and leaf maintenance respiration, are prescribed as a function of landcover classes from the moderate resolution imaging spectroradiometer (MODIS) landcover map (Friedl et al., 2002). The soil hydraulic parameters are prescribed as a function of soil classes from the State Soil Geographic/Food and Agriculture Organization soil database (FAO, 1991). The initial condition of the land surface state is deduced from the ERA5 reanalysis.

The WRF-Hydro routing modules are activated for the inner domain through a coupling between the inner domain grid at 10 km resolution and the subgrid at 1 km resolution displayed in Figure 1B, generated with the WRF-Hydro Pre-processing Tool and the digital elevation data from the hydrological data and maps based on Shuttle Elevation Derivatives at Multiple Scales (HydroSHEDS) database (Lehner et al., 2008). At each time step, the surface water and soil moisture variables from Noah-MP are disaggregated on the subgrid, routed overland, in the subsurface, and in the river channels according to diffusive wave formulations detailed in Gochis et al. (2020), and aggregated back to the Noah-MP grid. The disaggregation between the coarse grid and fine subgrid is obtained with a so-called disaggregation factor, which is defined for each variable as the ratio between the variable value on the fine subgrid and its corresponding value on the coarse grid and is updated at the end of each time step (Gochis et al., 2020). The aggregation is done as a linear averaging, so that the soil moisture field on the coarse grid is updated with the averaged effect of terrestrial water transport resolved on the fine subgrid. This fully coupled formulation allows for a comprehensive description of the water cycle at basin scale.

As in the WRF-Hydro setup of Arnault et al. (2016) for a West African Sahelian river basin, no groundwater bucket option is considered, because the baseflow produced with this approach would bring more water into the stream compared to what is observed. The groundwater bucket option being disabled, the water percolating below the 2-m soil column is considered as pure groundwater recharge, which can be calibrated with the so-called SLOPE parameter (e.g., Rummler et al., 2019) referred to as the percolation parameter S.

A WRF-Hydro simulation using the setup detailed in section 2.2 with default parameters produces too high daily discharge peaks according to gauge observation, as displayed in Figure 2B, with a Nash-Sutcliffe efficiency (NSE, Nash and Sutcliffe, 1970) of −2.68 and a percentage bias of −5%. Furthermore, the ratio of total discharge to the total volume of precipitation is 23% according to the observational data, and 20% according to the simulation. This shows that the total amount of discharge obtained with the default WRF-Hydro setup is approximately well calibrated, although the surface runoff generation occurs too quickly.

The usual approach to improve the WRF-Hydro simulated discharge is to calibrate the most sensitive parameters, such as the percolation parameter, the runoff-infiltration parameter, the lateral hydraulic conductivity factor, the overland roughness factor, the retention depth factor, and the river roughness Manning coefficients (Senatore et al., 2015; Rummler et al., 2019; Zhang et al., 2019; Camera et al., 2020; Fersch et al., 2020; Gochis et al., 2020; Li et al., 2020). However, the preliminary sensitivity experiments summarized in Table 1 show that the tuning of these parameters keeps on producing the high discharge peaks characterized by a negative NSE. These unrealistically high discharge peaks can be only minorly improved at the expense of larger underestimation of the total discharge. This suggests a process is missing in WRF-Hydro to slow down the simulated runoff generation in the case of the Nzoia basin.

As mentioned by Onencan et al. (2016), the Nzoia basin is prone to floods as the river flow pushes over its banks. Such an overbank flow is not considered in WRF-Hydro, which could be a reason for the much smoother observed discharge peaks compared to the simulated ones in Figure 2B. In order to circumvent this issue, the WRF-Hydro source code has been updated with an overbank flow option to allow for water flow between land and river in a two-way manner. In default WRF-Hydro, the surface water enters a channel if it exceeds the retention depth and never comes back to the land (Gochis et al., 2020). With the overbank flow option, a new parameter, referred to as the overbank flow parameter H_thres in meters, is introduced. When the water head in a channel pixel exceeds H_thres, then it does not receive the water flow originating from the upstream channel pixel, and this upper channel pixel water flow is moved toward the land surface instead.

To achieve this in a numerically balanced manner, the upstream channel pixel water flow calculated in the channel flow module is saved at each time step in an additional model variable. In the overland flow routing module, when the river head exceeds H_thres, the upstream channel pixel water flow is added to the surface water and removed from the channel inflow. In this numerical framework the channel inflow can be negative due to the removed water from the channel which has been brought back to the land. However, since this overbank flow process occurs only when the channel pixel is completely filled with water, that is with a water head of H_thres, this additional process intrinsically cannot generate negative streamflow in the channels. As shown in Figures 2B, 3 the tuning of H_thres results in a clear improvement of the simulated discharge results, which is detailed in the following sections.

A four-year period, beginning on 1^st January 2010, is chosen for the calibration of the WRF-Hydro discharge, using the setup of Section 2.2 enhanced with the overbank flow option of Section 2.3. This period is considered to be long enough to obtain robust calibration results with WRF-Hydro. No spin-up period is considered as no apparent spurious effects were detected at the beginning of the simulation, such as a spurious discharge peak that would be triggered by a too wet initial soil condition for example, as demonstrated by the daily timeseries shown in Figure 2. It is emphasized that the basin-averaged daily timeseries of precipitation simulated by WRF-Hydro is remarkably close to the IMERG observational dataset in Figure 2A, which justifies the suitability of this coupled setup for discharge calibration.

The simulated discharge is calibrated manually by tuning three sensitive parameters, namely the percolation parameter S, the overbank flow parameter H_thres, and the river roughness Manning coefficients. The reason behind this calibration strategy is that reducing H_thres smooths the discharge peaks but also removes too much water from the channels, which can be corrected by decreasing S to partially seal the soil column bottom and force more water to exfiltrate back to the surface. Reducing the Manning coefficients also tends to reduce water accumulation in the streams so that H_thres is less often reached, thus allowing to further modulate the overbank flow effect. The channel network displayed in Figure 1B has five stream orders in the Nzoia basin, which means that five Manning coefficients need to be calibrated in this case, as shown in Table 2.

For this calibration exercise the four-year WRF-Hydro simulation is run 8 × 6 = 48 times for all combinations among eight H_thres values (3, 4, 5, 6, 7, 8, 9, 10, in meters) and six S values (0.1, 0.05, 0.04, 0.03, 0.02, 0.01) with the default Manning coefficients Mdefault given in Table 2. Additionally, 8 × 2 × 4 = 64 runs are conducted for all combinations among the eight H_thres values, two S values (0.02, 0.01), and four sets of Manning coefficients M050, M045, M040, M035 given in Table 2. This makes a total of 112 calibration runs to identify a robust set of parameters' values for the Nzoia river discharge simulation with the WRF-Hydro setup of Sections 2.2 and 2.3.

The simulated discharge timeseries are assessed with three skill scores: the percentage bias, the Nash-Sutcliffe efficiency (NSE, Nash and Sutcliffe, 1970), and the Kling-Gupta efficiency (KGE, Gupta et al., 2009), as shown in Figure 3 and discussed in Section 3.1, with the best discharge result being displayed in Figure 2B. In this study the KGE formulation employs a bias scaling factor s_β set to 3 (see Eq. 11 of Gupta et al., 2009) to enhance the discrimination of biased results.

To conduct afforestation experiments for the Nzoia basin with WRF-Hydro, the calibrated WRF-Hydro from section 2.4 is used to generate a reference 20-year simulation from 1^st January 2001. The realism of the reference result is evaluated with the observational datasets presented in Section 2.1, as displayed in Figures 4–6 and discussed in Section 3.2. This evaluation is important to assess to which extent the model is suitable for a landcover change numerical experiment (de Noblet-Ducoudré et al., 2012; Sy et al., 2017).

The default landcover map considered in the reference simulation, which is displayed in Figure 7A, is modified in a synthetic manner according to three extreme afforestation scenarios. All grasslands and croplands areas in the Nzoia basin, which represents about 54% of the basin area, are replaced with a more natural vegetation category: savannas in scenario 1 (Figure 7B), woody savannas in scenario 2 (Figure 7C), and evergreen broadleaf forest in scenario 3 (Figure 7D). The 20-year simulation is reconducted for each of these three scenarios, and the differential results with respect to the reference are displayed in Figures 8–11, Table 3, and discussed in Section 3.3. The robustness of the areal climate change signal induced by the afforestation scenarios is assessed in a bootstrap approach by calculating for each model pixel the percentage of 17-year subset mean differences, among the 1,140 possible 17-year subset combinations out of the 20 simulated years, that have the same sign as the 20-year mean difference. The areal climate change signal is considered to be robust at the model pixels where this percentage exceeds 90%. The choice of 17 years for the subset size is motivated by the fact it is the highest subset size allowing to generate more than one thousand combinations.

The result of these extreme and purely synthetic afforestation experiments aims at giving an indication of the largest climatic impact that could be expected with an afforestation policy in the region of the Nzoia basin in tropical Africa, and to which extent this climatic impact depends on the vegetation cover type.

The results of the three skill scores in terms of percentage bias (Figure 3A), NSE (Figure 3B), and KGE (Figure 3C) from the experimental discharge simulations of 2010–2013 are summarized in Figure 3. Figure 3A confirms that reducing S (H_thres) increases (reduces) discharge. With the tested sets of Manning coefficients M050, M045, M040, and M035 (see Table 2), it also appears that discharge increases with decreasing river roughness, thus reducing the river-emptying effect of overbank flow. According to the NSE in Figure 3B, the best-simulated discharge performances are generally obtained for the configurations with a pronounced dry bias, such as for example with the four following parameters values sets (S = 0.01, H_thres = 4 m, M035), (S = 0.01, H_thres = 5 m, M040), (S = 0.01, H_thres = 5 m, M045), (S = 0.01, H_thres = 6 m, M050), with NSEs of 0.42, 0.42, 0.42, 0.39, respectively, and with biases of −24, −15, −20, −10%, respectively.

To distinguish the least biased best-performing configuration, a KGE formulation with an enhanced bias scaling factor is considered in Figure 3C. Accordingly, the best KGE is obtained for (H_thres = 6, S = 0.01, M045), with an NSE of 0.30 and a bias of −5%, as displayed in Figure 2B, so that (H_thres = 6, S = 0.01, M045) is the retained set of parameters values for the 20-year WRF-Hydro simulations presented in the following sections.

It is noted that an NSE of 0.30 remains relatively weak, but lies in the range of published discharge performances with WRF-Hydro for various river basins in the world (e.g., Senatore et al., 2015; Arnault et al., 2016, 2018; Kerandi et al., 2017; Rummler et al., 2019; Zhang et al., 2019; Camera et al., 2020; Fersch et al., 2020; Li et al., 2020; Sofokleous et al., 2023). As highlighted by Senatore et al. (2015), simulated discharge performance with WRF-Hydro is intrinsically limited by the quality of simulated precipitation, which is a drawback of the atmospheric-hydrological coupled modeling approach. Still, the WRF-Hydro results over the calibration period are realistic (Figure 2), which gives us confidence that the climate response to landcover change will be plausible.

Climatological results for the period 2001–2020 are evaluated with maps of annual precipitation and mean temperature in Figure 4, with maps of mean energy fluxes in Figure 5, and with mean monthly timeseries for the Nzoia basin in Figure 6.

The annual precipitation pattern from WRF-Hydro in Figure 4C resembles that from IMERG in Figure 4A, with a southwest-northeast gradient and enhanced amounts toward Lake Victoria. However, as can be seen in the bias map of Figure 4E, WRF-Hydro is drier in the lowlands and much wetter in the highlands. The fact that IMERG does not display clearly enhanced precipitation in mountainous areas is a well-known limitation of satellite-derived products (Pradhan et al., 2022). Furthermore, the 1,000 mm of annual precipitation in the parts of the Nzoia basin close to Lake Victoria described by Githui et al. (2010) cannot be seen in IMERG, but do appear in WRF-Hydro. Still, according to the timeseries in Figure 6A, IMERG captures the three precipitation peaks, exactly in April-May, in August, and in October as mentioned by Githui et al. (2010), whereas WRF-Hydro only exhibits two precipitation peaks in May and September with an overall precipitation bias of +13%. This precipitation overestimation is also reflected in the mean monthly discharge timeseries of Figure 6B, with a discharge bias of +16%, despite an underestimation of low flow between January and March. These results indicate that the simulated climate in the Nzoia basin with the calibrated WRF-Hydro setup presented in Section 2 is too wet, even though this positive bias was not so pronounced during the period 2010–2013 used for discharge calibration (Figure 2A). As highlighted by the multi-physics ensembles of Otieno et al. (2018, 2020), the WRF precipitation bias in the region northeast of Lake Victoria largely depends on the choice of cumulus convection scheme, which suggests that future research efforts should focus on further adapting the atmospheric modeling setup of Section 2 to improve WRF-Hydro precipitation results.

The mean temperature patterns from WRF-Hydro in Figure 4D and from CRU in Figure 4B are remarkably similar, however, there is a cold bias over much of the region (Figure 4F). This cold bias is even more pronounced in the mean monthly timeseries (Figure 6C), between 2 and 3 K, which is potentially due to unresolved high-elevation cold temperatures in the CRU data. The large amounts of simulated precipitation in September–October may also have enhanced the cold bias at that time.

The mean values and patterns of net radiation, sensible, and latent heat flux from WRF-Hydro in Figures 5D–F and from FLUXCOM in Figures 5A–C are relatively similar, especially in the region of the Nzoia basin. This is confirmed by the timeseries in Figure 6D which reveals only small discrepancies. The overestimation of latent heat flux from October to April could be related to the large precipitation overestimation in September–October and too wet soils in the subsequent months. The underestimation of latent heat flux from June to August could be related to the underestimation of net radiation during that period, and potentially also inaccuracy in the vegetation parameters. Nevertheless, these results indicate that the WRF-Hydro setup captures the surface energy partitioning, and especially the evapotranspiration flux, realistically.

The climatological differences induced by the extreme and synthetic afforestation scenarios presented in Section 2.5 and Figure 7 are visualized as differential maps of annual precipitation, mean temperature, mean energy fluxes, mean evaporation fluxes, and mean runoff fluxes in Figures 8–10, and as differential mean monthly timeseries of the terrestrial water budget terms for the Nzoia basin in Figure 11. The surface runoff flux shown in Figures 10, 11 is calculated as the sum of overland flow and channel inflow and is added to the WRF-Hydro outputs following the procedure detailed in Arnault et al. (2019).

Remarkably, all savannas, woody savannas and evergreen broadleaf forest scenarios result in a robust precipitation increase in the Nzoia basin and near surrounding areas toward the South, West and East, as displayed in Figures 8A, C, E. There is also a robust but small precipitation decrease over the highest mountains of the Nzoia basin and also outside of the basin toward the Southwest and Southeast, which indicates that the tested afforestation scenarios in the Nzoia basin increase precipitation locally at the expense of other regions, to some extent.

Concerning temperature, only the evergreen broadleaf forest scenario in Figure 8F robustly induces a cooling on the annual scale in the Nzoia basin, up to 0.5 K, whereas the savannas and woody savannas scenarios in Figures 8B, D rather robustly induce a very slight warming in the Nzoia basin. The reason for this model behavior can be partially elucidated with the maps of mean energy flux differences in Figure 9.

Remarkably, savannas, woody savannas, and evergreen broadleaf forest scenarios result in a robust slight increase in net radiation flux in Figures 9A, D, G, a robust large increase in sensible heat flux in Figures 9B, E, H, and a robust large decrease in latent heat flux in Figures 9C, F, I in the Nzoia basin. The slight increase in net radiation flux, which is related to the cloud cover change associated with the precipitation increase, the surface albedo change associated with the landcover change and a skin temperature change, is much smaller in amplitude compared to the changes in sensible and latent heat fluxes and is not discussed further.

The large decrease in latent heat flux, especially for the savannas and woody savannas scenarios in Figures 9C, F, is surprising, but consistent with the large sensible heat increase in Figures 9B, E, and the slight temperature increase in Figures 8B, D. For the case of the evergreen broadleaf forest scenario, the latent heat-decrease in Figure 9H and sensible heat increase in Figure 9I are much less pronounced. Temperature still decreases in this case (Figure 8F), which is related to the fact that the energy fluxes in Noah-MP are calculated above the canopy, whereas the temperature includes a below-canopy cooling effect that is much more efficient for the evergreen broadleaf forest type (Noah-MP, Niu et al., 2011). The in-situ measurements of Ceperley et al. (2017) in the West African Sahelian region confirm that savannas can increase the sensible heat flux in comparison to agricultural land, although Ceperley et al. (2017) attributed this heating enhancement to a larger number of exposed rocks rather than a diminution of latent heat as in our modeled case.

Among our afforestation modeling experiments the sensible heat flux increase is largest in the case of the woody savannas scenario in Figure 9E, and is associated with the largest precipitation increase in Figure 8C. This is an indication that simulated precipitation increase in these experiments is the result of enhanced atmospheric heating and enhanced convective instability induced by a reduction of the latent heat, which corresponds to the negative feedback process described for example by Taylor et al. (2012).

The latent heat flux decrease triggered by the afforestation scenarios is related to a robust large decrease in soil evaporation flux in the Nzoia basin (see Figures 10A, E, I), that is only partially balanced by a robust smaller increase in plant evapotranspiration flux (see Figures 10B, F, J). Indeed, the calculation of the soil evaporation in Noah-MP also considers the soil evaporation which happens below the canopy. Thus, according to Noah-MP the savannas and woody savannas vegetation types are more efficient in reducing soil evaporation loss compared to grasslands and croplands. This soil water loss-reducing effect is almost completely canceled in the case of the evergreen broadleaf forest scenario, where the plant evapotranspiration flux in Figure 10J is very much enhanced as expected. It is worth noting that a comprehensive panel of in-situ observations would be needed to verify these detailed model behaviors (Ceperley et al., 2017; Bliefernicht et al., 2018).

The overall robust increases in surface runoff flux (Figures 10C, G, K), and underground runoff flux (Figures 10D, H, L) in the Nzoia basin mainly arise from the precipitation increase triggered by the negative feedback process mentioned above. To a lesser extent the evapotranspiration flux decrease plays a role in this runoff increase, following the above-discussed water balance-based mechanism proposed for example by Achugbu et al. (2022) and more specifically for the Nzoia basin by Githui et al. (2010). The mean-monthly timeseries of Figure 11, Table 3 further illustrate the runoff increase response to a dual precipitation increase—evapotranspiration decrease, a process much exacerbated in the case of the woody savannas' scenario in Figure 11B and much attenuated in the case of the evergreen broadleaf forest scenario in Figure 11C. In the case of the woody savannas' scenario, water resources would be much increased, by about 13% for precipitation, 8% for surface runoff and 12% for underground runoff or groundwater recharge, as summarized in Table 3. The relatively narrow error margins displayed in Figure 11 and indicated in Table 3 demonstrate the significance of the simulated climate change signals presented in this article.

It is recognized that the tested afforestation scenarios are idealized and difficult to apply in the real world, but they demonstrate the potential of landcover change to modulate climate change. More generally, these results highlight the importance of the coupled modeling approach for a comprehensive assessment of landcover change impacts on climate (Sy and Quesada, 2020).