Leveraging soil diversity to mitigate hydrological extremes with nature-based solutions in productive catchments: an application and insights into the way forward

Nature-based solutions (NbSs) are increasingly explored to mitigate floods and agricultural droughts in productive catchments. However, the influence of local factors, such as soil characteristics, on the effectiveness of these interventions is often overlooked, as commonly used modelling approaches oversimplify interactions across soil–water processes. This study presents an innovative approach to model NbS scenarios, revealing the influence of the spatial variability of soil properties on their performance. A fully distributed, physically based hydrological model was used to represent NbSs at the catchment scale, explicitly simulating soil–water fluxes (e.g. infiltration, evapotranspiration, runoff). This model uses measurable local parameters (e.g. topography, soil properties, vegetation) to capture small-scale hydrological processes and represents NbS scenarios via their adjustments. Simulations were conducted for two catchments: one agricultural and one forested. In the agricultural catchment, measures included hedgerows, reduced tillage and soil pitting in maize crops. In the forested catchment, NbSs focused on forest diversification, practices aimed at limiting soil compaction and the restoration of peatlands. NbS performance in mitigating floods and agricultural droughts was assessed using spatial indicators. The models accurately reproduced discharge and saturated zone dynamics, capturing natural soil drainage characteristics and their interaction with NbS effectiveness. Results highlight that NbS performance strongly depends on natural soil drainage characteristics. In well-drained soils, hedgerows enhanced infiltration by improving soil hydraulic properties and boosting evapotranspiration. In waterlogged soils, more limited gains were observed. Well-drained soils offer co-benefits to agricultural drought mitigation, as they are prone to water deficits that NbSs can mitigate. Future evaluation of NbS effectiveness should recognize the spatial variability in performance to better inform their selection and placement. This study may serve as a basis for discussion and action, supporting decision-makers in implementing these measures in a coherent territory-based approach. It also reveals current knowledge gaps and identifies avenues for future research to refine NbS effectiveness assessments, such as strengthening the availability of spatially distributed data and advancing uncertainty analysis.

1 Introduction

Central Europe faced consecutive abnormally dry and hot summers in 2018, 2019 and 2020 (van der Wiel et al., 2023). The drought event of 2018–2019 was considered unprecedented in the last 250 years (Hari et al., 2020), adversely impacting agricultural productivity (Toreti et al., 2019). Droughts are classified as meteorological, agricultural, hydrological or socioeconomic. Meteorological drought refers to a deficiency in rainfall, while agricultural drought, driven by soil moisture deficits, occurs when soil water availability does not meet plant requirements. It is therefore influenced by vegetation characteristics, as species differ in their sensitivity to the same soil moisture deficit. Extended rainfall shortages may result in hydrological drought, affecting streamflow and aquifers, eventually causing socioeconomic drought when water supply fails to meet demand (Mishra and Singh, 2010).

Conversely, the same region, comprising western Germany, eastern Belgium, Luxembourg, and the Netherlands, experienced a very wet 2021 summer and an extreme precipitation event between 13 and 16 July, with unprecedented accumulations. Combined with already nearly saturated soil conditions, this extreme rainfall event led to major floods, making it one of the most severe natural catastrophes in Europe of the past half-century (Journée et al., 2023).

Due to the anthropogenic global warming and the demographic expansion, the likelihood of occurrence and the impact of droughts and floods are expected to be exacerbated (Aalbers et al., 2023; Dottori et al., 2023; Hari et al., 2020), consequently increasing the associated damage costs (Naumann et al., 2021), currently estimated at EUR 9 billion for droughts (Cammalleri et al., 2020) and EUR 7.6 billion for floods (Dottori et al., 2023) annually in the European Union.

In order to mitigate and potentially reverse these adverse effects, investments in a transformative adaptation of human systems are required (Pörtner et al., 2022). Given the uncertainty of future climate conditions, these investments could prioritize soft strategies, “no regret” approaches that yield benefits irrespective of climate change, as well as reversible and flexible options, providing significant social, economic and environmental benefits (Hallegatte, 2009). In line with this philosophy, the recently introduced umbrella concept of nature-based solutions (NbSs) is becoming increasingly popular among funders, researchers, policy makers and practitioners (Nesshöver et al., 2017). NbSs could be referred to as solutions rooted in natural processes and ecosystems, designed to address a spectrum of societal and environmental challenges, including hydrometeorological risk reduction (Ruangpan et al., 2020). NbS is a relatively broad concept that can be further specified through more targeted concepts, depending on the goals addressed, the type of ecosystem considered and the scale of intervention. For instance, natural water retention measures (NWRMs) and natural/small water retention measures (NSWRMs) refer to NbSs specifically aimed at water management. Although all these concepts overlap considerably and their diverse terminology can be confusing, Magnier et al. (2024) offer a structured, targeted analysis that clarifies distinctions between them.

While NbS literature is abundant on runoff and flood risk reduction in urban areas (Ruangpan et al., 2020), there has been limited exploration of NbSs as drought mitigation strategies (Yimer et al., 2024). Furthermore, potential combined effects, such as synergies and trade-offs of NbSs on floods and droughts, remain largely unexplored (Fennell et al., 2023b; Penning et al., 2023). NbSs are frequently considered effective in addressing hydrological droughts (Fennell et al., 2023a) and less so for agricultural droughts. However, NbSs are increasingly recognized as a useful concept extending beyond the urban/riverine context, such as in agricultural or forest settings (Hanson et al., 2020) to improve the drought resilience of these ecosystems while helping to reduce flood risks. In these contexts, NbSs may include soil–vegetation solutions aimed at enhancing soil health, functions and ecosystem services with agronomic and/or forest management practices, and/or land restoration (Keesstra et al., 2018). However, in these alternative land use contexts, such as agriculture or forestry, the effectiveness of NbSs might highly depend on local characteristics such as the infiltration capacity of the soil (Fennell et al., 2023a; Penning et al., 2023). Indeed, soil hydraulic properties and their spatial distribution play a crucial role in controlling small-scale (soil-profile-scale) hydrological processes within a watershed (Vereecken et al., 2022) and, consequently, the small-scale effectiveness of NbSs influencing hydrological processes within a watershed.

This fine-scale spatial variability raises important questions about the appropriate methods and scales to assess and monitor the effectiveness of these interventions before and after implementation. One popular approach for assessing the potential effectiveness of NbSs before implementation for floods and droughts involves modelling hydrological, hydraulic and water balance processes. This approach often involves comparing the modelled outcomes of a baseline scenario of the current conditions to a target scenario (Possantti and Marques, 2022). Several hydrological models have been used to simulate the effectiveness of NbSs, including, among others, SWAT+ (European Commission, 2020), LISFLOOD (Institute for Environment and Sustainability: Joint Research Centre et al., 2012), HEC-HMS (Guido et al., 2023) and MIKE SHE (Fennell et al., 2023b). These models generally integrate different modules to calculate different water fluxes (infiltration, runoff, evapotranspiration) and can be based on empirical, conceptual or physically based approaches. For example, SWAT and HEC-HMS use the well-known Soil Conservation Service (SCS) method to calculate water infiltration into the soil based on an empirical formula and the soil curve number (CN). Other models, such as MIKE SHE, are process-based and use physical equations, such as the Richards equation, to model infiltration. Some models allow for fully distributed simulations, while others are lumped models or adopt a non-distributed approach. Each model has its own advantages and limitations, and the choice depends on the specific needs of the study in terms of modelling capabilities, as well as the required spatial and temporal scale and resolution, and the available data (Kumar et al., 2021b). In the urban context, the hydrological effectiveness of NbSs is often assessed downstream through simulated discharge series using empirical or conceptual (sometimes lumped) rainfall runoff models (Possantti and Marques, 2022; Kumar et al., 2021b). However, this type of evaluation does not differentiate the effectiveness of individual upstream interventions, ignoring the significance of the spatial arrangement of NbSs within the watershed. This oversight presupposes that it has a minimal impact on their effectiveness downstream, as identified in the case of green roofs (Qiu et al., 2021).

In other land uses such as agricultural, this consideration might not be true and might require spatially distributed, physically based hydrological modelling approaches, able to rank the most effective NbS–location combinations in a given context (Brauman et al., 2022). Indeed, these types of models enable the virtual implementation of a set of NbSs at the catchment scale with a specific spatial arrangement. Water fluxes, calculated on a variable-resolution grid, are derived from physically meaningful equations that rely on physically meaningful local parameters that can be measured (topography, geology, soil hydraulic properties, surface roughness and vegetation). These factors influence small-scale (grid-cell-scale) hydrological processes, such as infiltration. Through the spatial and temporal adjustment of these parameters, scenarios of NbSs can be represented. Thus, the physical response (and effectiveness) of an intervention implemented at a specific location that might affect the global catchment's behaviour can be modelled.

