Regional water cycle sensitivity to afforestation: synthetic numerical experiments for tropical Africa

allows the assessment of potential regional climate impacts of three afforestation strategies. In all three cases, the afforestation-induced decrease in soil evaporation is larger than the afforestation-induced increase in plant transpiration, thus increasing sensible heat ﬂux and triggering a localized negative feedback process leading to more precipitation and more runoff. This effect is more pronounced with the woody savannas experiment, with 7% less evapotranspiration, but 13% more precipitation, 8% more surface runoff, and 12% more underground runoff predicted in the Nzoia basin. This study demonstrates a potentially large impact of afforestation on regional water resources, which should be investigated in more detail for policy making purposes.

A orestation as a climate change mitigation option has been the subject of intense debate and study over the last few decades, particularly in the tropics where agricultural activity is expanding.However, the impact of such landcover changes on the surface energy budget, temperature, and precipitation remains unclear as feedbacks between various components are di cult to resolve and interpret.Contributing to this scientific debate, regional climate models of varying complexity can be used to test how regional climate reacts to a orestation.In this study, the focus is on the gauged Nzoia basin ( , km ) located in a heavily farmed region of tropical Africa.A reanalysis product is dynamically downscaled with a coupled atmospheric-hydrological model (WRF-Hydro) to finely resolve the land-atmosphere system in the Nzoia region.To overcome the problem of Nzoia river flooding over its banks we enhance WRF-Hydro with an overbank flow routing option, which improves the representation of daily discharge based on the Nash-Sutcli e e ciency and Kling-Gupta e ciency (from − .to ., and − .to ., respectively).Changing grassland and cropland areas to savannas, woody savannas, and evergreen broadleaf forest in three synthetic numerical experiments allows the assessment of potential regional climate impacts of three a orestation strategies.In all three cases, the a orestation-induced decrease in soil evaporation is larger than the a orestation-induced increase in plant transpiration, thus increasing sensible heat flux and triggering a localized negative feedback process leading to more precipitation and more runo .This e ect is more pronounced with the woody savannas experiment, with % less evapotranspiration, but % more precipitation, % more surface runo , and % more underground runo predicted in the Nzoia basin.This study demonstrates a potentially large impact of a orestation on regional water resources, which should be investigated in more detail for policy making purposes.
A better understanding of how changes to landcover modify climate at the regional scale is crucial for developing and assessing climate change mitigation strategies (e.g., Wulfmeyer et al., 2014;Harper et al., 2018;Roe et al., 2019;Duveiller et al., 2020).For afforestation, besides a greater carbon sequestration potential, local temperature changes would be dominated by two competing processes that are enhanced evaporative cooling and enhanced radiative heating due to a lower albedo (Bonan, 2008).As pointed out by Breil et al. (2021), the evaporative cooling effect triggered by afforestation is also modulated by the saturation deficit of the near-surface atmosphere, which adds additional complexity to afforestation effects on temperature.Wulfmeyer et al. (2014), Liu et al. (2023), and Wang et al. (2023) brought evidence that afforestation in an arid region can trigger more precipitation, although enhanced precipitation can be associated with deforestation as well (Achugbu et al., 2023).
Land-atmosphere processes are represented in a simplified manner in regional climate models, which allow investigating the climatic response to afforestation at a regional scale (Davin et al., 2020).However, this modeling approach is hindered by large uncertainties in the numerical results that likely arise from the simplifying assumptions that are among its assets (e.g., de Noblet-Ducoudré et al., 2012;Lejeune et al., 2017;Davin et al., 2020).In the case of a model ensemble experiment for Europe, Davin et al. (2020) could not identify a robust temperature change related to afforestation.The effect of landcover change on precipitation in such numerical experiments is even more uncertain with low signal-to-noise ratios (e.g., Laux et al., 2017).These model uncertainties could be partly alleviated by improving the realism of the physical processes included in the model itself (de Noblet-Ducoudré et al., 2012).
Regional climate models' representation of land-atmosphere interactions is usually limited to a vertical transfer of energy and water at the land-atmosphere interface.This is a crude approximation of terrestrial hydrology at the river basin scale.Studies with coupled atmospheric-hydrological models, such as with the hydrologically enhanced Weather Research and Forecasting model WRF-Hydro (Gochis et al., 2020), have shown that the consideration of surface and subsurface lateral flow improves the simulated runoff and modifies the atmospheric branch of the simulated water cycle, with generally more evapotranspiration and more precipitation (e.g., Rummler et al., 2019;Zhang et al., 2019;Fersch et al., 2020;Arnault et al., 2021;Furnari et al., 2022).
Applying the uncoupled version of the WRF-Hydro hydrological model, Achugbu et al. (2022) investigated the impact of landcover change scenarios on evapotranspiration and runoff in a Sahelian West African river basin.Achugbu et al. (2022) confirmed that the increase of evapotranspiration associated with afforestation leads to a decrease in river discharge.
In another hydrological model experiment on landcover change for the Nzoia basin in tropical Africa, Githui et al. (2010) attributed the runoff increase between 1973 and 2001 to an increase of agricultural area and a reduction of forested area.However, the uncoupled hydrological modeling approach does not account for potential feedbacks between the terrestrial hydrological system and the overlying atmosphere, especially for precipitation.As such, the landcover change effect on climate is only partially addressed in such uncoupled model studies (Githui et al., 2010).
A more comprehensive effect of landcover change on both the terrestrial and atmospheric branches of the water cycle can in principle be evaluated with the atmospheric-hydrological coupled version of WRF-Hydro, although no such study is known to us so far.Our study therefore aims at presenting the first landcover change experiment with the coupled WRF-Hydro, exemplarily investigated for the hypothetical case of afforestation in the Nzoia basin.The choice of the study region is motivated by (1) its location in the tropics, (2) an expanding agricultural area which raises the question of the potential benefits of an afforestation policy, and (3) the availability of river discharge data to assess the ability of the model to represent water resources in this region.The objective of the study is to evaluate, from a coupled atmospheric-hydrological modeling point of view, whether a climate change mitigation policy such as afforestation in a tropical region would be hydrologically significant, particularly in terms of precipitation and river discharge.
In the following, Section 2 details the study region and available datasets, the WRF-Hydro model, a model extension to account for overbank flow, the calibration and evaluation strategy and the landcover change experiments.The results are discussed in Section 3, and concluding remarks are finally given in Section 4.
. Materials and method