This article introduces a reproducible methodology for developing a physically based surface hydrological model designed to evaluate the impact of NbS scenarios at the catchment scale. Serving as a foundational step, this modelling framework facilitates a series of innovative post-processing analyses to assess NbS hydrological effectiveness. The paper addresses four key objectives: (i) building a modelling approach that represents both the spatial variation of soil characteristics (focusing on soil natural drainage) and spatialized NbS scenarios, as well as their interactions; (ii) moving beyond a simple analysis of outflow discharges by emphasizing the influence of spatial variability in soil characteristics on the effectiveness of hydrological NbS measures inside the watershed; (iii) evaluating the impact of NbSs on different hydrological extremes, specifically floods and agricultural droughts; and (iv) investigating synergies and trade-offs between floods and agricultural drought mitigation through NbS implementation.

2 Materials and methods

2.1 Study area

The Vesdre catchment was selected as the study area because it was deeply impacted during the 2021 mid-July flood event in Belgium. More than 20 people lost their lives in this valley (Dewals et al., 2021), and many residents lost some or all of their belongings. This tragic event remains deeply imprinted in the memory of the region. Unfortunately, climate projections indicate that such extreme events could occur two to three times by 2050 in this area. At the same time, these projections suggest that summers will become increasingly hotter and drier on average (Fettweis, 2024). This situation is forcing local authorities to rethink the valley's future to enhance its resilience to the growing frequency of extreme events, including floods and droughts. A series of studies have been conducted, including the “Schéma Stratégique Vesdre” (Service public de Wallonie, 2024a), which establishes a vision of resilience and solidarity for the catchment. It proposes a coherent set of actions focused on adaptation to climate change, addressing both floods and droughts. In this study, these actions, tailored to the predominant land uses in the region, are tested through hydrological simulations of selected sub-catchments.

Figure 1 Study area with maps of land use/land cover.

Two sub-catchments with contrasting land uses belonging to the Vesdre catchment in Belgium were selected (Fig. 1). The first catchment (C1) is primarily agricultural (69 %) and peri-urban (17 %), covering 39.9 km² with altitudes ranging from 125 m to 325 m. The climate is temperate oceanic with an annual precipitation rate over the catchment of 952 ± 96 mm (mean and standard deviation calculated between 2004 and 2021). The second catchment (C2) is dominated by forests (45 %), meadows (35 %) and peatlands (10 %). It spans 71 km² with elevations ranging from 350 to 687 m. The terrain varies from relatively flat upstream on the High Fens Plateau (dominated by peatlands, wet moorlands and drained plantations of coniferous trees) to deep valleys downstream (dominated by coniferous plantations, deciduous forests and meadows). Since the latter half of the 19th century, extensive surface drainage efforts have been undertaken on the upstream plateau to facilitate the establishment of conifer plantations (mainly Picea abies). This has resulted in the degradation of the region's open wetlands and peat soils. The annual precipitation rate over the catchment is 1158 ± 123 mm (mean and standard deviation calculated between 2004 and 2021).

2.1.1 Soil and geological setting

The Vesdre catchment is characterized by a complex geological history, spanning from the Cambrian to the Holocene, marked by two orogenic cycles (Caledonian and Variscan) (Fig. 2). This evolution has resulted in significant lithological and pedological diversity. The nature of the soils in the studied catchments is shaped by two main contributions.

A large portion of the soils is developed on loess deposited during various glaciation phases, particularly the Würm. These loess deposits, rich in silt, have covered the entire region, especially on plateaus and areas with moderate relief. On plateaus and gentle slopes, loess is often preserved, forming a homogeneous cover that promotes good permeability and, consequently, effective drainage. These silts, relatively uniform and free of coarse fragments, are still found today in catchment 1 (C1), where they generally reach a thickness of 1 to 2 m. In contrast, in catchment 2 (C2), these deposits persist only in isolated areas and are generally associated with weathering residues from the underlying bedrock.

On steeper slopes, erosion tends to displace these deposits, reducing their thickness or even removing them completely in certain locations. In valley bottoms, colluvium and alluvium represent the parent material for most soils.

In parallel, the degradation and weathering of the underlying rocks contribute to soil formation. The resulting fragments mix with loess deposits and often reflect the lithological nature of the bedrock. The intensity of this process varies depending on topography and rock resistance. In areas where loess is preserved, such as plateaus and moderate slopes, the contribution of rock fragments tends to be less significant, resulting in soils that are less stony. Conversely, on steep slopes and in areas where loess has been largely eroded, the influence of the underlying rock becomes more pronounced, leading to soils with a higher coarse fraction and smaller depths of loose material.

Figure 2 Geological maps of C1 and C2, showing the absolute age of geological units (period, epoch or age) and the principal rock types within these units (Service public de Wallonie, 2025a, 2025b). Soil maps of C1 and C2, depicting textural classes, natural drainage characteristics and stoniness (Service public de Wallonie, 2023b).

Soil natural drainage primarily depends on several factors. The permeability of the substrate and the nature of the materials play a crucial role. Soils derived from recent deposits, such as colluvium, alluvium and loess, are generally permeable, facilitating infiltration. However, soils formed from more altered or clay-rich units, such as the Cretaceous Herve smectites (locally present in C1), exhibit reduced permeability, leading to slower drainage.

The presence of impermeable layers also influences drainage. In certain areas, particularly upstream of C2, alteration clays or compacted layers act as barriers, restricting percolation and leading to the formation of perched water tables. Under these conditions, the decomposition of fresh organic litter slows significantly, resulting in a peat-forming soil. If this process continues, it can eventually lead to the formation of sphagnum peat, as observed in the extreme eastern part of C2. Conversely, when soils lie directly on a more permeable substrate, drainage is improved. This is notably the case in the lower part of C2, where impermeable layers tend to be thinner and soils directly overlying schists benefit from the rock's schistosity, which enhances infiltration. Similarly, in the lower part of C1, limestone formations may also contribute to efficient drainage, despite the fact that the loess overlying these limestones tends to be more clay-rich.

Topography further influences drainage conditions. In valley bottoms, the water table is often permanent and reaches the surface, limiting infiltration in these soils. Locally, peat soils can be found in the valley bottoms of C2. In contrast, on plateaus and gentle slopes covered with loess, such as in C1, drainage is generally good, except in certain areas of the upper part of C2, where clay accumulations or the presence of impermeable layers may restrict infiltration.

2.2 Modelling framework

2.2.1 Hydrological model

The model selection was based on its capability to represent long time series of multiple hydrological processes (channel flow, overland flow, unsaturated and saturated flow, evapotranspiration, etc.), accounting for interactions between the saturated zone, vadose and surface, at the catchment scale using a fully distributed 3D lattice. Therefore, the coupled hydrologic/hydraulic MIKE SHE/MIKE 1D model was chosen. MIKE SHE/MIKE 1D solves partial differential physically based equations describing mass and energy transfers. The actual evapotranspiration is derived from the reference evapotranspiration using the Kristensen and Jensen method. The 2D overland flow is calculated using the diffusive wave approximation of the Saint-Venant equations. The 1D channel flow is derived from the fully dynamic higher-order approximation of the Saint-Venant equations. The 1D unsaturated flow is derived from the Richards equation. Saturated flow is described by the 3D Darcy's equation, numerically solved using the finite difference method (DHI, 2024). All input data are summarized in Table 1.

Table 1 Summary of model input data.

Hydrometeorological data: the climate data comprised a 5×5 km grid detailing hourly precipitation intensities and daily reference evapotranspiration. Reference evapotranspiration was estimated using the Penman–Monteith equation (Allan et al., 1998) applied to data collected by the observation network of the Belgian Royal Meteorological Institute (RMI), including temperature, humidity, wind speed and solar radiation (Royal Meteorological Institute of Belgium, 2010). Precipitations were determined by combining RMI's gridded daily data with hourly point measurements from the Wallonia Public Service's rain gauge network (https://hydrometrie.wallonie.be, last access: 24 March 2023), enabling the hourly distribution of RMI's daily rainfall data.

Hourly stream discharge data were obtained from 11 May 2011 at the outlet of C1 (gauging station 1), from 30 November 2015 at the outlet of C2 (gauging station 2) and from 1 January 2000 in an upstream branch of C2 (gauging station 3) with a drained area of 19.6 km² (https://hydrometrie.wallonie.be, last access: 6 October 2022).

Topographic data: the topographic data were a hydrologically corrected digital elevation model (DEM) at a 2 m resolution (Service public de Wallonie, 2023a). Resampling onto the model grid involves selecting the minimum elevation value to maintain the hydrological correction applied.

Land cover–land use data: land use/land cover (LULC) information was used to define spatially variable surface roughness and vegetation. The LULC map was a 5 m resolution layer of the land cover of 2018 (Service public de Wallonie, 2024b; Bassine et al., 2020) mapping several LULC classes. Resampling onto the model grid was performed by selecting the predominant LULC within each grid cell. Each LULC class was assigned an M value of Manning.

The surfaces were also categorized into seven vegetation classes based on the LULC map. Apart from permanent grasslands, which were considered open production surfaces, most croplands were cultivated with maize. Consequently, maize development dynamics were assigned to all croplands, and in the absence of information on cropping practices, no crop rotations were modelled. Each vegetation class was assigned annual dynamics (temporally varying) of the leaf area index (LAI), crop coefficient factor (K_c) and root depth (R_d) (Appendix Fig. A1). These parameters are used by the model to refine the reference evapotranspiration (Penman–Monteith) into the actual evapotranspiration following the approach of Kristensen and Jensen (DHI, 2024). In addition, detention storage, which is a parameter aimed at accounting for ponding in depressions at sub-cell scale, was fixed at 4 mm in both catchments.

Table 2 Equivalences between land use/land cover classes and vegetation classes.

Soil and subsoil data: the 1D Richards equation was used to model vertical water fluxes in the unsaturated zone, which required the determination of the soil hydrodynamic properties (soil water retention and soil hydraulic conductivity curves) using the Van Genuchten and Mualem functions (Van Genuchten, 1980) at each location within the unsaturated zone. Saturated flow was modelled with the 3D Darcy's equation requiring vertical and horizontal saturated conductivities to be specified at each location within the saturated zone.

The soil and subsoil were represented as a series of layers with homogeneous hydrodynamic properties within each layer but varying between layers in terms of thickness and lateral arrangement. Except for peat and degraded peat soils, the upper soil profile was discretized into two layers: topsoil, ranging from 0 to 0.4 m, and subsoil, extending from 0.4 m to a maximum of 2 m deep. These two layers could be truncated by the base at depths of 0.2, 0.3 or 0.6 m, considering the soil depth information from the Belgian soil map (Service public de Wallonie, 2023b).

The upper soil profile was delineated into homogeneous units based on the Belgian soil map (Fig. 1) (Service public de Wallonie, 2023b). This map is derived from approximately 6 million samplings conducted with a soil auger (max 1.25 m depth) between 1947 and 1991 of the Belgian cadastral plans, according to a square grid of 75 m per side (1 to 2.5 observations per hectare). This map delineates 6000 soil units defined by characteristics such as texture, natural drainage, diagnostic horizon and stone content (Legrain et al., 2011). For each homogeneous unit, except for peat, degraded peat and impermeable soils, retention and hydraulic conductivity curves (Van Genuchten and Mualem models, with m= 1–$\mathrm{1}/n$ and L= 0.5) for both soil layers were derived from the European hydraulic pedotransfer function (Weber et al., 2020) number 1 using depth and the averaged soil texture (“Textures et fractions granulométriques de référence des sols de Wallonie – Série”, Service public de Wallonie, 2023c) as predictors (Szabó et al., 2021). The predicted saturated water content, θ_s, was subsequently corrected to account for the loss of porosity due to stone content, s, as represented in the Belgian soil map.

$$\begin{array}{}\text{(1)}& {\mathit{\theta }}_{\text{s-pebble}}={\mathit{\theta }}_{\mathrm{s}}(\mathrm{1}-s)\end{array}$$

For the peat soil units, retention curves were obtained by measurements of water content on peat samples of local bogs using pressure plates at pF values of 1.00, 1.60, 1.85, 2.00, 2.30, 2.78, 2.95, 3.62 and 4.14. Van Genuchten functions were then adjusted to measure the retention curves. The saturated hydraulic conductivity of peat is spatially highly variable and challenging to measure in laboratory. Saturated hydraulic conductivity was initially determined based on values found in the literature (Wastiaux, 2008). As suggested by Szabó et al. (2024), it was subsequently refined through calibration using discharge measurements conducted between 2012 and 2015 at the outlet (gauging station 4) of a small peatland catchment of 14 ha, located 5 km away from C2. The thickness of the degraded peat soil unit was assumed to be 0.4 m, while the peat soil unit was assumed to be 2 m thick. Beneath peat or degraded peat, a low-permeability layer of 0.6 m thickness with 60 % clay, 35 % silt and 5 % sand was assumed, and its hydrodynamic parameters were retrieved using euptfv2.

Below the upper soil profile, multiple geological units extending to a depth of 18 m were delineated based on a 1×1 km grid. Each geological unit consisted of multiple layers, exhibiting variations in thickness and vertical hydraulic conductivity. These values were defined through a bibliographic study of the hydrogeological context of the Vesdre catchment (Sohier, 2011). Due to their imprecise nature and large spatial variability, vertical hydraulic conductivities were initially defined as ranges of plausible values and were subsequently refined during the model calibration process. Horizontal hydraulic conductivities were assumed to be 10 times greater than vertical hydraulic conductivities (DHI, 2024).

Model specification, calibration and validation: for both catchments, a model of the present situation (further referred to as BASELINE) was calibrated and validated. The total simulation period was from 1 January 2002 to 31 December 2021. At the start of the simulation process, the soil water content was set to field capacity. The first 2 years of simulation were considered a warming period and were not considered in the data analysis. Outputs were stored on an hourly basis. To ensure a fine representation of the catchments while maintaining reasonable computational times, the grid resolution of the C1 model was 20 m, whereas it was 40 m for the C2 model. Calibration was conducted in the period between 1 October 2017 and 1 October 2021, while validation encompassed the remaining available data.

The calibration was performed manually and limited to very few parameters, for which the available data were considered uncertain and had the most impact (the most sensitive parameters) on the modelled hydrographs. These parameters included the Manning's M roughness coefficients for overland and channel flows (which are mostly sensitive regarding the temporality of hydrographs and peak flows), as well as the saturated hydraulic conductivity of the geological layers (which are sensitive regarding baseflow, peak discharge values and groundwater head) (XEVI et al., 1997). Calibration of the Manning's M roughness and the saturated hydraulic conductivities of peat was performed using discharge measurements conducted between 2012 and 2015 at the outlet of the small peatland catchment of 14 hectares (gauging station 4). Calibration and validation were performed in relation to historical discharge measurements taken at the river gauging stations. The ability of the model to reproduce observed discharges was assessed against the Moriasi et al. (2007) model evaluation guidelines. This included the computation of the Nash–Sutcliffe efficiency (NSE), the ratio of the root mean square error to the standard deviation of measured data (RSR) and the percent bias (PBIAS).

$$\begin{array}{}\text{(2)}& \mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{{\sum }_i({{\mathrm{Q}}_{\mathrm{Obs}}}_i-{{\mathrm{Q}}_{\mathrm{Mod}}}_i{)}^{\mathrm{2}}}{{\sum }_i({{\mathrm{Q}}_{\mathrm{Obs}}}_i-\frac{{\sum }_i{{\mathrm{Q}}_{\mathrm{Obs}}}_i}{n}{)}^{\mathrm{2}}}}\end{array}$$

NSE varies from -∞ to 1. The closer to 1, the better the simulation performance. Q_Obs and Q_Mod are the observed and simulated discharges, respectively. i is the time step.

$$\begin{array}{}\text{(3)}& \mathrm{RSR}={\displaystyle \frac{\mathrm{RMSE}}{{\mathrm{STDEV}}_{\mathrm{obs}}}}={\displaystyle \frac{\sqrt{{\sum }_i({{\mathrm{Q}}_{\mathrm{Obs}}}_i-{{\mathrm{Q}}_{\mathrm{Mod}}}_i{)}^{\mathrm{2}}}}{\sqrt{{\sum }_i({{\mathrm{Q}}_{\mathrm{Obs}}}_i-\frac{{\sum }_t{{\mathrm{Q}}_{\mathrm{Obs}}}_i}{n}{)}^{\mathrm{2}}}}}\end{array}$$

RSR varies from 0 to +∞. The closer to 0, the better the simulation performance.

$$\begin{array}{}\text{(4)}& \mathrm{PBIAS}={\displaystyle \frac{{\sum }_i({{\mathrm{Q}}_{\mathrm{Obs}}}_i-{{\mathrm{Q}}_{\mathrm{Mod}}}_i)\cdot \mathrm{100}}{{\sum }_t{{\mathrm{Q}}_{\mathrm{Obs}}}_i}}\end{array}$$

PBIAS varies from −∞ to +∞. The closer to 0, the better the simulation performance. Positive values indicate model underestimation bias, and negative values indicate model overestimation bias.

In the absence of piezometric data at the study sites, a graphical comparison (Fig. 4) was conducted between the simulated vertical dynamics of the saturated zone (from October 2017 to October 2020) at each grid cell with the corresponding drainage class from the Belgian soil map (Table 3). This ensured that the modelled saturated zone dynamics match those described by the soil map, which was derived from field observations. Drainage class is a concept used in different soil classifications to describe the frequency and duration of soil wetness under natural conditions (Ditzler, 1999). In the Belgian soil classification, drainage classes are determined by the depth at which mottling and reduction processes occur in the soil profile, often corresponding to the levels of saturated zone fluctuations (Table 3). This method facilitates semi-quantitative spatial validation across all areas of the catchments.

Table 3 Definition of the drainage class in the Belgian soil classification for loam and clay soils.

2.2.2 Landscape planning scenarios

The impact of NbSs on catchment hydrology was evaluated through the comparison of responses between BASELINE and POST configurations, under consistent meteorological forcing. The POST landscape configuration represented a scenario where soil–vegetation NbSs were implemented at the scale of the entire watershed. These landscape planning POST scenarios were developed in accordance with the “Schéma Stratégique Vesdre” (Service public de Wallonie, 2024a), which outlines a long-term vision for sustainable territorial development and proposes a range of potential measures to mitigate flood and drought risks. The modelled scenarios (POST) incorporate as many feasible measures as possible within specific contexts (agriculture-dominated and forest-dominated) to assess their potential for flood and drought mitigation. The scenarios compare the current situation to a long-term projection (2050 horizon), assuming full implementation and maturity of the measures. The transition process between current and future states is not modelled, although it may significantly influence the actual effectiveness over time. It should be emphasized that accounting for the transition phase may be essential, as some measures require several years or even decades to achieve full efficiency, and the speed at which a measure becomes effective may serve as an important criterion for prioritizing NbSs.

Figure 3 Modelled NbSs in POST configurations of C1 and C2.

The following NbS measures were implemented in the agricultural catchment (C1) (Fig. 3a):

• Hedgerows: a dense network of hedgerows was systematically implemented around all agricultural parcels, representing around 700 km of barriers within the entire watershed (≈ 175 m ha⁻¹) (Table 4a).

• Agricultural practices: water and soil conservation practices were modelled on cultivated land, including soil pitting for maize crops (Table 4c) to form small depressions between rows (0.73 km², or 1.8 %) (Clement et al., 2023) and the adoption of reduced tillage practices (Table 4b) for the other crops (0.30 km², or 0.8 %).

In the forest catchment (C2), the following NbS measures were implemented (Fig. 3b):

• Restoration of peatlands: this involved restoring peatlands and open wetlands on degraded peaty soils, currently occupied by conifer plantations and moors (4.6 km², or 6.5 %). The scenario simulated the plugging of the existing surface drainage network with 275 plugs distributed along concentrated flow paths (drainage area from 2 to 5 ha) and the creation of 32 ponds of 640 m³ each (Table 4f).

• Forest diversification: this scenario was complemented by forest diversification, involving the conversion of monospecific conifer plantations into irregular mixed stands on temporarily waterlogged soils (5.4 km², or 7.6 %) (Table 4e).

• Practices aimed at limiting soil compaction: for all forest areas, the potential effect of forest practices aimed at limiting soil compaction, notably through compartmentalization with designated skid trails, was incorporated into this scenario (16.8 km², or 23.7 %) (Table 4d).

These NbS measures are incorporated into the model by adjusting some key parameter values (Table 4). It should be noted that although some literature exists on how NbSs may influence these parameters, significant uncertainties remain regarding their exact values, as they cannot be assumed to apply directly to this specific case. These parameter values should be considered as hypotheses, while remaining within plausible ranges.

Table 4 Summary of hypotheses for parameter modifications in models to represent NbS scenarios.

2.2.3 Integrated hydrological analysis of model outputs

A complete processing chain of results focusing on the quantification of the hydrologic impact of NbS measures and scenarios was conducted (Table 5). Hydrologic indicators were calculated or developed to capture the impact of various measures within the watershed, emphasizing the spatial variability of their effectiveness. These indicators are computed to assess how NbSs could influence hydrological extremes, particularly focusing on flooding events and agricultural droughts.

Table 5 Summary of hydrological indicators aimed at assessing the effectiveness of NbSs in mitigating floods and drought.

Peak flow and flood volume percent reduction: these indicators were integrative metrics permitting the evaluation of the overall effect of a given scenario on flood reduction without discriminating against the impact of the various measures that comprise the scenario. To this end, these two indicators were calculated for a set of rainfall events generating a detectable signal on the hydrograph. A focus was also made on the rainfall event between 13 and 18 July 2021:

Peak flow percent reduction – the evolution of peak discharges between the BASELINE and the POST scenarios.

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \mathrm{Peak}\mathrm{flow}\mathrm{percent}\mathrm{reduction}\\ & =\mathrm{100}\cdot {\displaystyle \frac{\mathrm{Peak}{\mathrm{discharge}}_{\mathrm{POST}}-\mathrm{Peak}{\mathrm{discharge}}_{\mathrm{BASELINE}}}{\mathrm{Peak}{\mathrm{discharge}}_{\mathrm{BASELINE}}}}\end{array}\end{array}$$

Flood volume percent reduction – The evolution of the total amount of water discharged at the outlet between the BASELINE and the POST scenarios.

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& \mathrm{Flood}\mathrm{volume}\mathrm{percent}\mathrm{reduction}\\ & =\mathrm{100}\cdot {\displaystyle \frac{\mathrm{Total}{\mathrm{volume}}_{\mathrm{POST}}-\mathrm{Total}{\mathrm{volume}}_{\mathrm{BASELINE}}}{\mathrm{Total}{\mathrm{volume}}_{\mathrm{BASELINE}}}}\end{array}\end{array}$$

Spatialized indicators: these indicators aimed to assess the local hydrological impact of specific NbS measures where they were located within the watershed. This allowed for discrimination of measures based on their relative impact within the same scenario. These indicators presupposed that the hydrological effects observed at a given location are exclusively attributed to the measure implemented at that site, irrespective of the presence of other measures elsewhere within the watershed.

Spatial variation of infiltration rates: as a spatialized indicator to assess the effectiveness of NbSs in mitigating flood, the cumulative infiltration, I_cum, between 13 and 18 July 2021 was calculated for each model grid cell.

$$\begin{array}{}\text{(7)}& I_{\mathrm{cum}}={\sum }_i({\mathrm{Infiltration}}_i-{\mathrm{Seepage}}_i)\end{array}$$

Infiltration is the amount of water (in mm) flowing into the soil. Seepage is the amount of water (in mm) outflowing from the soil. i is the time step.

Spatial vegetation drought stress duration frequency: this is an indicator to assess the effect of NbSs on vegetation water stress was (newly) developed, utilizing a frequency analysis. Consecutive days where plants at specific locations underwent water stress beyond a designated threshold against a return period were mapped. Thresholds were determined based on the soil water potential within the root zone. An extreme stress threshold was defined when the water potential in the root zone drops below −150 m ($\approx -\mathrm{1.5}$ MPa), closely approaching the permanent wilting point, which is classically recognized by soil scientists as the soil water potential beyond which most plants cannot survive (Wiecheteck et al., 2020). Additionally, a moderate stress threshold for vegetation was established at a water potential of −30 m ($\approx -\mathrm{0.3}$ MPa), corresponding in our study sites to approximately 30 % of the relative extractable water content in soil (Gourlez de la Motte et al., 2020; Granier et al., 2007). The methodology to construct this indicator is briefly described as follows: for each cell of the model and for each year (2004–2021), the maximum consecutive duration during which the water potential in the soil profile explored by the roots (rhizosphere) was below the stress threshold was extracted. Then, for each cell, a Gumbel distribution was used to establish the relationship between the duration (y) and the probability of occurrence (return period: TM) of the stress.

$$\begin{array}{}\text{(8)}& y=-c\left(\ln \ln {\displaystyle \frac{\mathrm{TM}}{\mathrm{TM}-\mathrm{1}}}\right)-a\end{array}$$

c and a are adjustment parameters. TM corresponds to the return period and is calculated according as

$$\begin{array}{}\text{(9)}& \mathrm{TM}={\displaystyle \frac{n+\mathrm{1}}{m}}.\end{array}$$

n is the rank in descending order, and m is the total number of years of observations. Once the parameters a and c are fitted for a given stress threshold, y can be retrieved using Eq. (8) (Maidment, 1992).

3 Results and discussion

3.1 Calibration and validation of the baseline model

The simulated discharges evaluated at each gauging station, against Moriasi criteria, showed satisfactory to good performance (Table 6). Among the three indicators, the NSE values were the least satisfactory according to Moriasi's criteria. This could be explained by the fact that NSE is very sensitive to peak flows, whereas measured flow rates during these events are more uncertain (well above the maximum gauged discharge on the rating curve). RSR and PBIAS were considered satisfactory/good (0.50 < RSR < 0.60 for good, 0.60 < RSR < 0.70 for satisfactory) and very good (PBIAS < ± 10).

Table 6 Summary of model performance at gauging stations.

The modelled dynamics of the saturated zone depth agreed well with the expected dynamics for each drainage class on the Belgian soil map (Fig. 4). The range of variation, as well as the frequencies of the modelled saturated zone depths, were slightly higher in the degraded peat soils compared with the intact peat (Wastiaux, 2008). Soils designated as “without waterlogging” on the soil map were seldom, if ever, subject to waterlogging. In most grid cells, the modelled saturated zone remained well below the surface, with few exceptions in late winter. As expected in temporarily waterlogged soils, the modelled saturated zone often reached the surface during the winter months and receded during the summer. Overall, the model successfully represents the variability in the natural drainage characteristics of soils in the BASELINE scenario, thus addressing our first objective, which was a prerequisite for assessing the effectiveness of the NbS measures on various soils with differing natural drainage characteristics.

Figure 4 Modelled saturated zone depth (watersheds C1 and C2) for each drainage class from the Belgian soil map. The thick line represents the median values calculated from each model grid cell for each drainage class. The thin lines indicate the 25th and 75th percentiles. The stacked bar chart represents the relative proportion of the surface of the natural drainage classes of soils found in C1 and C2.

3.2 Effectiveness of NbS scenarios against floods at the outlets

The modelled specific peak discharges under the BASELINE scenario (Fig. 5) indicate that, for an equal surface area, C2 generates more runoff than C1, despite being predominantly forested, a land cover generally associated with lower runoff generation compared to pasture (Bathurst et al., 2020) or urbanized areas. This result may seem counterintuitive, as it does not align with the common assumption that forest-dominated catchments produce less runoff than agriculturally dominated ones. However, C1 and C2 differ in terms of topography, morphometry, precipitation rate, geology and soil characteristics (Figs. 1 and 2), all of which can influence runoff production and transfer within a catchment. Figure 4 also highlights that C2 has, compared to C1, a higher proportion of waterlogged soils, especially in the upstream part (Fig. 2), which hinder infiltration and promote runoff.

Figure 5 Modelled specific peak discharges at the outlets of C1 and C2 before and after the implementation of landscape planning scenarios. The top section shows scatter plots illustrating the evolution of peak discharges for a large number of runoff events between 2004 and 2021. The regression slope coefficients indicate the average peak flow percent reduction. The squared dot highlights the peak discharge of the event occurring between 13 and 18 July 2021. The bottom section presents hydrographs of this specific event, comparing conditions before and after the implementation of the landscape planning scenario.

The extreme rainfall event between 13 and 18 July 2021 was one of the events generating the highest peak discharge (represented by the squares in the upper part of Fig. 5 in the modelled time series). The modelled hydrographs of this event in both catchments are presented in the lower part of Fig. 5.

Landscape planning scenarios showed contrasting modelled effectiveness in terms of reducing peak discharges between the two catchments. In C1 (mostly agricultural), modelled peak discharges were reduced by approximately 30 %, while in C2 (forest-dominated), the reduction was approximately 10 %. Moreover, for both catchments, the effectiveness in reducing peak discharges seemed not influenced by the magnitude of the events, as indicated by the linear trend observed in the scatter plot between the BASELINE and POST peak discharges. The response remained relatively consistent for both low- and high-discharge events. It is important to acknowledge that these results reflect the overall outcomes of two distinct scenarios applied in catchments that differ in land use, geology, topography and precipitation. They should not be interpreted as a generalizable measure of the effectiveness of the specific interventions implemented within each scenario.

3.3 Effectiveness of individual NbSs in mitigating floods and agricultural droughts

The following sections explore the effectiveness of individual measures. To this end, the analysis focuses on the rainfall event that occurred between 13 and 18 July 2021.

As highlighted in Fig. 6, model outcomes indicate that NbSs generally exhibit positive effects on infiltration (and runoff reduction) and the mitigation of agricultural drought. However, the magnitude of these effects seems to vary depending on the specific type of NbS and context, as observed in the previous section.

Figure 6 Variation of cumulative infiltration (from 13 to 18 July 2021) and of vegetation drought stress duration (matric potential of −30 m and return period of 25 years) on model grid cells before (BASELINE) and after (POST) implementing NbSs (a: hedgerows; b: reduced tillage; c: soil pitting; d: forest practices aimed at limiting soil compaction; e: forest diversification with practices aimed at limiting soil compaction; and f: restoration of peatlands).

3.3.1 Hedgerows

In the agricultural context (Fig. 6a), model results suggest that hedgerows can be highly effective in improving infiltration during an extreme rainfall event, as reviewed by Rajbanshi et al. (2023). This modelled increase in infiltration during extreme events could be attributed to two primary factors: the enhancement of the soil's hydrodynamic properties and the generally drier soil conditions under hedgerows, which increase the unfilled soil porosity and the hydraulic gradient responsible for infiltration. Due to the lack of farm traffic, the greater incorporation of organic matter or a root network promoting macropores, soils under hedgerows are often found to be less compacted and more permeable than those in adjacent arable or pasture fields (Holden et al., 2019; Wallace et al., 2021). Drier soil conditions under hedgerows could be attributed to increased actual evapotranspiration (Benhamou et al., 2013) and canopy interception. A counterintuitive observation from our simulations is that, when considering a broader temporal window, the cumulative annual modelled infiltration decreased under hedgerows due to this increased canopy interception; a larger portion of the rainfall was directly evaporated before reaching the soil for smaller and more common events (results not shown). Hedgerows have also been modelled to reduce agricultural droughts, as their root systems extend deeper than those of crops or meadows (Table 2). This allows them to access water in deeper and wetter soils layers even though soils in the surface were generally drier.

3.3.2 Agricultural practices

Soil pitting in maize crop also appeared to enhance modelled infiltration during extreme rainfall events (Fig. 6c), as observed by Clement et al. (2023). However, model outcomes showed a minimal effect of soil pitting on agricultural drought mitigation. This may suggest that extreme rainfall events, when a surplus of infiltration is observed, may be temporally disconnected from droughts events; hence, they do not follow close enough together in time for the excess infiltrated water to remain in the root zone until the agricultural drought begins. Another potential factor could be the implementation of hedgerows along field margins in the same scenario, drying soils at the edges and negatively affecting water availability for adjacent crops (Forman and Baudry, 1984).

For reduced tillage practices (Fig. 6b), we did not model large changes in infiltration during extreme events, as reviewed by Clement et al. (2024), nor did we model mitigation of agricultural drought. A well-known positive effect of conservation tillage is the concentration of organic matter into the uppermost soil layer (Haddaway et al., 2017), improving, in turn, aggregate formation (Chen et al., 2020; Obalum et al., 2019), pore size distribution and soil water retention characteristics (Chandrasekhar et al., 2019). In our model scenario, reduced tillage was virtually implemented by reducing the Van Genuchten alpha parameter of the topsoil by 20 %. This adjustment slightly shifted the soil porosity towards smaller pores, specifically from the range of 0 to −330 hPa to the range of −330 to −15 000 hPa, consequently increasing the soil's available water holding capacity. Despite expectations that this improvement in soil physical quality would buffer agricultural drought, our simulation results showed no significant positive effect on mitigating it. Wittwer et al. (2023) recently suggested that while conservation tillage practices could improve soil characteristics, they do not significantly affect root water uptake patterns under severe drought conditions. Moreover, they emphasized that plant ecophysiology and cropping system diversification in space (crop mixture) and time (crop rotation) (Sanford et al., 2021) have a greater impact on crop performance under drought than soil quality. This is a conclusion that our simulation results seem to support. NbSs that were modelled to solely improve the retention characteristics of surface soils without affecting root depth, such as reduced tillage practices or even forest practices aimed at limiting soil compaction, showed a minimal effect on reducing vegetation water stress. One possibility is that, during dry periods, soil desiccates from the surface downward over time. Therefore, if the plant root system extends deeper than the surface soil layer, enhancing the retention properties of this superficial layer, where desiccation occurs first, does not seem to be very impactful for increasing plant resilience to prolonged droughts.

3.3.3 Forest practices

In the forest context, model results indicate that forest diversification could be effective in improving infiltration during extreme rainfall events and in decreasing days of agricultural droughts (Fig. 6e). Mixed forests are often shown to outperform monocultures regarding agricultural drought resistance, especially in dry areas (Liu et al., 2022; Pardos et al., 2021). In our model scenario, forest diversification was implemented primarily by increasing rooting depth, reflecting the observation that mixed forests often develop deeper and more structurally developed root systems than monocultures, especially spruce even-aged plantations that are common in the C2. This results from the inherent differences in rooting patterns between species, leading to generally deeper root extensions or denser fine root systems in mixed stands compared to pure stands (Jose et al., 2006; Meinen et al., 2009), as well as from the competitive environment in mixed stands, which may encourage greater root investment (Reubens et al., 2007). However, the drought resilience of mixed forests may also be attributed to other mechanisms, such as resource partitioning (e.g. root stratification), facilitation (e.g. active hydraulic redistribution), and selection (e.g. a higher likelihood of including drought-tolerant species) (Pardos et al., 2021). However, these complex mechanisms were not explicitly represented in the model.

Practices aimed at limiting forest soil compaction (Fig. 6d) also appeared to be effective in improving infiltration but exhibited low efficiency in reducing agricultural drought, as rooting depth was not increased with this measure.

In accordance with the current literature (Wastiaux, 2008), our model showed almost no improvement of infiltration and agricultural drought mitigation after peatland restoration in moors or in drained spruce plantations on degraded peat soil (Fig. 6f).

3.4 Spatial variability of effectiveness of NbSs against flood

This section examines spatial variability and the role of soil characteristics in the effectiveness of modelled NbSs against floods, with the aim of identifying the most effective NbS–location combinations. The displayed cumulative infiltration data from 13 to 18 July 2021, before (BASELINE) and after (POST) implementing NbSs, are modelled for each natural drainage class (Table 3) from the Belgian soil map (Fig. 7). Almost all infiltration distributions exhibited a characteristic bimodal shape: one peak near zero, where the soils do not or barely infiltrate, and another peak where infiltration approaches the total cumulative rainfall. Based on these distributions, we can identify runoff production zones (below the cumulative precipitation) and runoff interception/re-infiltration zones (above). Grid cells with negative infiltration rates indicate areas where seepage flow is occurring.

Figure 7 Comparison of cumulative infiltration (from 13 to 18 July 2021) before (BASELINE) and after (POST) implementing NbSs (a: hedgerows; b: reduced tillage; c: soil pitting; d: forest practices aimed at limiting soil compaction; e: forest diversification with practices aimed at limiting soil compaction; and f: restoration of peatlands) for each drainage class from the Belgian soil map. Each data point (n) represents the infiltration of one model grid cell. Significance levels – ns (p-value > 0.05), ^* (0.01 < p-value < 0.05), ${}^{**}$ (0.001 < p-value < 0.01) and ${}^{***}$ (p-value < 0.001) – indicate the results of Kolmogorov–Smirnov tests comparing BASELINE and POST distributions based on 200 randomly sampled data points. Significance levels are omitted for sample sizes smaller than 200. The centred vertical black lines depict the difference in kernel density estimates between different BASELINE and POST distributions. The numbers below each boxplot and black point are the means. Horizontal dotted lines represent zero infiltration and the mean cumulated rainfall between 13 and 18 July 2021.

Observing the density distributions of infiltration, a shift is evident between non-saturated soils and more saturated soils, with soils without waterlogging generally exhibiting higher infiltration rates than waterlogged soils.

When comparing the BASELINE and POST infiltration distributions and their significance levels, it becomes clear that the effectiveness of NbSs in enhancing infiltration diminishes as soil water saturation increases. This trend is observed across a gradient from well-drained soils to temporarily waterlogged areas, degraded peat soils and permanently waterlogged zones. It should be acknowledged that other correlated local factors, such as topography, may also contribute to the observed pattern and could influence NbS effectiveness. However, these findings still suggest that NbSs located on well-drained soils may offer greater potential for improving infiltration, implying that natural soil drainage characteristics should be considered when prioritizing areas of NbS implementation.

This trend might also explain why forest practices aimed at reducing soil compaction have been slightly more effective than forest diversification combined with these practices (Fig. 6); the latter were mostly implemented on soils prone to temporary waterlogging, whereas the former were predominantly applied to well-drained soils (Fig. 7). In waterlogged soils, the enhanced permeability at the surface did not result in increased infiltration rates, probably because the soil was saturated or close to saturation.

Conversely, in the absence of waterlogging constraints, enhancing soil surface properties seems to significantly increase infiltration. This effect appeared to be further amplified by changes in initial soil moisture conditions before rainfall events, as in the cases of hedgerows and forest diversification (Fig. 8). The transition to more evaporative vegetation generally leads to drier surface soils locally. This drying reduces the shallow soil water matric potential, thereby increasing the hydraulic gradient responsible for infiltration. However, the magnitude of this effect on the variation of initial moisture conditions is not uniformly distributed within the watershed and seems to vary with soils' natural drainage or the average saturated zone depth (Fig. 8).

In areas where soils are permanently or temporarily waterlogged, the implementation of hedgerows or forest diversification did not affect initial soil moisture conditions or had only a minimal effect. This may be due to the anoxic conditions that might limit root development and transpiration (Caubel, 2001). Additionally, in these areas, increased evapotranspiration may not result in drier soils, as water consistently flows towards these areas from adjacent waterlogged soils.

In initially dry soils, increased evapotranspiration results in somewhat drier soil conditions, but the effect on enhanced infiltration is moderate. The reason for this is that the modelled infiltration capacity of these soils was barely surpassed by rainfall intensity on the extreme event of July 2021 (Fig. 8).

Figure 8 Impact of hedgerows and forest diversification on infiltration and antecedent soil moisture conditions as a function of the average saturated zone depth. (a, b) Cumulative infiltration from 13 to 18 July 2021. The horizontal line is the spatially averaged cumulated rainfall. (c, d) Variation of cumulative infiltration from 13 to 18 July 2021. (e, f) Variation of unfilled macroporosity (porosity defined between 0 and −100 cm of the matric potential) in the top 30 cm of soil on 12 July 2021. The abscissa is the average modelled saturated zone depth (from 1 January 2004 to 31 December 2021). Data points represent average observations for fixed saturated zone depth intervals. The dotted line represents a non-linear bell-shaped regression (Bragg: $y=c+\left(d{\exp }^{\left(-b{\left(x-e\right)}^{\mathrm{2}}\right)}\right))$ fitted to observed data, where b is the slope around the inflection points of the bell, c is the constant adjusting the minimum, d+c is the maximum value at the top of the bell and e is the x value at which such maximum occurs.

In contrast, soils with an average groundwater table depth between 0.5 and 2 m seem to exhibit optimal conditions for improving infiltration, as evidenced by a significant bell-shaped function (Fig. 8). In these areas, increased evaporation enhances unfilled soil macroporosity, also represented by a bell-shaped function (Fig. 8). The resulting increased pore space, coupled with a higher hydraulic gradient, may significantly boost infiltration. These soils corresponded to soils without waterlogging with moderate or imperfect drainage in the Belgian soil classification. Our findings suggest that such areas may be the most favourable for implementing NbSs that enhance evapotranspiration, such as hedgerows, to improve infiltration within a watershed.

Additionally, in soils without waterlogging, re-infiltration (also called runon) of runoff generated uphill can occur. This phenomenon is clearly illustrated in Fig. 7, where infiltration rates exceed precipitation rates. Typically, runoff–runon processes occur where areas of low infiltration capacity lie upslope of areas with higher infiltration capacity (Woolhiser et al., 1996). Hence, runoff–runon processes highly depend on the spatial heterogeneity of soil infiltration capacity. In our simulations, runoff–runon processes appeared to mostly coincide with areas with good infiltration capacity just downstream of impermeabilized soils, such as roads or urbanized zones. Further research would be valuable in determining whether these areas should be prioritized for implementing NbSs aimed at retaining overland runoff and enhancing re-infiltration.

3.5 Spatial variability of effectiveness of NbSs against agricultural drought

The density distributions of days experiencing vegetation stress across different soil drainage classes reveal that such stresses are most prevalent in soils without waterlogging (Fig. 9). The vegetation in temporarily waterlogged areas, degraded peat soils and permanently waterlogged zones was far less susceptible to these stresses, primarily due to the naturally higher soil water availability in these environments, even after prolonged droughts. Therefore, NbSs installed on non-waterlogged/water-limited soils were the most effective in limiting agricultural drought. As a concrete example, hedgerows implemented on temporarily waterlogged soils showed a limited effect on drought compared to soils without waterlogging (Fig. 9a).

Figure 9 Comparison of days of drought (matric potential of −30 m and return period of 25 years) before (BASELINE) and after (POST) implementing NbSs (a: hedgerows; b: reduced tillage; c: soil pitting; d: forest practices aimed at limiting soil compaction; e: forest diversification with practices aimed at limiting soil compaction; and f: restoration of peatlands) for each drainage classification from the Belgian soil map. Each data point (n) represents the infiltration of one model grid cell. Significance levels – ns (p-value > 0.05), ^* (0.01 < p-value < 0.05), ${}^{**}$ (0.001 < p-value < 0.01) and ${}^{***}$ (p-value < 0.001) – indicate the results of Kolmogorov–Smirnov tests comparing BASELINE and POST distributions based on 200 randomly sampled data points. Significance levels are omitted for sample sizes smaller than 200. The centred vertical black lines depict the difference in kernel density estimates between different BASELINE and POST distributions.

3.6 Synergies and trade-offs between flood and drought mitigation with NbSs

The tested NbSs primarily had a positive impact on flood reduction, with 80.6 % of managed grid cells showing improved infiltration capacity (Fig. 10). In contrast, the positive impact on agricultural drought was less significant, with only 27.3 % of managed grid cells exhibiting improvement. However, we observed a minimal negative impact from NbSs on agricultural droughts, with only 4.1 % of managed grid cells being affected. In most cases (68.6 % of managed grid cells), there was no discernible impact of NbSs on agricultural droughts.

Figure 10 Scatter plot and regression lines of the variation cumulative infiltration (from 13 to 18 July 2021) against days of drought (matric potential of -30 m and return period of 25 years) on model grid cells before (BASELINE) and after (POST) implementing NbSs (a: hedgerows; b: reduced tillage; c: soil pitting; d: forest practices aimed at limiting soil compaction; e: forest diversification with practices aimed at limiting soil compaction; and f: restoration of peatlands). Vertical and horizontal lines represent thresholds between negative and positive impacts on flood and drought, respectively. Percentages represent the proportion of total points (model grid cells) inside and between each quadrant.

Managed grid cells that were positively affected against both floods and agricultural droughts (25.1 %) were predominantly those with hedgerows and forest diversifications (Fig. 10). Nevertheless, we believe that the observed relationship between increased infiltration during extreme rainfall events and reduced vegetation water stress is not causal, just correlated. This correlation may be primarily explained by changes in vegetation and modelled root depth (Tables 2 and 4) rather than by an increase in infiltration. Notably, NbSs that substantially reduced the number of days experiencing water stress, such as hedgerows or forest diversification combined with practices aimed at limiting soil compaction, also resulted in drier soils due to higher initial interception and evapotranspiration rates. However, a common characteristic of these NbSs is the modification of the vegetation and, consequently, root depth, counterbalancing drier upper soil layers. Increased root depth benefits against agricultural drought in two primary ways: it expands the volume of soil reachable to roots, thereby increasing the total available water, and deeper root systems enable plants to access deeper and typically wetter soil layers, which can be critical during periods of prolonged drought. It has been reviewed that a deeper root system can significantly enhance drought resistance, particularly in water-limited environments where subsoil water is available (Li et al., 2022).

On the other hand, the implementation of NbSs that extend root depth can enhance infiltration, as an extensive root network facilitates the formation of macropores within the soil matrix, even after root senescence (Beven and Germann, 1982; Wallace et al., 2021). Additionally, a robust root system supports substantial water uptake to facilitate transpiration, which consequently leads to soil drying (Caubel, 2001). Soil desiccation promotes infiltration by increasing the hydraulic gradient between surface and soil, by creating more air filled spaces for water to infiltrate, and by inducing structural cracks that facilitate rapid preferential flow (Zhu et al., 2018). However, this very specific phenomenon was not represented in the model.

From this perspective, promoting root growth and sustaining plant transpiration could be considered important to help reduce both vegetation water stress, caused by prolonged drought, and flood risks. In an agricultural context, this can be achieved by limiting soil compaction, which might also hamper infiltration (Unger and Kaspar, 1994). For instance, biopores have been found to substantially mitigate transpiration deficits during droughts by promoting root systems as preferential growth paths, permitting crops to access water from deeper and moister soil layers (Landl et al., 2019). Compaction caused by machinery traffic in forest soils may also reduce the rooting density of trees due to a lack of soil aeration (von Wilpert and Schäffer, 2006) and an increased ground penetration resistance (Goutal, 2012).

Finally, it is important to note that the balance between synergies and trade-offs is not solely dependent on the type of NbS but also on the location where they are implemented. As discussed in previous sections, NbSs were more effective in soils without waterlogging, both in mitigating the risks of flooding and addressing agricultural droughts. As a direction for future work, this opens the door to incorporating allocation change experiments and multi-objective optimization into our modelling framework. Because different NbSs or combinations of NbSs can create both benefits and trade-offs across multiple (often conflicting) objectives, it is essential to identify one or several optimal implementation strategies that maximize benefits while minimizing trade-offs. Among others, one approach is based on Pareto-optimality (Deb and Jain, 2014), where a set of optimal solutions is generated such that no objective can be improved without worsening another. From this set of optimal alternatives, decision-makers can explore and select appropriate solutions based on their priorities. Among the available tools for such analyses, the Constrained Multi-Objective Optimization of Land Use Allocation (CoMOLA) software (Strauch et al., 2019) has been developed to facilitate these assessments and has been applied to other NbS modelling frameworks using SWAT+ (European Commission, 2020).

3.7 Study limitations and knowledge gaps

3.7.1 Data requirements

Representing the hydrological functioning of a watershed at a fine spatial scale remains a challenge for modellers. Physically based distributed hydrological models offer great potential for capturing detailed processes, but they require a large number of parameters to be determined or calibrated. In many cases, the data required for model calibration and validation are unavailable or highly uncertain at the necessary scale. For instance, calibrating and validating models at the watershed outlet does not necessarily ensure accurate representation of the spatialized hydrological processes within the watershed. Multiple combinations of spatialized model input parameters can yield identical modelled hydrographs at the catchment outlet. To address these challenges, a semi-quantitative spatial validation was carried out by comparing the spatially distributed average depths of the modelled saturated zone with the drainage classes from the Belgian soil map. This exercise demonstrated that the modelling approach effectively captures the fine-scale interactions between soil and water processes. Although this method is specifically tailored to the Belgian context, as this map is specific to Belgium, it provides a valuable framework that could be adapted to other regions. This calls for future research on spatially distributed and cross-scale hydrological observations to further enhance model calibration and validation.

Our results indicated that soil infiltration capacity may be a crucial factor in determining the most suitable locations for implementing NbSs. However, accurately representing the spatial and temporal (Pirlot et al., 2024) variability of soil hydrodynamic properties remains challenging with current methodologies. The maps depicting the soils' hydrodynamic properties are inherently limited in their representativeness due to uncertainties in the applied pedotransfer functions, which are primarily derived from agricultural soils, overlooking much of the soil diversity, such as stony soils. Additional limitations arise from uncertainties and the spatial resolution of the underlying soil property maps (e.g. soil depth, texture, organic carbon and bulk density), as well as the suitability of the continuous function used to describe soil hydrodynamic properties (Weber et al., 2024). Although characterizing the spatiotemporal variability of soil hydrodynamic properties offers valuable benefits in many fields such as hydrology, agriculture, irrigation, ecology or environmental sciences, this remains poorly explored and warrants further research.

Similarly, vegetation plays a crucial role in the soil hydrological functioning at fine scales, significantly influencing watershed dynamics. A striking example is the impact of planting hedgerows, which locally decreases soil moisture, enhancing infiltration rates, resulting in substantial reductions in peak flow at the watershed outlet. While these effects highlight the importance of vegetation, key aspects of its characteristics (also used as model input), particularly tree rooting systems, remain surprisingly underdeveloped. Despite their critical role, our understanding of roots, including their growth, development and interaction with complex physical, chemical and biological soil properties, is still limited (Centenaro et al., 2018).

3.7.2 Uncertainties

Three major sources of uncertainty affect estimates of NbS effectiveness in our study: (i) the choice of the hydrological model and its ability to represent key processes, (ii) the model parameterization and modelling of hydrological functioning across the entire catchment in the BASELINE scenario, and (iii) the parameterization of NbSs (Brauman et al., 2022). The first source arises from the inherent limitations of any modelling approach, as different models emphasize different processes and may not fully capture key local dynamics. Unlike many other models, the proposed approach explicitly accounts for the spatial variability of local conditions, which may influence NbS effectiveness (Lalonde et al., 2024). However, certain key processes remain unrepresented, for instance, the recycling of atmospheric moisture (Keys et al., 2018; Theeuwen et al., 2023) and the temporal variability of soil hydraulic properties (Pirlot et al., 2024), which may introduce additional uncertainties. Furthermore, it is critical that the model's spatial resolution corresponds to the scale at which NbSs are implemented. For example, the simulation of hedgerows is limited by the model's 20 m resolution, which may not accurately capture their true configuration, introducing uncertainty. The second source is partially addressed through our model validation process, which integrates all available data, including data from gauging stations and drainage classes from the Belgian soil map, enhancing confidence in the model's representation of catchment-scale hydrological dynamics. The third source has not been explicitly quantified. Instead, we compare our results with reported NbS effectiveness in the literature to assess whether our findings align with expected trends. While this provides a useful reference, variations in study contexts may introduce further uncertainties. This calls for more field-based evaluation of these measures on various conditions. A sensitivity analysis of NbS parameterization, spanning a plausible range of parameter values, could also offer valuable insights into how parameter input uncertainties influence model output uncertainties. Improving the robustness in the estimation of NbS effectiveness remains a critical area for future research.

3.7.3 Multi-scale assessment of NbSs

Currently, there is no globally accepted or standardized approach for monitoring NbSs (Kumar et al., 2021a). The fine-scale spatial variability of NbS effectiveness demonstrated in our study highlights the need to carefully consider the appropriate methods and scales for monitoring these interventions after implementation. These approaches should be determined based on the type, extent and expected outcomes of NbSs. Macro-scale indicators at the watershed level, such as reductions in peak flow, can be highly useful for assessing the overall, cumulative effects of a combination of NbSs implemented upstream. While these indicators are useful, they do not distinguish the individual impacts of specific interventions. Therefore, fine-scale indicators are also necessary to optimize the benefits of NbSs. Nevertheless, these assume that individual interventions act independently and do not interact at small scales, an assumption that may not be always true.

Monitoring approaches often rely on ground-based measurement stations, which are often limited in number and scattered sparsely. Encouraging, recent advances in airborne and satellite-based remote sensing technologies may enable the possibility to bridge this gap between scales, improving the capacity to systematically monitor NbS performance (Kumar et al., 2021a). For instance, monitoring the effectiveness of NbSs in mitigating agricultural drought can be achieved through the comparison of vegetation indices such as the normalized difference vegetation index (NDVI) or the vegetation health index (Patel et al., 2012). Additionally, flood and estimated damage costs might be evaluated by GIS-based monitoring (Haq et al., 2012).

It is important to recognize that the proposed “Spatial variation of infiltration rates” indicator represents just one aspect of NbS effectiveness in reducing flood risks. While it provides valuable insights into how an NbS mitigates the hazard itself, a comprehensive assessment should also account for downstream exposure or vulnerability, which are crucial factors to consider when planning flood risk management (Klijn et al., 2015). For instance, interventions may reduce flood risks locally, but their impacts can also be felt much farther downstream (Lane, 2017). For a comprehensive assessment of the effectiveness of NbSs, it is crucial to adopt cross-scale approaches because the effects of NbSs are not necessarily confined to the immediate area of implementation. Therefore, in-depth evaluation of NbSs requires models capable of capturing complex interactions across spatial scales.

Furthermore, the effectiveness of NbSs may vary over time, depending on (i) meteorological variability (e.g. antecedent moisture conditions) and (ii) the time required for a measure to become fully effective. For instance, some measures, such as forest diversification, may take several years or even decades to reach full maturity. In the current study, while the effect of antecedent moisture conditions was effectively modelled, the transition phase between the BASELINE and POST scenarios was omitted. However, time-variable effects of NbSs may be important to consider when defining priorities in a specific local catchment context (Fennell et al., 2023b). These time-variable effects also raise the question of incorporating future climate scenarios into the modelling of NbS effectiveness (Gómez Martín et al., 2021).

3.7.4 Multi-objective assessment of NbSs

In the present article, the potential synergies or trade-offs that can arise with NbSs for mitigating floods and agricultural droughts were evaluated. Beyond these findings, other co-benefits or trade-offs may emerge, which should be assessed according to the specific context and interests of the region in question. For example, afforestation is known to increase watershed evapotranspiration, reducing river flow and groundwater levels. While this trade-off may be limited in humid climates (Hou et al., 2023; Zhang et al., 2017), which are also less prone to hydrological droughts, such practices, like planting alien trees, could be problematic in arid or semi-arid regions that are more sensitive to long-term hydrological and socioeconomic droughts (Rebelo et al., 2022). Hence, the impacts of these actions extend beyond hydrological effects, encompassing important social, cultural, economic and environmental aspects (Sowińska-Świerkosz and García, 2021). NbSs are inherently integrative and multifunctional; focusing solely on a hydrological performance without considering the broader social, economic, political (governance) and environmental context can undermine the success and legitimacy of projects (Brauman et al., 2022). This underscores the value of multi-objective optimization procedures, which help identify optimal scenarios that balance competing objectives (Strauch et al., 2019; Castro, 2022; Yang et al., 2023).

3.7.5 Bridging the gap between simulated scenarios and real-world implementation

As demonstrated in this study and many others, NbSs offer significant potential for mitigating floods and droughts. To translate this potential into real-world implementations, several challenges must be carefully addressed. Raška et al. (2022) identified four key aspects that influence the successful implementation of NbSs: (i) uncertainties regarding their effectiveness, (ii) decisions on optimal locations, (iii) institutional frameworks, and (iv) resource availability. The first two aspects primarily arise from the inherent variability of NbS effectiveness across spatial and temporal scales, as well as their multi-objective nature – issues that this study, somehow, aims to elucidate. By providing a quantitative assessment of NbS effectiveness, this research contributes to a more informed decision-making process for stakeholders and policymakers, creating a foundation for discussions on implementation strategies. Strong stakeholder engagement plays a key role in ensuring successful NbS deployment. NbS deployment may generate conflicting interests between different groups, creating both winners and losers (Bark et al., 2021). An interdisciplinary approach that integrates hydrology with social sciences such as sociology and economics can help promote transparency and fairness, particularly in designing effective compensation mechanisms that encourage participation. Moreover, aligning NbS implementation with governance structures is essential. While NbSs are most effective when they are planned at the watershed scale, decision-making often occurs within administrative boundaries that may not fully reflect hydrological realities. Strengthening coordination between authorities and clearly defining responsibilities will enhance the feasibility of NbS adoption.

4 Conclusion

We presented a process-based modelling approach, making full use of the available data, to evaluate the effectiveness of NbSs implemented in productive forest and agricultural catchments for flood and agricultural drought mitigation. Rather than focusing solely on the riverbed, the measures considered in this study target catchment slopes, aiming to enhance vertical water fluxes such as infiltration, recharge and evapotranspiration. These solutions are designed for multifunctionality and are generally easier to implement than hard engineering measures. The modelling approach allowed for the representation of fine-scale interactions between soil and water processes across different spatial scales, which may play a crucial role in the functioning of measures. It also enabled the computation and development of indicators to assess NbS potential for flood and agricultural drought mitigation at multiple scales, focusing on catchment slopes and outlet. Our results emphasized the critical role of spatial variability in NbS effectiveness, influenced by soil characteristics. A key finding is that the effectiveness of NbSs in reducing flood risk is significantly influenced by the soil's natural drainage characteristics, which vary across a watershed. In this study, soils with moderate drainage characteristics appeared to show the highest potential for enhanced infiltration through NbSs that increase evapotranspiration, such as hedgerows. Additionally, NbSs that enhanced root depth were most effective in mitigating agricultural drought, especially in water-limited environments. This study shows that NbSs can influence the hydrological functioning of a watershed as a whole, simultaneously affecting both floods and droughts. This underscores the need to move away from siloed evaluations of NbSs and rather adopt multi-scale and multi-objective assessment approaches. The study acknowledges the challenges of current modelling and monitoring techniques, particularly in the need for measured data at various spatial scales and uncertainty estimates, showing the way for future research. This study encourages decision-makers to envision more hydrologically resilient landscapes – an increasingly urgent necessity in our latitudes.

Appendix A

Figure A1 Vegetation development trajectories (LAI, rooting depth and Kc) used as model input to calculate evapotranspiration.