. . Study region and datasets
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 2 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. . .Coupled WRF-Hydro modeling setup 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 sixhourly 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 FIGURE (A) Daily timeseries of Nzoia basin-averaged precipitation P (in mm/day) for a four-year period spanning from st January to st December and filtered with a -day Gaussian filter, as derived from the observational product IMERG and from the inner domain of the default and calibrated WRF-Hydro simulations.(B) Daily timeseries of discharge Q (in m /s) at the Nzoia basin outlet for a four-year period spanning from st January to st December , as derived from the gauge measurement and from the default and calibrated WRF-Hydro simulations.The calibrated simulation employs the modified WRF-Hydro source code presented in Section .with the best set of parameters values according to the calibration exercise illustrated in Figure , that is (H thres = m, S = ., M ).
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 The first column gives a short description of the experiments conducted.The next six columns give the values of the percolation parameter (SLOPE), the runoff-infiltration parameter (REFKDT), the lateral hydraulic conductivity factor (LKSATFAC), the overland flow roughness (OVROUGHRTFAC), the retention depth factor (RETDEPRTFAC), and the river roughness Manning coefficients (MannN).The two last columns give the Nash-Sutcliffe Efficiency (NSE) and percentage bias of simulated discharge at the Nzoia outlet gauge with respect to observation.
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.

. . Two-way land-river water flow model extension
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 ).
Frontiers in Climate frontiersin.org ./fclim. .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.

. . River discharge calibration strategy
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.

. . Landcover change numerical experiment strategy
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. .Results and discussion
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., 2016Arnault et al., , 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.

. . Climatology evaluation
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 to st December , (B) a nine-year period spanning from st January to st December , (C) a -year period spanning from st January to st December , and (D) for a thirteen-year period spanning from st January to st December , when the corresponding observational data is available.
so pronounced during the period 2010-2013 used for discharge calibration (Figure 2A).As highlighted by the multi-physics ensembles of Otieno et al. (2018Otieno et al. ( , 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.

. . A orestation impacts
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 FIGURE Maps of (A, C, E) annual precipitation di erence P displayed as a percentage of the annual IMERG precipitation and (B, D, F) mean temperature di erence T (in K) between the reference WRF-Hydro simulation and the (A, B) savannas, (C, D) woody savannas, (E, F) evergreen broadleaf forest WRF-Hydro simulations.Mean di erential values are calculated as the mean value from the a orestation experiments minus the mean value from the reference experiment, for a -year period spanning from st January to st December .The annual di erential precipitation maps (A, C, E) have the same color scale indicated by the color bar on the right side of (E), and the mean di erential temperature maps (B, D, F) also have the same color scale indicated by the color bar on the right side of (F).In all maps, the crossed areas correspond to the areas where the climate change signal is not robust, that is where < % of the , subsets of -year mean di erences have the same sign as the -year mean di erence.
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  (G-I) evergreen broadleaf forest WRF-Hydro simulations.Mean di erential values are calculated as the mean value from the a orestation experiments minus the mean value from the reference experiment, for a -year period spanning from st January to st December .All di erential energy flux maps have the same color scale indicated by the color bar on the right side of (G).In all maps, the crossed areas correspond to the areas where the climate change signal is not robust, that is where < % of the , subsets of -year mean di erences have the same sign as the -year mean di erence.
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 increaseevapotranspiration decrease, a process much exacerbated in the FIGURE Maps of (A, E, I) annual soil evaporation di erence E soil (in % of IMERG annual precipitation), (B, F, J) annual plant evapotranspiration di erence E plant (in % of IMERG annual precipitation), (C, G, K) annual surface runo di erence R surface (in % of IMERG annual precipitation), and (D, H, L) annual underground runo di erence R ground (in % of IMERG annual precipitation), between the reference WRF-Hydro simulation and the (A-D) savannas, (E-H) woody savannas, (I-L) evergreen broadleaf forest WRF-Hydro simulations.In (A, E, I) the soil evaporation is calculated as the sum of the soil evaporations outside and below the canopy.In (B, F, J) the plant evapotranspiration is calculated as the sum of the plant transpiration and intercepted canopy water evaporation.Mean di erential values are calculated as the mean value from the a orestation experiments minus the mean value from the reference experiment, for a -year period spanning from st January to st December .All di erential water flux maps have the same color scale indicated by the color bar on the right side of (I).In all maps, the crossed areas correspond to the areas where the climate change signal is not robust, that is where < % of the , subsets of -year mean di erences have the same sign as the -year mean di erence.

FIGURE
Climatological monthly timeseries of di erential water flow terms of the terrestrial water budget spatially integrated in the area of the Nzoia basin, namely the di erential precipitation P (in m /s), the di erential total evapotranspiration E (in m /s), the di erential surface runo R surface (in m /s), the di erential underground runo R ground (in m /s), and the di erential soil water storage change StorageChange (in m /s), between the reference WRF-Hydro simulation and the (A) savannas, (B) woody savannas, (C) evergreen broadleaf forest WRF-Hydro simulations.Climatological monthly di erential values are calculated as the climatological monthly value from the a orestation experiments minus the climatological monthly value from the reference experiment.In each panel, the lines indicate the -year average monthly di erences for the period spanning from st January to st December .The shaded areas represent an error margin calculated as the spread obtained from the , subsets of -year average monthly di erences among the simulated years.The legend is provided in the upper-right side of (A).
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).

. Summary and concluding remarks
The first application of the coupled atmospheric-hydrological model WRF-Hydro to synthetic landcover change numerical experiments was presented, exemplarily realized for the case of the Nzoia basin in Kenya.As part of the calibration, WRF-Hydro was enhanced with an overbank flow option.With this additional realism, the calibrated discharge performance was greatly improved TABLE Simulated annual water flow changes in km /year, and relative annual water flow changes as a percentage of the reference annual precipitation, induced by a orestation type for precipitation P, total evapotranspiration E, surface runo R surface , underground runo R ground , as derived from the di erences between the reference WRF-Hydro simulation and the savannas, woody savannas, evergreen broadleaf forest WRF-Hydro simulations during the period spanning from st January to st December .1.9 (±0.2) 7.9 (±1.0) 3.In this table a differential value is defined as the value from an afforestation experiment minus the value from the reference experiment.The numbers in parenthesis give the error margin estimated as the range of the results from the 1,140 combinations of 17-year subsets among the 20 simulated years.

Savannas
in the study basin.Then, the calibrated model was run for a 20year historical period using the default landcover map, and three modified landcover maps according to the following afforestation scenarios: grasslands and croplands of the Nzoia basin replaced by (1) savannas, (2) woody savannas, (3) evergreen broadleaf forest.These synthetic afforestation scenarios led to robust results in the area of the Nzoia basin.A large reduction of the soil evaporation was only partially balanced by a smaller increase in plant evapotranspiration, which reduced the latent heat flux, increased the sensible heat flux, and finally triggered more precipitation, surface runoff, and groundwater recharge.This mechanism was particularly pronounced with the woody savannas scenario, but much attenuated with the evergreen broadleaf forest scenario, which suggests that an afforestation policy can have a large impact on the regional climate and water resources availability, depending on the planting strategy.For decision-makers, the woody savannas scenario would imply to build more infrastructures to store water and bring it to less humid parts of the country, protect the population from increased flood risks, and relocate the farmers having lost their land for cropping.
To increase the value of such idealized synthetic numerical results for decision-makers, observation and model improvement efforts should be jointly conducted.The establishment of insitu measurement stations, as realized for example in West Africa (Bliefernicht et al., 2018), is crucial for validating the representation of physical processes in the model for a specific region.A parameter sensitivity study would be beneficial to better understand the physical mechanisms associated with the modeled impacts of afforestation on climate.Enhancing the horizontal resolution of the model grid would enable the implementation of more realistic landcover change scenarios, although this would require several additional tests to identify suitable sets of physics schemes for the atmospheric modeling, and to recalibrate the hydrological parameters.Additional model developments, such as soil parameterization (Zhang et al., 2023), groundwater coupling (Rummler et al., 2022), vegetation dynamics (Glotfelty et al., 2021;Warrach-Sagi et al., 2022), carbon-nitrogen biogeochemistry (Chang et al., 2018), should also be considered to further enhance the realism of simulated land-atmosphere interactions with coupled WRF-Hydro.

FIGURE
FIGURETopography (in m above sea level) of (A) the outer and inner domains D and D at and km resolution, respectively, and (B) the sub-domain D sub at km resolution coupled with D for water routing computations.The black rectangle in (A) gives the location of D .The red contour line, blue curved lines, and black cross in (B) give the locations of the Nzoia basin, main rivers, and outlet gauge station EEF , respectively.The topography color scale provided by the color bar on the left side of (A) is the same for both panels.

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 biosphereatmosphere 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.

FIGURE
FIGUREPerformance of the simulated discharge at the Nzoia river gauge location with the WRF-Hydro setup of Sections .and ., in terms of (A) percentage bias, (B) Nash-Sutcli e e ciency NSE and (C) Kling-Gupta e ciency KGE, as a function of the overbank flow parameter H thres (in m) on the x-axis, percolation parameter S (-) with the colored lines, and river roughness Manning coe cients M (-) of Table with the plain lines and symbolled dashed lines.The three configurations with a bias lower than +/-% and a positive NSC, namely the configurations using (H thres = m, S = ., M), (H thres = m, S = ., M ), and (H thres = m, S = ., M ), are highlighted with bolded symbols.KGE in (C) is calculated with a bias scaling factor set to (see Section .) to favor the less biased good result, in this case with (H thres = m, S = ., M).

FIGURE
FIGUREMaps of (A) annual precipitation P (in mm/year) from the observational product IMERG, (B) mean temperature T (in K) from the observational product CRU, (C, D) annual precipitation and mean temperature from the inner domain of the WRF-Hydro reference simulation, (E) precipitation bias (in %) between IMERG and WRF-Hydro reference, and (F) temperature bias (in K) between CRU and WRF-Hydro reference.Mean values are calculated for a -year period spanning from st January to st December .The annual precipitation maps (A, C) have the same color scale indicated by the color bar on the right side of (C), and the mean temperature maps (B, D) also have the same color scale indicated by the color bar on the right side of (D).

FIGURE
FIGURE Maps of (A) mean net radiation flux R net (in W/m ), (B) mean sensible heat flux H sensible (in W/m ) and (C) mean latent heat flux H latent (in W/m ) from the observational product FLUXCOM, (D-F) mean net radiation, sensible and latent heat fluxes from the inner domain of the WRF-Hydro reference simulation, (G-I) net radiation, sensible and latent heat fluxes biases (in %) between FLUXCOM and WRF-Hydro reference.Mean values are calculated for a thirteen-year period spanning from st January to st December , when the FLUXCOM data is available.The energy flux maps (A-F) have the same color scale indicated by the color bar on the right side of (D), and the energy flux bias maps (G-I) also have the same color scale indicated by the color bar on the right side of (G).

FIGURE
FIGUREClimatological monthly timeseries of (A) Nzoia basin-averaged precipitation P (in mm/month), (B) Nzoia basin outlet discharge Q (in m /s), (C) Nzoia basin-averaged temperature T (in K), (D) Nzoia basin-averaged net radiation, sensible and latent heat fluxes R net , H sensible , H latent (in W/m ), as derived from (A) the observational product IMERG, (B) the gauge measurement, (C) the observational product CRU, (D) the observational product FLUXCOM, and from the inner domain of the WRF-Hydro reference simulation.Climatological monthly values of (A) precipitation, (B) discharge, (C) temperature, (D) energy flux values are calculated from (A) a -year period spanning from st January to st December , (B) a nine-year period spanning from st January to st December , (C) a -year period spanning from st January to st December , and (D) for a thirteen-year period spanning from st January to st December , when the corresponding observational data is available.

FIGURE
FIGURE Maps of landcover distribution in the inner domain following the MODIS classification, as used in (A) the reference WRF-Hydro simulation, (B) the savannas WRF-Hydro simulation, (C) the woody savannas WRF-Hydro simulation, and (D) the evergreen broadleaf forest WRF-Hydro simulation.The black contour line in all panels indicates the location of the Nzoia basin where the landcover change experiment is conducted.Landcover outside the Nzoia basin is maintained as in the reference.The color classification of MODIS landcover classes is provided in the upper right side of the figure.

FIGURE
FIGUREMaps of (A, D, G) mean net radiation flux di erence R net (in W/m ), (B, E, H) mean sensible heat flux di erence H sensible (in W/m ), and (C, F, I) mean latent heat flux di erence H latent (in W/m ) between the reference WRF-Hydro simulation and the (A-C) savannas, (D-F) woody savannas, (G-I) evergreen broadleaf forest WRF-Hydro simulations.Mean di erential values are calculated as the mean value from the a orestation experiments minus the mean value from the reference experiment, for a -year period spanning from st January to st December .All di erential energy flux maps have the same color scale indicated by the color bar on the right side of (G).In all maps, the crossed areas correspond to the areas where the climate change signal is not robust, that is where < % of the , subsets of -year mean di erences have the same sign as the -year mean di erence.

P
TABLE List of model parameter sets employed in a series of preliminary experiments with the WRF-Hydro setup of section ., and their corresponding Nzoia discharge results for a -year simulation period beginning on st January .
TABLE List of river roughness Manning coe cients as a function of stream order tested in the WRF-Hydro calibration experiment of Section . .