the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Wflow_sbm v0.7.3, a spatially distributed hydrological model: from global data to local applications
Willem J. van Verseveld
Albrecht H. Weerts
Martijn Visser
Joost Buitink
Ruben O. Imhoff
Hélène Boisgontier
Laurène Bouaziz
Dirk Eilander
Mark Hegnauer
Corine ten Velden
Bobby Russell
The wflow_sbm hydrological model, recently released by Deltares, as part of the Wflow.jl (v0.7.3) modelling framework, is being used to better understand and potentially address multiple operational and water resource planning challenges from a catchment scale to national scale to continental and global scale. Wflow.jl is a free and opensource distributed hydrological modelling framework written in the Julia programming language. The development of wflow_sbm, the model structure, equations and functionalities are described in detail, including example applications of wflow_sbm. The wflow_sbm model aims to strike a balance between lowresolution, lowcomplexity and highresolution, highcomplexity hydrological models. Most wflow_sbm parameters are based on physical characteristics or processes, and at the same time wflow_sbm has a runtime performance well suited for largescale highresolution model applications. Wflow_sbm models can be set a priori for any catchment with the Python tool HydroMTWflow based on globally available datasets and through the use of pointscale (pedo)transfer functions and suitable upscaling rules and generally result in a satisfactory (0.4 ≥ Kling–Gupta efficiency (KGE) < 0.7) to good (KGE ≥ 0.7) performance for discharge a priori (without further tuning). Wflow_sbm includes relevant hydrological processes such as glacier and snow processes, evapotranspiration processes, unsaturated zone dynamics, (shallow) groundwater, and surface flow routing including lakes and reservoirs. Further planned developments include improvements on the computational efficiency and flexibility of the routing scheme, implementation of a water demand and allocation module for water resource modelling, the addition of a deep groundwater concept, and computational efficiency improvements through for example distributed computing and graphics processing unit (GPU) acceleration.
 Article
(10076 KB)  Fulltext XML
 BibTeX
 EndNote
Hydrological models have proven to be useful tools in better understanding multiple operational and water resource planning challenges including drought (e.g. Trambauer et al., 2015) and flood forecasting (e.g. Alfieri et al., 2013), the assessment of water availability (e.g. van Beek et al., 2011), and analysing the impact of food production on river systems (e.g. Jägermeyr et al., 2017). An advantage of spatially distributed (gridded) hydrological models, in contrast to spatially lumped models, is the ability to directly use the spatially varying information contained in spatial datasets for model setup, forcing and validation. Highresolution spatial datasets have become increasingly available, often at a global scale, and can be used to represent land cover, vegetation (e.g. leaf area index (LAI)) and soil properties in spatially distributed hydrological models. For example, SoilGrids provides gridded soil information at 250 m resolution globally (Hengl et al., 2017). With regard to forcing, the release of the fifthgeneration ECMWF atmospheric reanalysis of the global climate (ERA5) (Hersbach et al., 2020) dataset (1959–present), with a spatial resolution of ∼ 31 km × 31 km and a temporal resolution of 1 h, and ERA5Land with a spatial resolution of ∼ 9 km × 9 km are worth mentioning. Recently, it has been argued that the development of a hyperresolution global hydrological model at 1 km^{2} or finer is a “grand challenge for hydrology” and is needed to address the water problems facing society (Wood et al., 2011; Bierkens et al., 2014).
Notwithstanding the advantages of and need for (hyperresolution) spatially distributed hydrological models, parameterization of these models is not straightforward because of overparameterization and as a result overfitting (Jakeman and Hornberger, 1993; Beven, 1993, 2006). Furthermore, transferability of hydrological parameters across spatial and temporal scales is important for reducing calibration time and the application of hydrological models in ungauged or poorly gauged basins. However the impact of transferring model parameters across spatial and temporal resolutions on model performance is unequivocal, and high parameter transferability across spatial resolution may also be the result of inadequate representation of spatial variability in (largescale) hydrological models (Melsen et al., 2016). Finally, there is the scientific debate on the “best” approach to processbased hydrological modelling leading to appropriate physical realism, especially related to model structure and model solutions (Kirchner, 2006; Clark et al., 2016, 2017).
Concerning hydrological model structure and solutions, Hrachowitz and Clark (2017) classified hydrological models along two dimensions, process complexity and spatial resolution. Hydrological models with high process complexity and spatial resolution are for example ParFlow (Kollet and Maxwell, 2006), HydroGeoSphere (Brunner and Simmons, 2012) and HYDRUS (2D/3D) (Ŝimůnek et al., 2008), while for example HBV (Bergström, 1992), SUPERFLEX (Fenicia et al., 2011) and FLEXTopo (Gao et al., 2014) are characterized by low spatial resolution and low process complexity. For the highresolution, highcomplexity hydrological models, the majority of the parameters are based on physical characteristics and may be estimated directly or by upscaling from field or remotely sensed observations, depending on the model resolution. For lowresolution, lowcomplexity hydrological models, the majority of parameters are effective parameters at the catchment scale and calibration is required to identify parameter values. Generally, highresolution, highcomplexity hydrological models are computationally demanding, which limits their application to smaller domains or requires either a reduction in model resolution or highperformance computing resources for largescale applications. Free and opensource spatially distributed hydrological models that require a low calibration effort (parameters based on physical characteristics) and have fast runtimes applicable to largescale highresolution modelling (medium complexity) are not available or have very limited availability to our knowledge.
Wflow_sbm represents a family of spatially distributed hydrological models that have the vertical hydrological simple bucket model (SBM; Vertessy and Elsenbeer, 1999) concept in common but can have different lateral concepts that control how water is routed for example over the land or river domain. This paper presents the wflow_sbm model configuration that makes use of the kinematicwave approach for river, overland and lateral subsurface flow. It is part of the opensource distributed hydrological model platform Wflow.jl (van Verseveld et al., 2024) developed at Deltares. Wflow_sbm strikes a balance between lowresolution, lowcomplexity and highresolution, highcomplexity hydrological models, giving an answer to most of the aforementioned challenges. In this model, the soil part is largely based on Topog_SBM (Vertessy and Elsenbeer, 1999), with gravitybased infiltration and vertical flow through the soil column as well as capillary rise representing a simplified version of Richards' equation. Furthermore it uses a 1D kinematicwave approach for channel, overland and lateral subsurface flows that is similar to TOPKAPI (Benning, 1995; Todini and Ciarapica, 2002), G2G (Bell et al., 2007), 1KDHM (Tanaka and Tachikawa, 2015) and Topog_SBM (Vertessy and Elsenbeer, 1999), as an approximation for dynamic waves and variably saturated subsurface flow (Richards' equation). Its advantage is that most wflow_sbm parameters are based on physical characteristics and at the same time wflow_sbm has a runtime performance well suited for largescale highresolution modelling.
Furthermore, in line with the need to improve the transparency, reproducibility and ease of setting up hydrological models (Clark et al., 2017; Knoben et al., 2021), we use the wflow plugin (HydroMTWflow; Eilander et al., 2022) of the HydroMT Python package (Eilander and Boisgontier, 2022) to set up wflow_sbm models for any catchment based on globally available datasets, e.g. SoilGrids (Hengl et al., 2017), GlobCover 2009 (Arino et al., 2010) and MERIT Hydro (Yamazaki et al., 2019). Pointscale (pedo)transfer functions (PTFs) from the literature are used to derive model parameters at the highest available resolution of the data and are scaled with suitable upscaling operators (Imhoff et al., 2020) to the desired model resolution. The advantage of this is that transfer functions are only constrained by field and laboratory measurements, although we acknowledge that the scale at which these PTFs can be applied remains uncertain (van Looy et al., 2017; Samaniego et al., 2017). Nevertheless, the application of this method to the Rhine Basin resulted, for most discharge gauging stations in the central and northern part of the basin, in Kling–Gupta efficiency (KGE; Gupta et al., 2009) values between 0.6 and 0.9 (Imhoff et al., 2020). In the meantime, wflow_sbm and the aforementioned approach were used and tested to model the basins in the upper region of the greater Chao Phraya River in Thailand (Wannasin et al., 2021a, b) and the Citarum River in Indonesia (Rusli et al., 2021). Meijer et al. (2021) used wflow_sbm to rapidly develop a water resource model for the upper Niger Basin using global online data. Sperna Weiland et al. (2021) used wflow_sbm to assess climate change impacts in nine river basins across Europe, while Aerts et al. (2022) used wflow_sbm to assess the impact of various model resolutions (200 m, 1 km, 3 km) on wflow_sbm performance for the CAMELSUS dataset.
The objective of this paper is to describe the wflow_sbm model in detail (model structure and equations) and to present some applications and envisaged future developments. Section 2 describes the development of the wflow_sbm model within the Wflow.jl framework and its model structure, model equations and functionalities. In Sect. 3 we describe the computational performance of wflow_sbm. Several applications of wflow_sbm are demonstrated in Sect. 4, followed by conclusions and foreseen future work in Sect. 5.
2.1 Overview
Wflow.jl (v0.7.3) (van Verseveld et al., 2024) is an opensource modelling framework for distributed hydrological modelling, containing multiple distributed hydrological model concepts implemented in the programming language Julia (Bezanson et al., 2017). It is a continuation of the wflow framework (Schellekens et al., 2020), which is based on the PCRaster Python framework (Karssenberg et al., 2010). The switch to the programming language Julia was made because Julia offers high performance (speed of C), required for largescale highresolution hydrological model applications, and is an “easytouse” language. Julia also opens up opportunities to parallelize the code for further improved computational performance. Wflow.jl provides several different vertical and lateral concepts that can be used for hydrological modelling and is compliant with the Basic Model Interface (BMI). Three vertical hydrological concepts are available within Wflow.jl: HBV96 (wflow_hbv), FLEXTopo (wflow_flextopo) and SBM (wflow_sbm).
Wflow_sbm is the main hydrological model concept of the Wflow.jl framework and represents a family of hydrological models that have the vertical SBM concept in common. Wflow_sbm can have different lateral concepts that control how water (river, overland and subsurface flow) is routed, easily enabled by the modular structure of Wflow.jl. The wflow_sbm model presented here (Fig. 1) consists of the vertical SBM concept, and for the lateral components, the kinematicwave approach is used for river, overland and lateral subsurface flow, similarly to TOPKAPI (Benning, 1995; Todini and Ciarapica, 2002), G2G (Bell et al., 2007), 1KDHM (Tanaka and Tachikawa, 2015) and Topog_SBM (Vertessy and Elsenbeer, 1999). The vertical SBM concept is largely based on Topog_SBM (Vertessy and Elsenbeer, 1999), which considers the soil to be a “bucket” with a saturated and unsaturated store. While Topog_SBM is specifically designed to simulate fastrunoff processes during discrete storm events in small catchments (<10 km^{2}) as evapotranspiration losses are ignored, wflow_sbm can be applied to a wider variety of catchments. The main differences between wflow_sbm and Topog_SBM are as follows:

the addition of evapotranspiration and interception losses;

the addition of a root water uptake reduction function (Feddes et al., 1978);

the addition of capillary rise;

the addition of glacier, snow buildup and melting processes and an avalanche option for downhill snow transport;

water being routed downstream over an eightdirection (D8) network, instead of the element network being based on contour lines and trajectories, used by Topog_SBM;

the introduction of an option to divide the soil column into different layers to allow for transfer of water within the unsaturated zone.
Wflow_sbm has been applied in various catchments around the world showing satisfactory (0.4 ≥ KGE < 0.7) to good (KGE ≥ 0.7) performance (e.g. López López et al., 2016; Hassaballah et al., 2017; Giardino et al., 2019; Gebremicael et al., 2019; Imhoff et al., 2020; LaverdeBarajas et al., 2020; Wannasin et al., 2021a, b; Rusli et al., 2021; Meijer et al., 2021).
Figure 1 presents the different processes and fluxes in the wflow_sbm model. Precipitation enters each grid cell through the interception routine (total precipitation is first intercepted), based on the Gash model (Gash, 1979) or a modified Rutter model (Rutter et al., 1971, 1975) depending on the time stamp the model is using. Throughfall and stemflow from the interception routine are transferred to the optional snow (based on the HBV96 hydrological model concept; Bergström, 1992) and glacier routines (based on the HBVlight degreedaybased model; Seibert et al., 2018). The soil in every grid cell is considered a single bucket, divided into a saturated and unsaturated store, with the option to divide the soil column into different layers. Available infiltration (stemflow and throughfall not converted into snow, including meltwater) infiltrates into the soil or becomes direct runoff based on the river fraction or openwater (excluding rivers) fraction. Soil infiltration is determined separately for the paved and nonpaved areas, as these have different infiltration capacities. Naturally, only the water that can be stored in the soil can infiltrate. If not all water can infiltrate, this is added as saturation excess water to the runoff routing scheme for overland flow. Infiltration excess occurs when the infiltration capacity is smaller than the available infiltration rate, and this amount of water is also added to the runoff routing scheme for overland flow. An exponential decay of the saturated hydraulic conductivity with soil depth is assumed. Transfer of water to the unsaturated store and to the saturated store is based on Brooks and Corey (1964) when the soil column is divided into different layers, and in the case of one soil layer the original Topog_SBM vertical transfer formulation can also be used. Part of the water evaporates through soil evaporation, transpiration that is first derived from the saturated store if roots intersect with the saturated store and then from the unsaturated store, and openwater (excluding rivers) and river evaporation. Besides transpiration, capillary rise and leakage result in a flux from the saturated store to the unsaturated store and outside of the model domain, respectively. The kinematicwave approach is used to route subsurface flow laterally. Saturation excess water occurs when the water table of lateral subsurface flow reaches the surface, and exfiltration of water from the unsaturated store to the surface because of a changing water table is added as saturation excess water to the runoff routing scheme for overland flow. For overland and river routing, the kinematicwave approach is also used. Reservoir and lake models (optional) can be included within the kinematicwave river routing.
The wflow_sbm model is described in more detail, including equations, in Sect. 2.2–2.7. These sections link to the main routines of wflow_sbm (Fig. 1):

“Interception” (Sect. 2.2)

“Snow and glaciers” (Sect. 2.3)

“The soil module and evapotranspiration” (Sect. 2.4)

“Lateral subsurface flow” (Sect. 2.5)

“Surface flow routing” (Sect. 2.6)

“Reservoirs and lakes” (Sect. 2.7).
Table A1 lists wflow_sbm state and flux variables (nonexhaustive). Additionally, wflow_sbm model inputs and parameters are listed in Table A2, including default values. Tables A1 and A2 both list the symbols that are used in Sect. 2.2–2.7 as well as the corresponding Wflow.jl names. It is possible to provide wflow_sbm with cyclic model parameters, for example monthly LAI maps (part of the HydroMTWflow default model setup), to estimate timedependent interception parameters (see also Sect. 2.2.3) or to provide model parameters as part of forcing (besides precipitation, potentialevapotranspiration and temperature fields). For symbols that represent model parameters in Sect. 2.2–2.7 (except in Sect. 2.2.3), a timedependent notation is not used because providing timedependent model parameters is optional.
The equations of most hydrological processes of the wflow_sbm model are solved using the explicit Euler method for each model time step t of length Δt [s] (model time step size) specified in the configuration file. For the unsaturated flow through multiple soil layers, variable internal time stepping is used, as a maximum change in soil water per unsaturated soil layer is allowed to prevent the overestimation of vertical unsaturated flows (Sect. 2.4.2). The kinematicwave equations (Sect. 2.5 and 2.6) are solved using Newton's method. For the river and overland flow kinematic wave, variable internal time stepping based on the Courant number or a fixed time step size [s] specified in the configuration file is used when iterations of the kinematic wave for surface flow are enabled in the configuration file; otherwise these equations are solved at model time step t of length Δt [s]. Reservoir and lake equations (Sect. 2.7) are solved using the explicit Euler method, except for the lake model with a parabolic weir, where the modified Puls approach is used. The reservoirs and lakes are part of the kinematicwave river routing, and these equations are solved at the internal time step of the kinematicwave river routing. The lateral subsurface flow is solved for each model time step using Newton's method, and depending on the model time step size, resolution and model parameters related to lateral subsurface flow, this may result in loss of accuracy. With a daily model time step, we estimate that using a minimum grid resolution of 200 m gives generally accurate lateralsubsurfaceflow results. The use of the explicit Euler method for most equations means that results may differ between simulations that use a different model time step (for example daily vs. hourly). Most flux variables in wflow_sbm are defined per model time step t, and external flux parameters (expressed per day (the model base time step size)) are converted to the userdefined model time step size, during model initialization. For model equations in this section handling these flux variables, the model time step t is implicitly embedded (e.g. for subtracting a flux variable from a storage variable, the flux variable is not multiplied by t in the equation), which is identical to the code implementation of these equations.
To run a wflow_sbm model, several files are required: (1) a configuration file in the Tom's Obvious Minimal Language (TOML) format; (2) a netCDF file containing static and (optional) cyclic data, for example model parameters, flow direction, river network and gauges; and (3) a netCDF file containing forcing data (precipitation, potentialevapotranspiration and temperature fields). Storage and rating curves for lakes should be provided in CSV format. The static and forcing maps should have the same spatial domain and resolution; i.e. the regridding of forcing data is not supported. The focus of Wflow.jl is on the computations (computational engine), and the modular structure of the code simplifies extending the base code for pre and postprocessing purposes.
In the TOML file the following aspects are defined: simulation period and model time step size, modelspecific settings like the model type (e.g. “sbm” for wflow_sbm) or whether to include snow modelling, locations and names of input and output files, mapping of internal model variables and parameters to external netCDF variables, optional modification of input model parameters and forcing, and output options. Glacier and snow modelling and lake and reservoir modelling are optional (specific model settings). Wflow_sbm typically runs at a daily time step (recommended maximum model time step size) and a spatial resolution of ∼ 1 km (we recommend a maximum grid resolution of ∼ 5 km). Subdaily model time steps are supported, for example for flow forecasting purposes or small (fastresponding) catchments. Output options consist of gridded data (netCDF) and scalar data (netCDF or CSV). Scalar data can be generated for individual grid cells or areas (e.g. subcatchment).
For users that mainly want to run simulations without installing Julia, Wflow.jl is available as a compiled executable (crossplatform). Users that want to explore and modify the code and want to extend Wflow.jl (e.g. writing your own Julia scripts around the Wflow.jl package), we recommend installing Wflow.jl as a Julia package. The Wflow.jl documentation provides more details about the Wflow.jl installation and usage in the “User guide” section (van Verseveld et al., 2024).
2.2 Interception
For interception the Gash model (Gash, 1979) for daily (or larger) time steps and a modified Rutter model (Rutter et al., 1971, 1975) for subdaily time steps are available within the Wflow.jl framework. The Gash model is a stormbased interception model, assuming one precipitation event per model time step and is applied when wflow_sbm runs at a daily (or larger) time step. For subdaily time steps, a modified Rutter model is used.
2.2.1 Gash model
The original Gash model considers precipitation input to be a series of discrete storm events, where each storm event is divided into three sequential phases: (1) wetting phase during which precipitation saturates the canopy, (2) saturation phase during which the canopy is saturated and the precipitation intensity is higher than the average evaporation rate of the saturated canopy, and (3) drying phase during the period after precipitation has ceased. The precipitation rate to completely saturate the canopy ${P}_{\mathrm{sat},\phantom{\rule{0.125em}{0ex}}\mathrm{max}}^{t}$ [mm t^{−1}] at time step t is defined as
where ${P}_{\mathrm{sat}}^{t}$ [mm t^{−1}] is the average precipitation intensity and ${E}_{\mathrm{sat}}^{t}$ [mm t^{−1}] is the average evaporation rate during saturation of the canopy, S_{canopy, max} [mm] is the canopy storage capacity, f_{canopygap} [–] is the canopy gap fraction, and f_{stemflow} [–] is the stemflow fraction. When wflow_sbm is not provided with the leaf area index (LAI), the ratio $\frac{{E}_{\mathrm{sat}}^{t}}{{P}_{\mathrm{sat}}^{t}}$ [–] is used as the model parameter (Table A2); otherwise ${E}_{\mathrm{sat}}^{t}=(\mathrm{1}{f}_{\mathrm{canopygap}}){E}_{\mathrm{pot},\phantom{\rule{0.125em}{0ex}}\mathrm{total}}^{t}$, where ${E}_{\mathrm{pot},\phantom{\rule{0.125em}{0ex}}\mathrm{total}}^{t}$ [mm t^{−1}] is potential evapotranspiration (for the computation of f_{canopygap} as a function of LAI, see Sect. 2.2.3), ${P}_{\mathrm{sat}}^{t}$ is then equal to the total precipitation input rate P^{t} [mm t^{−1}], and f_{canopygap} is time dependent when LAI is provided as a cyclic (e.g. monthly) parameter or as part of forcing. The stemflow fraction f_{stemflow} in wflow_sbm is defined as a fixed fraction (0.1) of f_{canopygap} limited by the canopy fraction (1−f_{canopygap}), and stemflow ${P}_{\mathrm{stemflow}}^{t}$ [mm t^{−1}] at time step t is calculated as follows:
Interception during the wetting phase ${I}_{\mathrm{wet}}^{t}$ [mm t^{−1}], saturation phase ${I}_{\mathrm{sat}}^{t}$ [mm t^{−1}] and dry phase ${I}_{\mathrm{dry}}^{t}$ [mm t^{−1}] at time step t is given by Eqs. (4)–(6):
The total interception ${I}_{\mathrm{total}}^{t}$ [mm t^{−1}] at time step t, assuming that trunk interception can be neglected, is the sum of the interception in all three phases, bounded by ${E}_{\mathrm{pot},\phantom{\rule{0.125em}{0ex}}\mathrm{total}}^{t}$:
Throughfall ${P}_{\mathrm{throughfall}}^{t}$ [mm t^{−1}] at time step t is the remainder after subtracting the total interception and stemflow from the precipitation:
The remaining potential evaporation ${E}_{\mathrm{pot},\phantom{\rule{0.125em}{0ex}}\mathrm{remainder}}^{t}$ [mm t^{−1}] at time step t is given by
2.2.2 Modified Rutter model
For subdaily time steps, a modified Rutter interception model is used, which compared to the Gash model keeps track of the canopy storage ${S}_{\mathrm{canopy}}^{t}$ [mm] and is updated in two steps. Stemflow is calculated in the same way as the Gash model; see Eq. (3). The precipitation rate on the canopy ${P}_{\mathrm{canopy}}^{t}$ [mm t^{−1}] at time step t is a function of the total precipitation input rate and the canopy gap and stemflow fractions:
The initial drainage rate ${D}_{\mathrm{canopy},\mathrm{s}\mathrm{1}}^{t}$ [mm t^{−1}] from the canopy storage at time step t is the surplus of canopy storage at the previous time step compared to the canopy storage capacity S_{canopy, max} [mm]:
This check is required because S_{canopy, max} (and f_{canopygap}) can change over time (see Sect. 2.2.3). The canopy storage is then updated based on the initial canopy drainage rate and precipitation that falls on the canopy (Eq. 12); then the evaporation rate from the canopy storage ${E}_{\mathrm{canopy}}^{t}$ [mm t^{−1}] is determined (Eq. 13) and subtracted from ${S}_{\mathrm{canopy}}^{t}$ (Eq. 14):
The remaining potential evaporation ${E}_{\mathrm{pot},\phantom{\rule{0.125em}{0ex}}\mathrm{remainder}}^{t}$ [mm t^{−1}] at time step t is given by
The canopy storage ${S}_{\mathrm{canopy}}^{t}$ is drained again if required with drainage rate ${D}_{\mathrm{canopy},\mathrm{s}\mathrm{2}}^{t}$ at time step t,
and ${S}_{\mathrm{canopy}}^{t}$ is updated to get the final canopy storage ${S}_{\mathrm{canopy}}^{t}$ [mm]:
Throughfall ${P}_{\mathrm{throughfall}}^{t}$ [mm t^{−1}] at time step t is given by adding the total drainage rate from the canopy and the precipitation rate directly on the ground:
The total interception ${I}_{\mathrm{total}}^{t}$ [mm t^{−1}] at time step t is given by
2.2.3 Interception model parameters from leaf area index (LAI)
Within wflow_sbm it is possible to estimate interception model parameters based on LAI [m^{2} m^{−2}] as static input (generally not used), as timedependent input such as cyclic time series (climatology) or as part of forcing. It is assumed that the canopy capacity for leaves ${S}_{\mathrm{leaf},\phantom{\rule{0.125em}{0ex}}\mathrm{max}}^{t}$ is linearly related to LAI^{t} at time step t through the specific leaf storage S_{leaf} [mm] (van Dijk and Bruijnzeel, 2001):
The specific leaf storage is related to land cover type. Also the storage for the woody part of the vegetation S_{wood, max} [mm] is required to estimate total canopy capacity ${S}_{\mathrm{canopy},\phantom{\rule{0.125em}{0ex}}\mathrm{max}}^{t}$ [mm] at time step t. The relations between land cover and S_{leaf} and S_{wood, max} are based on Pitman (1989) and Liu (1998). The canopy gap fraction ${f}_{\mathrm{canopygap}}^{t}$ [–] at time step t is determined using the extinction coefficient k based on van Dijk and Bruijnzeel (2001), and the value of k is related to land cover type:
2.3 Snow and glaciers
2.3.1 Snow
Snow processes are adopted from the HBV96 hydrological model concept (Bergström, 1992). The effective precipitation rate ${P}_{\mathrm{effective}}^{t}$ [mm t^{−1}] (throughfall and stemflow) occurs as snowfall ${P}_{\mathrm{snow}}^{t}$ [mm t^{−1}] at time step t if the air temperature ${T}_{\mathrm{air}}^{t}$ [°C] at time step t is below a userdefined temperature threshold s_{fall, T threshold} [°C]. An interval parameter s_{fall, T interval} [°C] defines the range over which precipitation falls partly as snow and partly as rain, with 100 % snow at the lower end and decreasing linearly to 0 % at the upper end. The fraction of precipitation that occurs as rainfall ${f}_{\mathrm{rain}}^{t}$ [–] at time step t is calculated as
This fraction is used to calculate snowfall ${P}_{\mathrm{snow}}^{t}$ and rainfall ${P}_{\mathrm{rain}}^{t}$ [mm t^{−1}] at time step t as follows:
For snowmelt, HBV96 uses a degreeday approach, an empirical relationship between melt and air temperature. If ${T}_{\mathrm{air}}^{t}$ is above a melting temperature threshold s_{melt, T threshold} [°C], snowmelt occurs. The potential snowmelt rate ${M}_{\mathrm{snow},\mathrm{pot}}^{t}$ [mm t^{−1}] at time step t, using the degreeday factor s_{ddf} [mm t^{−1} °C^{−1}], is calculated as follows:
The actual snowmelt rate ${M}_{\mathrm{snow},\mathrm{act}}^{t}$ [mm t^{−1}] at time step t is limited by the snow storage ${S}_{\mathrm{snow}}^{t\mathrm{1}}$ [mm] at the previous time step and is calculated by taking the minimum of ${M}_{\mathrm{snow},\mathrm{pot}}^{t}$ and ${S}_{\mathrm{snow}}^{t\mathrm{1}}$. The snowpack retains water that can refreeze if ${T}_{\mathrm{air}}^{t}$ is below s_{melt, T threshold}. The potential refreezing rate ${M}_{\mathrm{refreeze},\mathrm{pot}}^{t}$ [mm t^{−1}] at time step t is controlled by s_{ddf}, a coefficient of refreezing s_{refreeze} [–] (fixed: 0.05), ${T}_{\mathrm{air}}^{t}$ and s_{melt, T threshold}, as follows:
The actual refreezing rate ${M}_{\mathrm{refreeze},\mathrm{act}}^{t}$ [mm t^{−1}] is based on the liquid water content of snow ${S}_{\mathrm{snow},\mathrm{liquid}}^{t\mathrm{1}}$ [mm] at the previous time step and the potential refreezing rate ${M}_{\mathrm{refreeze},\mathrm{pot}}^{t}$, by taking the minimum of ${M}_{\mathrm{refreeze},\mathrm{pot}}^{t}$ and ${S}_{\mathrm{snow},\mathrm{liquid}}^{t\mathrm{1}}$. Snow storage ${S}_{\mathrm{snow}}^{t}$ [mm] at time step t is then a function of snow storage ${S}_{\mathrm{snow}}^{t\mathrm{1}}$ at the previous time step, snowfall, the actual refreezing rate and actual snowmelt at time step t:
The liquid water content of snow ${S}_{\mathrm{snow},\mathrm{liquid}}^{t}$ at time step t is a function of the liquid water content of snow ${S}_{\mathrm{snow},\mathrm{liquid}}^{t\mathrm{1}}$ at the previous time step; the actual refreezing rate, actual snowmelt and rainfall rate at time step t; and the maximum amount of water that the snowpack can hold. This maximum amount is controlled by the waterholding capacity s_{whc} [–] of snow and snow storage at time step t:
The amount that does exceed the waterholding capacity of snow ($\mathrm{max}({S}_{\mathrm{snow},\mathrm{liquid}}^{t}{S}_{\mathrm{snow}}^{t}{s}_{\mathrm{whc}},\mathrm{0})$) is available as rainfall at time step t.
To control unlimited buildup of the snowpack at high altitude where temperature rarely reaches above the melting temperature, the optional avalanche routine can be used to transport snow, based on the local drain network, downhill, where it becomes available for snowmelt. The fraction of snow that can be transported downhill is calculated as
with ${f}_{\text{snow transport}}^{t}$ [–] the fraction of snow at time step t that is available for transport downhill, c_{land slope} [m m^{−1}] the slope of the land surface and s_{max} [mm] the maximum snowpack with a fixed value of 10 000 mm. The fraction of snow that can be transported downhill is multiplied with the snowpack storage and represents the transport capacity of snow ${M}_{\mathrm{snow},\mathrm{downhill}}^{t}$ [mm t^{−1}] at time step t:
Snow is then transported downhill based on the local drain network and the transport capacity of snow, which limits the snow transport, updating the snow storage ${S}_{\mathrm{snow}}^{t}$ [mm] and liquid water content of snow ${S}_{\mathrm{snow},\mathrm{liquid}}^{t}$ [mm] in each grid cell at time step t.
2.3.2 Glaciers
Glacier modelling considers two main processes: glacier buildup from snow turning into firn (adopted from the HBVlight model; Seibert et al., 2018) and glacier melt (using a temperature degreeday model). This glacier modelling approach considers firn to be a part of the accumulated glacier mass. First, a fixed fraction g_{snow to firn} [t^{−1}], which typically ranges between 0.001 and 0.006 for a daily time step, of snow storage ${S}_{\mathrm{snow}}^{t}$ [mm] on top of the glacier is converted into firn for each time step t of length Δt [s]:
where ${S}_{\text{snow to firn}}^{t}$ [mm t^{−1}] is the snowintofirn conversion rate at time step t, with a maximum conversion rate of 8 mm d^{−1}. This maximum conversion rate is scaled by Δt and the model base time step size Δt_{b} of 86 400 [s]. The snow storage from the snow module (Sect. 2.3.1) ${S}_{\mathrm{snow}}^{t}$ [mm] at time step t is then updated as follows:
with f_{glacier} [–] the fraction of a grid cell covered by a glacier. When the snowpack on top of the glacier is almost all melted (${S}_{\mathrm{snow}}^{t}<\mathrm{10}$ mm), glacier melt is enabled and estimated with a degreeday model. If the air temperature ${T}_{\mathrm{air}}^{t}$ is above a melting temperature threshold g_{melt, T threshold} [°C], glacier melt occurs. The potential glacier melt ${M}_{\mathrm{glacier},\mathrm{pot}}^{t}$ [mm t^{−1}], using the degreeday factor g_{ddf} [mm t^{−1} °C^{−1}], is calculated as
The actual glacier melt ${M}_{\mathrm{glacier},\mathrm{act}}^{t}$ [mm t^{−1}] at time step t is limited by the glacier storage at the previous time step ${S}_{\mathrm{glacier}}^{t\mathrm{1}}$ [mm] (expressed in mm water equivalent) and ${S}_{\text{snow to firn}}^{t}$ as follows:
The glacier storage ${S}_{\mathrm{glacier}}^{t}$ [mm] at time step t is then updated as follows:
A map with S_{glacier} values can be provided as an initial state (default: 5500 mm) when wflow_sbm is initialized with default values in the code (“cold” start); see also Table A2.
2.4 The soil module and evapotranspiration
2.4.1 Infiltration
The infiltration rate of available water ${F}_{\mathrm{available}}^{t}$ [mm t^{−1}] at time step t into the soil (throughfall, stemflow, snowmelt and glacier melt) is first added to the saturated parts of the grid cell: the river flow and overland flow components of wflow_sbm. This is based on the river fraction f_{river} [–] (river flow component) and openwater fraction (excluding rivers) f_{open water} [–] (overland flow component) within a grid cell as follows:
where ${R}_{\mathrm{river}}^{t}$ [mm t^{−1}] is runoff from the river fraction in a cell at time step t and ${R}_{\text{open water}}^{t}$ [mm t^{−1}] is runoff from the openwater fraction in a cell at time step t. ${R}_{\mathrm{river}}^{t}$ and ${R}_{\text{open water}}^{t}$ are later added to the wflow_sbm river and overland flow components, respectively. The infiltration rate of the remaining available water ${F}_{\mathrm{available}}^{t}$ [mm t^{−1}] at time step t into the soil is determined as
The soil in wflow_sbm is considered a bucket with a depth z_{soil} [mm] and is divided into a saturated store S_{sat} [mm] and an unsaturated store S_{unsat} [mm]. The top of the saturated store forms a pseudowater table at depth z_{watertable} [mm] such that the value of S_{sat} is given by
where θ_{s} and θ_{r} are the saturated and residual soil water contents, respectively, both expressed in mm mm^{−1}. The amount of water that can infiltrate is a function of the infiltration capacity c_{infiltration, paved} [mm t^{−1}] of the compacted soil (or paved area) fraction (f_{paved} [–]) of each grid cell, the infiltration capacity c_{infiltration, unpaved} [mm t^{−1}] of the noncompacted soil fraction (or unpaved area) ((1−f_{paved})) of each grid cell, the initial storage capacity of the unsaturated zone ${S}_{\mathrm{unsat},\mathrm{max}}^{t}$ [mm] (Eq. 43) at time step t and an optional reduction factor f_{frozen} applied to the infiltration capacity when snow is modelled, and the model setting soilinfreduction is set to true. The parameter f_{frozen} depends on the nearsurface soil temperature, which is modelled based on the approach of Wigmosta et al. (2009):
where ${T}_{\mathrm{soil}}^{t}$ [°C] is the nearsurface soil temperature at time step t, ${T}_{\mathrm{air}}^{t}$ [°C] is the air temperature at time step t, ${T}_{\mathrm{soil}}^{t\mathrm{1}}$ [°C] is the nearsurface soil temperature at the previous time step and w_{soil} is a weighting coefficient [t^{−1}]. The optional infiltration capacity reduction factor ${f}_{\mathrm{frozen}}^{t}$ at time step t is based on the model parameter f_{red, frozen} [–] and the nearsurface soil temperature as follows:
where $b=\frac{\mathrm{1}}{(\mathrm{1}{f}_{\mathrm{red},\mathrm{frozen}})}$, a=0 and c=8. The initial storage capacity of the unsaturated zone ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}\mathrm{max}}^{t}$ [mm] at time step t is based on the saturated storage ${S}_{\mathrm{sat}}^{t\mathrm{1}}$ [mm] at the previous time step, the sum of unsaturated storage in each soil layer n for m unsaturated soil layers ${\sum}_{n=\mathrm{1}}^{m}{S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}n}^{t\mathrm{1}}$ [mm] at the previous time step, and the total soil water capacity of the wflow_sbm soil bucket. ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}\mathrm{max}}^{t}$ is calculated as follows:
The infiltration rate of remaining available water is split into two parts, the part that falls on compacted areas ${F}_{\mathrm{available}}^{t}{f}_{\mathrm{paved}}$ [mm t^{−1}] and the part that falls on noncompacted areas ${F}_{\mathrm{available}}^{t}(\mathrm{1}{f}_{\mathrm{paved}})$ [mm t^{−1}] at time step t. The maximum infiltration rate in these areas is calculated by taking the minimum of the infiltration capacity and the infiltration rate for these areas:
The actual total infiltration rate ${F}_{\mathrm{total}}^{t}$ [mm t^{−1}] is a function of the total infiltration rate (compacted and noncompacted areas) and the initial unsaturated storage capacity:
Finally, the infiltration excess water rate ${F}_{\mathrm{excess}}^{t}$ [mm t^{−1}] at time step t is determined as
2.4.2 Soil water accounting scheme
The soil bucket in wflow_sbm with a depth z_{soil} [mm] can be split up into different layers. Assuming a unit head gradient, the potential transfer of water ${Q}_{\mathrm{transfer},\phantom{\rule{0.125em}{0ex}}\mathrm{pot},\phantom{\rule{0.125em}{0ex}}n}^{t}$ [mm t^{−1}] from an unsaturated layer n at time step t is controlled by the vertical saturated hydraulic conductivity ${K}_{\mathrm{vz},\phantom{\rule{0.125em}{0ex}}n}^{t}$ [mm t^{−1}] at time step t (Eq. 50) the effective saturation degree of layer n at time step t, and a Brooks–Corey power coefficient c_{n} [–] based on the pore size distribution index λ_{n} [–] (Brooks and Corey, 1964) of layer n:
where ${\mathit{\theta}}_{n}^{t}$ [mm mm^{−1}] is the soil water content of unsaturated soil layer n at time step t and θ_{s} and θ_{r} are as previously defined. The vertical saturated hydraulic conductivity ${K}_{\mathrm{vz},\phantom{\rule{0.125em}{0ex}}n}^{t}$ of unsaturated soil layer n for m unsaturated soil layers (based on the pseudowatertable depth at the previous time step ${z}_{\mathrm{watertable}}^{t\mathrm{1}}$ [mm]; n=1 refers to the upper soil layer) at time step t is given by
where the model parameter K_{v0} [mm t^{−1}] is the vertical saturated hydraulic conductivity at the soil surface (that declines exponentially with depth), f_{Kv, n} is an optional (default: 1.0) multiplication factor [–] for each soil layer n to correct the vertical saturated hydraulic conductivity, f_{Kv} is a scaling parameter [mm^{−1}] and z_{bottom, n} [mm] is the soil depth at the bottom of soil layer n. The thickness ${z}_{n,\phantom{\rule{0.125em}{0ex}}\mathrm{thickness}}^{t}$ [mm] of unsaturated soil layer n for m unsaturated soil layers at time step t is given by
where z_{n, thickness} is the actual thickness of soil layer n and ${z}_{\mathrm{bottom},\phantom{\rule{0.125em}{0ex}}n\mathrm{1}}$ [mm] is the soil depth at the bottom of the soil layer above soil layer n. For each unsaturated soil layer n, the transfer of water at time step t for m unsaturated soil layers is calculated as follows:
where ${Q}_{\mathrm{in},\phantom{\rule{0.125em}{0ex}}n}^{t}$ for unsaturated soil layer n=1 (upper soil layer) is the actual total infiltration rate ${F}_{\mathrm{total}}^{t}$ (Eq. 46) and, for unsaturated soil layer n>1, ${Q}_{\mathrm{in},\phantom{\rule{0.125em}{0ex}}n}^{t}$ is the actual transfer of water from the soil layer above layer n (${Q}_{\mathrm{transfer},\mathrm{act},n\mathrm{1}}^{t}$ [mm t^{−1}]); ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}n}^{t\mathrm{1}}$ [mm] is the unsaturated storage of layer n at the previous time step; and ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}n}^{t}$ [mm] is the subsequent updated unsaturated storage of layer n at time step t. During each model time step t, internal time stepping (iterations) is applied to Eqs. (53), (54) and (55), under unsaturated conditions based on a maximum allowed change in S_{unsat, n} of 0.2 mm for each internal time step to prevent an overestimation of ${Q}_{\mathrm{transfer},\phantom{\rule{0.125em}{0ex}}\mathrm{pot},\phantom{\rule{0.125em}{0ex}}n}^{t}$.
When the soil bucket in wflow_sbm is not split up into different layers, it is possible to use the original Topog_SBM vertical transfer formulation. The transfer of water from the unsaturated store to the saturated store is in that case controlled by the vertical saturated hydraulic conductivity at depth ${z}_{\mathrm{watertable}}^{t\mathrm{1}}$, an optional multiplication factor f_{Kv, 1} [–] to correct the vertical saturated hydraulic conductivity, the ratio between the unsaturated storage ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}\mathrm{1}}^{t}$ [mm] (resulting from Eq. 52) and the saturation deficit ${S}_{\mathrm{deficit}}^{t}$ [mm], and the available unsaturated storage ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}\mathrm{1}}^{t}$ at time step t:
2.4.3 Evapotranspiration
Openwater evaporation from waterbodies (excluding rivers) ${E}_{\text{open water}}^{t}$ [mm t^{−1}] and rivers ${E}_{\mathrm{river}}^{t}$ [mm t^{−1}] at time step t is based on the fraction of open water f_{open water} [–], the fraction of rivers f_{river} [–], the water level in the kinematic reservoir of the river flow component ${S}_{\mathrm{wl},\phantom{\rule{0.125em}{0ex}}\mathrm{river}}^{t\mathrm{1}}$ [mm] and the overland flow component ${S}_{\mathrm{wl},\phantom{\rule{0.125em}{0ex}}\mathrm{land}}^{t\mathrm{1}}$ [mm] at the previous time step, and the remaining potential evaporation after interception ${E}_{\mathrm{pot},\phantom{\rule{0.125em}{0ex}}\mathrm{remainder}}^{t}$ [mm t^{−1}] as follows:
The potentialevaporation rate remaining after interception (Eq. 9 or 15) and openwater evaporation (rivers and waterbodies (excluding rivers)) ${E}_{\mathrm{pot},\phantom{\rule{0.125em}{0ex}}\mathrm{remainder}}^{t}$ [mm t^{−1}] at time step t is then
Potential soil evaporation ${E}_{\mathrm{pot},\phantom{\rule{0.125em}{0ex}}\mathrm{soil}}^{t}$ [mm t^{−1}] at time step t is based on ${E}_{\mathrm{pot},\phantom{\rule{0.125em}{0ex}}\mathrm{remainder}}^{t}$ and the canopy gap fraction f_{canopygap} [–] (assumed to be identical to the amount of bare soil and can vary in time; Sect. 2.2.3). When the soil bucket in wflow_sbm is not split up into different layers, soil evaporation ${E}_{\mathrm{act},\phantom{\rule{0.125em}{0ex}}\mathrm{soil}}^{t}$ [mm t^{−1}] is calculated as follows:
with ${S}_{\mathrm{deficit}}^{t}$, z_{soil}, θ_{s}, θ_{r} and ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}\mathrm{1}}^{t}$ (Eq. 55) as previously defined. As such, soil evaporation will be potential if the soil is fully wetted, and it decreases linearly with increasing soil moisture deficit, limited by ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}\mathrm{1}}^{t}$. When the soil bucket in wflow_sbm is split up into different layers, soil evaporation ${E}_{\mathrm{act},\phantom{\rule{0.125em}{0ex}}\mathrm{soil}}^{t}$ [mm t^{−1}] is restricted to the upper layer. As with the case of a single soil layer, potential soil evaporation is scaled according to the wetness of the soil layer, based on the unsaturated layer storage from Eq. (55), as follows:
where z_{1, thickness} [mm] is the actual thickness of the upper layer and θ_{s}, θ_{r} and ${z}_{\mathrm{watertable}}^{t\mathrm{1}}$ are as previously defined. Soil evaporation ${E}_{\mathrm{act},\phantom{\rule{0.125em}{0ex}}\mathrm{soil}}^{t}$ is subtracted from the unsaturated storage in the upper soil layer ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}\mathrm{1}}^{t}$ (Eq. 55), and the remaining potential soil evaporation ${E}_{\mathrm{remainder},\phantom{\rule{0.125em}{0ex}}\mathrm{soil}}^{t}$ is determined as follows:
When the soil bucket in wflow_sbm is split up into different layers, soil evaporation ${E}_{\mathrm{act},\phantom{\rule{0.125em}{0ex}}\mathrm{soil},\phantom{\rule{0.125em}{0ex}}\mathrm{sat}}^{t}$ [mm t^{−1}] from the saturated store is possible when the water table ${z}_{\mathrm{watertable}}^{t\mathrm{1}}$ is present in the upper soil layer with actual thickness z_{1, thickness} [mm], and it is calculated as follows:
with θ_{s} and θ_{r} as previously defined. In the case of a single soil layer or when ${z}_{\mathrm{watertable}}^{t\mathrm{1}}$ is not present in the upper soil layer, ${E}_{\mathrm{act},\phantom{\rule{0.125em}{0ex}}\mathrm{soil},\phantom{\rule{0.125em}{0ex}}\mathrm{sat}}^{t}$ is set at zero. ${E}_{\mathrm{act},\phantom{\rule{0.125em}{0ex}}\mathrm{soil},\phantom{\rule{0.125em}{0ex}}\mathrm{sat}}^{t}$ is subtracted from the saturated store (at the previous time step):
Potential transpiration ${E}_{\mathrm{pot}\phantom{\rule{0.125em}{0ex}}\mathrm{trans}}^{t}$ [mm t^{−1}] at time step t is based on the remaining potential evaporation after interception and openwater evaporation (Eq. 61) and the canopy gap fraction as follows:
In wflow_sbm, transpiration is first taken from the saturated store if the roots reach the water table ${z}_{\mathrm{watertable}}^{t\mathrm{1}}$ at the previous time step. The fraction of wet roots f_{wet roots} [–] (ranging between 0 and 1) is determined using a sigmoid function that defines the sharpness of the transition between fully wet and fully dry roots. Transpiration ${E}_{\mathrm{trans},\phantom{\rule{0.125em}{0ex}}\mathrm{sat}}^{t}$ [mm t^{−1}] from the saturated store at time step t is calculated as follows:
where c_{rd} [mm^{−1}] is a model parameter that controls the sharpness of the sigmoid function and z_{rooting} [mm] is the rooting depth. The saturated store is then updated as follows:
The remaining potential transpiration ${E}_{\mathrm{pot}\phantom{\rule{0.125em}{0ex}}\mathrm{trans},\phantom{\rule{0.125em}{0ex}}\mathrm{remainder}}^{t}$ [mm t^{−1}] available for transpiration from the unsaturated store is calculated as follows:
The maximum allowed water extraction rate by roots ${E}_{\mathrm{root},\phantom{\rule{0.125em}{0ex}}\mathrm{max},\phantom{\rule{0.125em}{0ex}}n}^{t}$ [mm t^{−1}] from unsaturated soil layer n at time step t is a function of the fraction of roots ${f}_{\mathrm{roots},\phantom{\rule{0.125em}{0ex}}n}^{t}$ [–] of unsaturated layer n and the available unsaturated storage [mm] of layer n:
Next, a root water uptake reduction model based on Feddes et al. (1978) is used to calculate a reduction coefficient as a function of soil matric suction. Soil matric suction is calculated following Brooks and Corey (1964):
where h_{n} is the soil matric suction [cm] of unsaturated soil layer n; h_{b} is the air entry value [cm]; and θ_{n}, θ_{r}, θ_{s} and λ_{n} are as previously defined. In wflow_sbm, soil matric suction ${h}_{n}^{t}$ for each unsaturated soil layer n with thickness ${z}_{n,\phantom{\rule{0.125em}{0ex}}\mathrm{thickness}}^{t}$ [mm] at time step t is calculated as follows:
where ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}\mathrm{1}}^{t}$ is given by Eq. (65) and, for unsaturated layers n>1, ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}n}^{t}$ is given by Eq. (55). The root water uptake reduction coefficient ${f}_{\text{root uptake,}\phantom{\rule{0.125em}{0ex}}n}^{t}$ at time step t with ${h}_{n}^{t}$ below or equal to h_{3} (400 cm) is set to 1, with ${h}_{n}^{t}$ above or equal to h_{4} (15 849 cm) is set to 0, and with ${h}_{n}^{t}$ between h_{3} and h_{4} declines linearly from 1 to 0. The values for h_{2} (100 cm), h_{3} and h_{4} are fixed, and h_{1} is set by h_{b} (default: 10 cm); h_{b} can be defined as input to the model. In the original transpiration reduction curve of Feddes et al. (1978), root water uptake below h_{1} is set to zero (oxygen deficit) and between h_{1} and h_{2} root water uptake is limited. The assumption that very wet conditions do not affect root water uptake too much is probably generally applicable to natural vegetation; however for crops this assumption is not valid. This could be improved in the Wflow.jl code by applying the reduction to crops only. While the h_{3} value is fixed, in the original transpiration reduction curve of Feddes et al. (1978), h_{3} varies with the potential transpiration rate, and this could also be improved in the code. For unsaturated soil layer n, transpiration ${E}_{\mathrm{trans},\phantom{\rule{0.125em}{0ex}}\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}n}^{t}$ [mm t^{−1}] is controlled by ${E}_{\mathrm{root},\phantom{\rule{0.125em}{0ex}}\mathrm{max},\phantom{\rule{0.125em}{0ex}}n}^{t}$, the remaining potential transpiration ${E}_{\mathrm{pot}\phantom{\rule{0.125em}{0ex}}\mathrm{trans},\phantom{\rule{0.125em}{0ex}}\mathrm{remainder}}^{t}$ [mm t^{−1}] (for soil layer n=1 see Eq. 73, and for layers n>1 see Eq. 79), the unsaturated storage ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}n}^{t}$ [mm] (for soil layer n=1 see Eq. 65, and for layers n>1 see Eq. 55) and ${f}_{\text{root uptake,}\phantom{\rule{0.125em}{0ex}}n}^{t}$ at time step t:
At the same time ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}n}^{t}$ and the remaining potential transpiration ${E}_{\mathrm{pot}\phantom{\rule{0.125em}{0ex}}\mathrm{trans},\phantom{\rule{0.125em}{0ex}}\mathrm{remainder}}^{t}$ [mm t^{−1}] are updated by subtracting ${E}_{\mathrm{trans},\phantom{\rule{0.125em}{0ex}}\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}n}^{t}$:
After the soil water transfer, evaporation and transpiration computations, a soil water balance check is performed. Unsaturated storage that exceeds the maximum storage per layer is transferred to the layer above (or surface), from the bottom to the top unsaturated soil layer, resulting in an excess water rate at the surface ${R}_{\mathrm{excess},\phantom{\rule{0.125em}{0ex}}\mathrm{unsat}}^{t}$ [mm t^{−1}]. The actual infiltration rate ${F}_{\mathrm{act}}^{t}$ [mm t^{−1}] is then calculated as follows:
with ${F}_{\mathrm{total}}^{t}$ as previously defined. The rate of water that cannot infiltrate due to saturated soil conditions ${F}_{\mathrm{excess},\phantom{\rule{0.125em}{0ex}}\mathrm{sat}}^{t}$ [mm t^{−1}] is determined as
with ${F}_{\mathrm{available}}^{t}$, ${F}_{\mathrm{total}}^{t}$ and ${F}_{\mathrm{excess}}^{t}$ as previously defined.
Capillary rise in wflow_sbm is determined when an unsaturated zone occurs in the soil column at time step t based on the water table depth at the previous time step (${z}_{\mathrm{watertable}}^{t\mathrm{1}}>\mathrm{0}$), using the following approach. First a maximum capillary rise ${C}_{\mathrm{max}}^{t}$ [mm t^{−1}] at time step t is determined from the minimum of ${K}_{\mathrm{vz},\phantom{\rule{0.125em}{0ex}}m}^{t}$ (Eq. 50); the actual transpiration rate from the unsaturated store ${\sum}_{n=\mathrm{1}}^{m}{E}_{\mathrm{trans},\phantom{\rule{0.125em}{0ex}}\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}n}^{t}$ for m unsaturated soil layers; ${S}_{\mathrm{sat}}^{t}$ (Eq. 72) and the unsaturated store capacity ${S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}\mathrm{max}}^{t}$, which is based on ${S}_{\mathrm{sat}}^{t}$ (Eq. 72); and the sum of unsaturated storage ${\sum}_{n=\mathrm{1}}^{m}{S}_{\mathrm{unsat},\phantom{\rule{0.125em}{0ex}}n}^{t}$ for m unsaturated soil layers (soil water balance check after Eq. 78):
with z_{soil}, θ_{s} and θ_{r} as previously defined. Then, the maximum capillary rate is scaled using the following empirical equation (e.g. Zammouri, 2001; Yang et al., 2011; Wang et al., 2016):
where ${C}_{\mathrm{act}}^{t}$ [mm t^{−1}] is the capillary rate at time step t; z_{cap, maxdepth} [mm] is the critical water depth beyond which capillary rise ceases; n_{cap} [–] is an empirical coefficient related to soil properties and climate, generally set between 1–3; and ${z}_{\mathrm{watertable}}^{t\mathrm{1}}$ and z_{rooting} are as previously defined. When the soil bucket in wflow_sbm is split up into different layers, ${C}_{\mathrm{act}}^{t}$ is divided over the different unsaturated soil layers, from the bottom to the top unsaturated soil layer, without exceeding the saturated water content θ_{s}.
2.4.4 Leakage
In wflow_sbm it is possible to have leakage L^{t} [mm t^{−1}] at time step t from the saturated store ${S}_{\mathrm{sat}}^{t}$ (Eq. 72) to deeper groundwater by setting the maximum leakage model parameter L_{max} [mm t^{−1}] > 0. This water is lost from the saturated store and runs out of the model domain. L^{t} is calculated as follows:
where ${f}_{\mathrm{Kv},\phantom{\rule{0.125em}{0ex}}{n}_{\mathrm{b}}}$ is the optional multiplication factor (to correct the vertical saturated hydraulic conductivity) for the bottom soil layer n_{b} and K_{v0}, f_{Kv} and z_{soil} are as previously defined.
2.5 Lateral subsurface flow
In wflow_sbm the kinematicwave approach is used to route subsurface flow laterally. The saturated store can be drained laterally by saturated downslope subsurface flow for a slope with width w [m] according to
where c_{land slope} is the land slope [m m^{−1}], Q_{subsurface} is subsurface flow [m^{3} d^{−1}], ${K}_{\mathrm{h}\mathrm{0}}=\mathrm{0.001}{K}_{\mathrm{v}\mathrm{0}}{f}_{\mathrm{Kh}\mathrm{0}}\frac{\mathrm{\Delta}{t}_{\mathrm{b}}}{\mathrm{\Delta}t}$ is the horizontal saturated hydraulic conductivity at the soil surface [m d^{−1}] based on the vertical saturated hydraulic conductivity at the soil surface K_{v0} [mm t^{−1}] and an optional multiplication factor f_{Kh0} [–] (default: 100) applied to K_{v0}, z_{ssf, watertable} [m] is the water table depth (set by z_{watertable} [mm] after unit conversion at the start of the lateralsubsurfaceflow computation), z_{ssf, soil} [m] is the soil depth (set by z_{soil} [mm] after unit conversion), f_{ssf, Kv} is a scaling parameter [m^{−1}] (set by f_{Kv} [mm^{−1}] after unit conversion), and Δt_{b} and Δt are as previously defined. This is combined with the following continuity equation:
where h is the water table height [m]; x is the distance downslope [m]; w is the flow width [m]; $R=\mathrm{0.001}\frac{\mathrm{\Delta}{t}_{\mathrm{b}}}{\mathrm{\Delta}t}{R}_{\mathrm{input}}$ is the net input rate [m d^{−1}] computed from the net input rate variable R_{input} [mm t^{−1}], which is part of the vertical SBM concept; and θ_{s} and θ_{r} are as previously defined. Substituting for $h\left(\frac{\mathit{\delta}{Q}_{\mathrm{subsurface}}}{\mathit{\delta}h}\right)$ gives
The kinematicwave equation for lateral subsurface flow is solved iteratively using Newton's method. In wflow_sbm, the flow width w is calculated for each grid cell by dividing the cell area with the distance downslope x, based on the flow direction and the grid cell dimensions. The land slope c_{land slope} needs to be provided for each grid cell in wflow_sbm. The net input rate R_{input} in wflow_sbm consists of transfer of soil water ${Q}_{\mathrm{transfer},\phantom{\rule{0.125em}{0ex}}\mathrm{act},\phantom{\rule{0.125em}{0ex}}m}^{t}$ from unsaturated soil layer m above the water table ${z}_{\mathrm{watertable}}^{t\mathrm{1}}$ [mm] at the previous time step and the losses through capillary rise ${C}_{\mathrm{act}}^{t}$, transpiration ${E}_{\mathrm{trans},\phantom{\rule{0.125em}{0ex}}\mathrm{sat}}^{t}$ from the saturated store, leakage L^{t} and soil evaporation ${E}_{\mathrm{act},\phantom{\rule{0.125em}{0ex}}\mathrm{soil},\phantom{\rule{0.125em}{0ex}}\mathrm{sat}}^{t}$ from the saturated store, with the unit mm t^{−1} for these terms. After the lateralsubsurfaceflow calculation, which is bounded by the maximum lateral subsurface flow rate based on z_{ssf, soil}, a check is made to determine if saturation of the entire soil column occurs and as a consequence saturation excess overland flow is triggered. Water exfiltrating under saturated conditions ${R}_{\mathrm{exfilt},\phantom{\rule{0.125em}{0ex}}\mathrm{sat}}^{t}$ [m t^{−1}] at time step t is calculated as follows:
where $\mathrm{\Delta}{S}_{\mathrm{subsurface}}^{t}$ [m] is the change in subsurface storage at time step t; ${Q}_{\mathrm{subsurface},\phantom{\rule{0.125em}{0ex}}\mathrm{in}}^{t}$ [m^{3} d^{−1}] is the subsurface flow into a cell at time step t; ${Q}_{\mathrm{subsurface},\phantom{\rule{0.125em}{0ex}}\mathrm{out}}^{t}$ [m^{3} d^{−1}] is the subsurface flow out of a cell at time step t; ${z}_{\mathrm{ssf},\phantom{\rule{0.125em}{0ex}}\mathrm{watertable}}^{t\mathrm{1}}$ [m] is the water table depth at the previous time step; Δt_{ssf} [d] is the time increment (same length as Δt [s], expressed in different units [d]) of the lateral subsurface flow component; and R_{input}, w, x, θ_{s} and θ_{r} are as previously defined. Additionally, after the lateralsubsurfaceflow calculation, wflow_sbm checks if exfiltration ${R}_{\mathrm{exfilt},\phantom{\rule{0.125em}{0ex}}\mathrm{unsat}}^{t}$ [mm t^{−1}] of the unsaturated store onto the land surface occurs because of a change in water table depth z_{watertable} [mm] (set by ${z}_{\mathrm{ssf},\phantom{\rule{0.125em}{0ex}}\mathrm{watertable}}^{t}$ after unit conversion). This check is performed from the bottom unsaturated layer (at the previous time step) to the top unsaturated layer, where the excess of unsaturated storage for each layer is transferred from the bottom to the top unsaturated layer, and can result in exfiltration of water onto the land surface.
2.6 Surface flow routing
The kinematicwave approach is used for river and overland flow routing. The kinematicwave equations (Chow et al., 1988) are
and they can be combined as
where Q is the surface flow in the kinematic wave [m^{3} s^{−1}], x is the length of the flow pathway [m], A is the crosssection area of the flow pathway [m^{2}], Q_{inflow} is the lateral inflow per unit length into the kinematic wave [m^{2} s^{−1}], and α and β are coefficients. These coefficients can be determined using Manning's equation (Chow et al., 1988), resulting in
where P_{w} [m] is the wetted perimeter, c_{slope} (c_{land slope} for overland flow and c_{river slope} for river flow) is the slope [m m^{−1}] and n (n_{land} for overland flow and n_{river} for river flow) is Manning's coefficient [s m${}^{\mathrm{1}/\mathrm{3}}$]. The wetted perimeter P_{w} for river flow is calculated by adding the river width (w_{river}) and 2 times half of the river bankfull depth (h_{bankfull}). For overland flow, P_{w} is set equal to the effective flow width, determined by dividing the grid cell area by the flow length and subtracting w_{river}. In wflow_sbm for river flow the parameters w_{river}, length (x_{river}) and c_{river slope} need to be provided, and for overland flow c_{land slope} needs to be provided. The lateral inflow per unit flow length for overland flow routing consists of infiltration excess water ${F}_{\mathrm{excess}}^{t}$, saturation excess water during infiltration ${F}_{\mathrm{excess},\phantom{\rule{0.125em}{0ex}}\mathrm{sat}}^{t}$, exfiltration water from the unsaturated store ${R}_{\mathrm{exfilt},\phantom{\rule{0.125em}{0ex}}\mathrm{unsat}}^{t}$, water exfiltrating under saturated conditions ${R}_{\mathrm{exfilt},\phantom{\rule{0.125em}{0ex}}\mathrm{sat}}^{t}$, runoff from open water ${R}_{\text{open water}}^{t}$ and openwater evaporation loss ${E}_{\text{open water}}^{t}$, converted to m^{2} s^{−1}. The lateral inflow per unit length of x_{river} for river flow routing consists of overland flow [m^{3} s^{−1}], lateral subsurface flow [m^{3} d^{−1}], runoff from the river ${R}_{\mathrm{river}}^{t}$ [mm t^{−1}] and river evaporation loss ${E}_{\mathrm{river}}^{t}$ [mm t^{−1}], converted to m^{2} s^{−1}. Like the lateral subsurface routing, these equations are solved in wflow_sbm using Newton's method. The number of iterations for surface flow in the kinematic wave at time step t of length Δt [s] defaults to the Courant number C:
where c_{k} [m s^{−1}] is the kinematicwave celerity, ${c}_{\mathrm{k}}=\frac{\mathrm{1}}{\mathit{\alpha}\mathit{\beta}{Q}^{\mathit{\beta}\mathrm{1}}}$, and Δx [m] is the space increment. The number of iterations within a time step t is calculated by multiplying the 95th percentile of C (to remove potential very high values (outliers)) for the wflow_sbm model domain with 1.25. The number of iterations can also be fixed to a specific subtime step [s] for both river and overland flow; this is a model setting in the wflow_sbm configuration file. For river cells in wflow_sbm, where overland and river flow can both be present, lateral subsurface and overland flow into the river cell is partitioned based on the land slope of the river cell c_{land slope, river} [m m^{−1}] and the land slope c_{land slope, upstream} [m m^{−1}] of the upstream cell:
where f_{to river} [–] is the fraction of lateral subsurface or overland flow from an upstream cell that flows into the river and f_{to land} [–] is the fraction of lateral subsurface or overland flow from an upstream cell that flows into the downstream kinematic reservoir of lateral subsurface and overland flow, respectively. In the case where a river cell has the same flow direction as the upstream cell, f_{to river}=0, and thus overland and lateral subsurface flow from the upstream cell does not contribute to flow into the river.
2.7 Reservoirs and lakes
2.7.1 Reservoirs
In wflow_sbm, reservoirs can be included in the kinematicwave routing for river flow. The first step in the reservoir module is to calculate the storage ${S}_{\mathrm{res}}^{{t}_{i}}$ [m^{3}] at kinematicwave time step t_{i} with length Δt_{i} [s], based on the storage ${S}_{\mathrm{res}}^{{t}_{i}\mathrm{1}}$ at the previous time step of the kinematic wave, inflow ${Q}_{\mathrm{in},\phantom{\rule{0.125em}{0ex}}\mathrm{res}}^{{t}_{i}}$ [m^{3} s^{−1}] at time step t_{i}, average precipitation ${P}_{\mathrm{res}}^{{t}_{i}}$ [mm ${t}_{i}^{\mathrm{1}}$] (converted from ${P}_{\mathrm{res}}^{t}$ [mm t^{−1}]) and potential evapotranspiration ${E}_{\mathrm{pot},\phantom{\rule{0.125em}{0ex}}\mathrm{res}}^{{t}_{i}}$ [mm ${t}_{i}^{\mathrm{1}}$] (converted from ${E}_{\mathrm{pot},\phantom{\rule{0.125em}{0ex}}\mathrm{res}}^{t}$ [mm t^{−1}]) on the reservoir area A_{res} [m^{2}] at time step t_{i}:
Then the storage fraction ${f}_{\mathrm{res},\phantom{\rule{0.125em}{0ex}}\mathrm{storage}}^{{t}_{i}}$ [–] is calculated based on the maximum storage of the reservoir S_{res, max} [m^{3}] (above this storage amount water is spilled):
The minimum release ${R}_{\mathrm{min}}^{{t}_{i}}$ [m^{3} ${t}_{i}^{\mathrm{1}}$] at the kinematicwave time step t_{i} is based on a sigmoid function, the minimum flow requirement downstream of the reservoir Q_{min req.} [m^{3} s^{−1}], the target minimum storage fraction (of S_{res, max}) f_{res, min} [–] and ${f}_{\mathrm{res},\phantom{\rule{0.125em}{0ex}}\mathrm{storage}}^{{t}_{i}}$ at time step t_{i},
and ${R}_{\mathrm{min}}^{{t}_{i}}$ is subtracted from the reservoir storage ${S}_{\mathrm{res}}^{{t}_{i}}$:
An additional release ${R}^{{t}_{i}}$ [m^{3} ${t}_{i}^{\mathrm{1}}$] occurs when the reservoir storage is above the target maximum storage fraction f_{res, max} [–], controlled by the maximum release capacity below the spillway Q_{max, res} [m^{3} s^{−1}]:
where ${R}_{\mathrm{pot}}^{{t}_{i}}$ [m^{3} ${t}_{i}^{\mathrm{1}}$] is the potential reservoir release and ${S}_{\text{res, above max}}^{{t}_{i}}$ [m^{3}] is the reservoir storage above the maximum reservoir storage. ${R}^{{t}_{i}}$ is subtracted from the reservoir storage ${S}_{\mathrm{res}}^{{t}_{i}}$:
The total reservoir inflow ${Q}_{\mathrm{in},\phantom{\rule{0.125em}{0ex}}\mathrm{res}}^{t}$ [m^{3} t^{−1}] and outflow ${Q}_{\mathrm{out},\phantom{\rule{0.125em}{0ex}}\mathrm{res}}^{t}$ [m^{3} t^{−1}] for model time step t of length Δt [s] are calculated as follows:
where n_{i} refers to the number of iterations within model time step t and ${Q}_{\mathrm{in},\phantom{\rule{0.125em}{0ex}}\mathrm{res}}^{{t}_{i}}$, ${R}_{\mathrm{min}}^{{t}_{i}}$, ${R}^{{t}_{i}}$ and Δt_{i} are as previously defined.
2.7.2 Lakes
Lakes can be included in the kinematicwave routing for river flow, and as with the reservoirs in wflow_sbm, a mass balance approach is used for modelling lakes:
where S_{lake} is lake storage [m^{3}], t_{i} is the kinematicwave time step of length Δt_{i} [s], Q_{in, lake} is the sum of inflows [m^{3} s^{−1}], Q_{out, lake} is the lake outflow at the outlet [m^{3} s^{−1}], ${P}_{\mathrm{lake}}^{{t}_{i}+\mathrm{\Delta}{t}_{i}}$ is the precipitation amount [mm] during Δt_{i} (converted from ${P}_{\mathrm{lake}}^{t}$ [mm t^{−1}]), ${E}_{\mathrm{lake}}^{{t}_{i}+\mathrm{\Delta}{t}_{i}}$ is lake evaporation [mm] during Δt_{i} (converted from ${E}_{\mathrm{lake}}^{t}$ [mm t^{−1}]) and A_{lake} is the lake surface [m^{2}]. Most of the terms in Eq. (108) are known at the current or previous time step, except ${S}_{\mathrm{lake}}^{{t}_{i}+\mathrm{\Delta}{t}_{i}}$ and ${Q}_{\mathrm{out},\phantom{\rule{0.125em}{0ex}}\mathrm{lake}}^{{t}_{i}+\mathrm{\Delta}{t}_{i}}$. For lakes characterized by a storage curve of the form S_{lake}=A_{lake}H_{lake} and the rating curve
where H_{0, lake} is the minimum water level under which the outflow is zero, α_{lake} [m s^{−1}] is a parameter that depends on lake outlet characteristics and the β_{lake} exponent has a value of 2 (parabolic weir), the modified Puls approach is used. Then, S_{lake} can be expressed as follows:
where $h={H}_{\mathrm{lake}}{H}_{\mathrm{0},\phantom{\rule{0.125em}{0ex}}\mathrm{lake}}$. Inserting this equation in the mass balance equation gives
The solution for ${Q}_{\mathrm{out},\phantom{\rule{0.125em}{0ex}}\mathrm{lake}}^{{t}_{i}+\mathrm{\Delta}{t}_{i}}$ is then
where ${f}_{\mathrm{lake}}=\frac{{A}_{\mathrm{lake}}}{\mathrm{\Delta}{t}_{i}\sqrt{{\mathit{\alpha}}_{\mathrm{lake}}}}$ and ${\mathrm{SI}}_{\mathrm{lake}}=\frac{{S}_{\mathrm{lake}}^{{t}_{i}}}{\mathrm{\Delta}{t}_{i}}+{Q}_{\mathrm{in},\phantom{\rule{0.125em}{0ex}}\mathrm{lake}}^{{t}_{i}+\mathrm{\Delta}{t}_{i}}+\frac{\mathrm{0.001}({P}_{\mathrm{lake}}^{{t}_{i}+\mathrm{\Delta}{t}_{i}}{E}_{\mathrm{lake}}^{{t}_{i}+\mathrm{\Delta}{t}_{i}}){A}_{\mathrm{lake}}}{\mathrm{\Delta}{t}_{i}}$.
The modified Puls approach is not applicable for lakes characterized by a rating curve (Eq. 109) with β_{lake}≠2 (nonparabolic weir; for a rectangular weir, usually a value of $\mathrm{3}/\mathrm{2}$ is used) or a rating curve from measurements (linear interpolation of Q_{out, lake} and H_{lake} values in a lookup table) in combination with a storage curve from measurements (linear interpolation of S_{lake} and H_{lake} values in a lookup table) or computed from the relationship S_{lake}=A_{lake}H_{lake}. For these lakes Q_{out, lake} is first computed for each time step based on H_{lake} at the previous time step. Then, S_{lake} is updated with Eq. (108), and H_{lake} is updated with the storage curve based on the updated S_{lake}. For nearby lakes which are connected, it is possible to link the lakes and return flow can be allowed from the downstream to the upstream lake. An average lake water level (H_{lake, avg} [m]) should be provided as an initial state when wflow_sbm is initialized with default values in the code (cold start); see also Table A2. The total lake inflow ${Q}_{\mathrm{in},\phantom{\rule{0.125em}{0ex}}\mathrm{lake}}^{t}$ [m^{3} t^{−1}] and outflow ${Q}_{\mathrm{out},\phantom{\rule{0.125em}{0ex}}\mathrm{lake}}^{t}$ [m^{3} t^{−1}] for model time step t are calculated as follows:
where n_{i} refers to the number of iterations within model time step t and ${Q}_{\mathrm{in},\phantom{\rule{0.125em}{0ex}}\mathrm{lake}}^{{t}_{i}}$, ${Q}_{\mathrm{out},\phantom{\rule{0.125em}{0ex}}\mathrm{lake}}^{{t}_{i}}$ and Δt_{i} are as previously defined.
One of the reasons to switch to the Julia programming language is that it offers high performance, which is required for largescale highresolution hydrological model applications. Here we compare the simulation times of wflow_sbm between the Julia (van Verseveld et al., 2024) and Python (Schellekens et al., 2020) versions for three large catchments: the Moselle, Meuse and Rhine (Fig. 2). We used HydroMTWflow (Eilander et al., 2022) to set up the models for the three catchments at a resolution of 30^{′′} (∼ 1 km × 1 km). The models were run at a daily time step for 5 years (2000–2005) with ERA5 forcing data. We excluded the input/output (I/O) operations to allow for a clean comparison between the Julia and Python version and ran the simulation on a machine with the following specifications: a desktop with an Intel Xeon Gold 6144 CPU (with four cores, four threads exposed to the user) and 16 GB RAM.
The switch to Julia results in substantially smaller simulation times, independent of the size of the catchment (Fig. 2). By enabling threads (Julia version), the simulation times decrease further, leading to a model that runs 4–5 times faster compared to the Python version. For the Rhine catchment, the simulation time for 1 year is 120 min for the wflow_sbm Python version; for the Julia version, this is 36 and 23 min with one and four threads, respectively. These simulation times take up most of the total computational time (e.g. ∼ 98 % for the Rhine model with Julia running on one and multiple threads). These results show that the wflow_sbm Julia version is suitable for largescale highresolution model applications.
The wflow_sbm model has been applied to a number of specific cases. Below we describe these specific applications and the a priori parameter estimation, including forcing, with HydroMTWflow for a variety of hydroclimates and hydrological processes.
4.1 Parameter estimation with HydroMTWflow
The estimation of wflow_sbm model parameters is based on earlier work by Imhoff et al. (2020) that focused on the Rhine Basin and on the development of the Iterative Hydrography Upscaling (IHU) method by Eilander et al. (2021) to derive flow direction and representative river length, slope and width parameters. Eilander et al. (2021) showed an improved accuracy with IHUupscaled flow direction maps, applied to MERIT Hydro, compared to other oftenused upscaling methods. Furthermore, for a case study of the Rhine Basin, Eilander et al. (2021) showed that with IHU applied to MERIT Hydro, errors in the timing and magnitude of simulated peak discharge compared to simulations at the native data resolution are minimized. Imhoff et al. (2020) used available pointscale (pedo)transfer functions (PTFs) from the literature to generate seamless parameter maps for the Rhine Basin. Following a multiscale parameter regionalization (MPR) technique (Samaniego et al., 2010), parameters were estimated at the original data resolution (“level 0”), and upscaled to the model resolution (“level 1”) with upscaling operators. Although universal scaling rules for hydrological model parameters are not available, the correct upscaling operator is found when model parameters are characterized by a constant mean and standard deviation across different spatial resolutions. Additionally, model fluxes and states should be consistent across different spatial resolutions. For the Rhine Basin, Imhoff et al. (2020) found that modelled actual evapotranspiration fluxes were consistent across different spatial resolutions (1.2, 2.4, 3.6 and 4.8 km). Routed discharge in headwater basins was not consistent across scales (KGE decreased from the finest to the coarsest resolution), while for the main Rhine River, routed discharge was consistent. For recharge fluxes, relatively large differences were found for regions with high drainage densities.
The transfer functions and upscaling operators to derive wflow_sbm model parameters for any region in the world are part of the HydroMTWflow software (Eilander et al., 2022) and are listed in Table 1. For two sensitive wflow_sbm model parameters, the temperature threshold (s_{fall, T threshold}) and the multiplication factor f_{Kh0}, a PTF is not available (Imhoff et al., 2020). For s_{fall, T threshold} and f_{Kh0}, a uniform default value of 0.0 °C and 100.0 is applied (Table 2), respectively. The a priori parameter estimation for wflow_sbm provides a model setup without the need for much further calibration; in most cases only the model parameter f_{Kh0} is tuned (e.g. Wannasin et al., 2021a; Sperna Weiland et al., 2021).
Table 1 also includes references to examples of global datasets that can be used to set up a wflow_sbm model with HydroMTWflow. For soil properties, SoilGrids (Hengl et al., 2017) at 250 m resolution is available. For land cover the datasets GlobCover 2009 (Arino et al., 2010) at 300 m resolution, VITO v2.0.2 (Buchhorn et al., 2019) at 100 m resolution and CORINE Land Cover (CLC) 2018 (European Environment Agency, 2018) at 100 m resolution are currently available. Leaf area index climatology is based on the MCD15A3H Version6 leaf area index product, at 500 m resolution (Myneni et al., 2015). For elevation and derived data, the MERITDEMbased hydrography map (MERIT Hydro; Yamazaki et al., 2019) at 90 m resolution is used. It contains a global flow direction map derived from the MultiErrorRemoved ImprovedTerrain DEM dataset (MERIT DEM; Yamazaki et al., 2017) and a synthetic water layer map that consists of a combination of the Global “1s” Water Body Map (G1WBM; Yamazaki et al., 2015), Global Surface Water Occurrence (GSWO; Pekel et al., 2016) and waterrelated features from OpenStreetMap. The fineresolution MERIT Hydro flow direction map is upscaled to the wflow_sbm model resolution with the Iterative Hydrography Upscaling (IHU) method (Eilander et al., 2021). River width (w_{river}) and bankfull depth (h_{bankfull}) are based on MERIT Hydro (river mask based on a minimum upstream area) and the global reachlevel bankfull river width dataset from Lin et al. (2019). River bankfull depth h_{bankfull} is estimated from bankfull discharge (Q_{bankfull}) data in Lin et al. (2019) with the following power law relationship:
with c=0.27 (default) and p=0.30 (default).
For glacierrelated model parameters, the Randolph Glacier Inventory 6.0 (Pfeffer et al., 2014) is available. Lakerelated parameters are derived from the HydroLAKES Version 1.0 (Messager et al., 2016) dataset, and reservoir parameters are based on a combination of the Global Reservoir and Dam Database (GRanD), Version 1, Revision 01 (v1.01) (Lehner et al., 2011); HydroLAKES Version 1.0 (Messager et al., 2016); and GSWO. For more details on the setup of a wflow_sbm model, from global (or regional/local) datasets with HydroMTWflow, we refer to the documentation and code of HydroMTWflow (Eilander et al., 2022).
Figure 3 shows examples of model parameter maps set up with HydroMTWflow for the Moselle catchment (see also Sect. 4.2.4). The parameters are related to elevation (slope), elevation and flow direction (Strahler stream order), land cover (rooting depth), and soil properties (vertical saturated hydraulic conductivity).
Rawls and Brakensiek (1989)van Dijk and Bruijnzeel (2001)Brakensiek et al. (1984)Myneni et al. (2015)Engman (1986)Kilgore (1997)Liu et al. (2005)Schenk and Jackson (2002)Fan et al. (2016)Pitman (1989)Liu (1998)Horn (1981)Lin et al. (2019)Lin et al. (2019)Hengl et al. (2017)ESDAC (2004)Pitman (1989)Liu (1998)Tóth et al. (2015)Tóth et al. (2015)Pfeffer et al. (2014)Messager et al. (2016)Lehner et al. (2011)Messager et al. (2016)Pekel et al. (2016)4.2 Wflow_sbm model cases
The wflow_sbm model cases have been set up with HydroMTWflow at a resolution of 30^{′′} based on the (pedo)transfer functions listed in Table 1 and HydroMTWflow constant default model parameters listed in Table 2. We present two model cases illustrating the model sensitivity to the model parameter f_{Kh0}, a multiplication factor applied to the vertical saturated hydraulic conductivity at the soil surface K_{v0} to calculate the horizontal saturated hydraulic conductivity (Whanganui catchment; see Sect. 4.2.1), and the parameter f_{Kv} (Crystal River catchment; see Sect. 4.2.2) that controls the exponential decline in vertical saturated hydraulic conductivity. We also illustrate how wflow_sbm performs based on the a priori parameter estimation with only changing the model parameter f_{Kh0} for the Umeälven catchment (Sect. 4.2.3), where reservoir operations and snow processes play an important role; for the Moselle catchment (Sect. 4.2.4) including discharge and catchmentaverage soil moisture as output; and for the Rhine catchment (Sect. 4.2.5) to demonstrate the ability of wflow_sbm to represent the spatial distribution of actual evapotranspiration E_{a} and snow storage S_{snow}. Finally, we present a model case for the Ouémé catchment (Sect. 4.2.6), where groundwater loss plays an important role. The locations of the wflow_sbm model cases on a global map are shown in Fig. 4a. For each model case, a map of the catchment with elevation and rivers is presented in Fig. 4b–f. For each model case, four soil layers are defined as (default) 0–100, 100–400, 400–1200 and 1200 mm up to the maximum soil depth z_{soil} (Table 1) to capture changes in soil hydraulic properties and roots, and thus soil moisture fluxes, with depth. Wflow_sbm determines the actual soil layer thickness for each layer per grid cell based on z_{soil}. For river and overland flow, the time step is set to a fixed value of 900 and 3600 s, respectively. When snow is enabled, a reduction factor to infiltration in soils because of frozen conditions is not applied. Other more specific model settings are described per model case in Sect. 4.2.1–4.2.6.
The model performance of the wflow_sbm applications is here assessed with the modified Kling–Gupta efficiency (KGE; Kling et al., 2012) for discharge and catchmentaverage soil moisture as
where r is the correlation coefficient between observed and simulated values, β is a bias term, and γ is the variability ratio:
where μ_{sim} is the mean of simulated values, μ_{obs} is the mean of observed values, σ_{sim} is the standard deviation of simulated values and σ_{obs} is the standard deviation of observed values. KGE = 1 means a perfect match between simulated and observed values, and KGE ≈ −0.41 indicates the model simulation is as accurate as the observed mean (Knoben et al., 2019). For the assessment of the reproduction of spatial patterns for the model case of the Rhine catchment, the spatial efficiency metric (E_{SP}; Dembélé et al., 2020) is used:
where r_{s} is the Spearman rankorder correlation coefficient; d_{i} is the difference between the two ranks (observed variable X_{obs} and simulated variable X_{sim}) of each cell i in n grid cells; α is a term that determines the spatial location matching, calculated as the root mean squared error (E_{RMS}) of the standardized values (z scores, Z_{X}) of X_{sim} and X_{obs}; and γ is the variability ratio (equal to the γ term of KGE). As with KGE, E_{SP} ranges from −∞ to 1 (a perfect match between X_{sim} and X_{obs}).
For each model case, we use the first year as a warmup period. For the simulations of all but three model cases, we use ERA5 forcing, temperature and potential evapotranspiration (using the de Bruin method; de Bruin et al., 2016), which are derived based on downscaled ERA5 fields using a fixed lapse rate of −0.0065 °C m^{−1}. For the model case of the Crystal River catchment (Sect. 4.2.2), forcing is based on the dataset by Maurer et al. (2002); for the model case of the Rhine catchment (Sect. 4.2.5), we use the MultiSource WeightedEnsemble Precipitation version2.8 (MSWEP V2.8) (Beck et al., 2019) global precipitation product instead of ERA5 rainfall; and for the case of the Ouémé catchment (Sect. 4.2.6), we use Climate Hazards group Infrared Precipitation with Stations (CHIRPS) rainfall (Funk et al., 2015) estimates instead of ERA5 rainfall.
4.2.1 New Zealand – Whanganui River – effect of model parameter f_{Kh0}
The wflow_sbm model parameter f_{Kh0} is a multiplication factor applied to the vertical saturated hydraulic conductivity at the soil surface K_{v0} to calculate the horizontal saturated hydraulic conductivity used for computing lateral subsurface flow. This parameter compensates for anisotropy, smallscale saturated hydraulic conductivity (soil core) measurements that do not represent largerscale hydraulic conductivity and smaller flow length scales (hillslope) in reality not represented by the model resolution. Landcoverderived model parameters are based on VITO v2.0.2 (Buchhorn et al., 2019). For this model case, the snow (including the snow avalanche routine) and glacier model are enabled. To illustrate the effect of different f_{Kh0} values (1, 20, 50, 100), Fig. 5 shows the discharge simulation (daily time step) and KGE values for Global Runoff Data Centre (GRDC) station 5865600 of the Whanganui River catchment in New Zealand, for the year 1996. Figure 5 clearly shows that higher f_{Kh0} values generally result in higher baseflow values and flattened peaks. The f_{Kh0} value of 20 results in the highest KGE of 0.71 for the year 1996. For the complete period of simulation 1979–2009, the KGE values were 0.63, 0.79, 0.68 and 0.55 for f_{Kh0} values of 1, 20, 50 and 100, respectively.
By changing the parameter f_{Kh0}, it is expected that the contribution of overland flow and lateral subsurface flow to river discharge will change. We show in Fig. 6a the effect of different f_{Kh0} values (1 and 100) on the average lateral inflow components subsurface flow Q_{subsurface, to river} and overland flow Q_{land, to river} to the river. With an f_{Kh0} value of 1, the contribution of Q_{subsurface, to river} is minimal (maximum contribution is 0.0011 m^{3} s^{−1}), and river inflow consists mainly of Q_{land, to river}, with high values during discharge peaks that quickly drop to low values under baseflow conditions (Fig. 6a). The average water table depth z_{watertable} is low without much variation with an f_{Kh0} value of 1 (Fig. 6b). With an f_{Kh0} value of 100, the contribution of Q_{subsurface, to river} is higher and Q_{land, to river} is lower during peaks and higher under baseflow conditions compared to an f_{Kh0} value of 1. On average z_{watertable} is higher and shows more variation with an f_{Kh0} value of 100 compared to an f_{Kh0} value of 1 (Fig. 6b). Thus, f_{Kh0} controls the distribution of Q_{subsurface, to river} and Q_{land, to river}, related to the overall wetness of the catchment and the magnitude of lateral subsurface flows (Q_{subsurface}), which has an effect on the peak discharges and baseflow values of the hydrograph.
4.2.2 USA – Crystal River – effect of exponential decline in K_{v0}
The wflow_sbm parameter f_{Kv} controls the exponential decline in the vertical saturated hydraulic conductivity K_{v0} at the soil surface with depth and thus vertical flow and lateral subsurface flow. A priori, f_{Kv} is estimated with HydroMTWflow through the use of two different fitting methods using nonlinear least squares (curve_fit) and a leastsquares solution (linalg.lstsq) applied to the estimated vertical saturated hydraulic conductivity K_{vz} at different depths from SoilGrids (see also Table 1). Figure 7 shows simulated (daily time step) discharge for the Crystal River near Redstone (Colorado, USA) in 2003, with landcoverderived model parameters based on VITO v2.0.2 (Buchhorn et al., 2019). For this model case, the snow (including the snow avalanche routine) model is enabled. For both fitting methods, the performance is good (KGE 0.9 or higher); the fitting method using nonlinear least squares (curve_fit) results in a higher KGE value. This fitting method captures the rising limb (partially caused by snowmelt) during the period 15 March–15 May and generally the falling limb (less overestimation) of the hydrograph better. The average f_{Kv} value for the catchment was 0.0027 and 0.0011 for fitting using nonlinear least squares and the leastsquares solution, respectively. A higher f_{Kv} value results in a more exponential decline in K_{v0} with depth.
As with the parameter f_{Kh0}, it is expected that by changing the f_{Kv} parameter the contribution of overland flow and lateral subsurface flow to river discharge will change. We show the effect of the different f_{Kv} values as a result of different fitting methods on the average lateral inflow components Q_{subsurface, to river} and Q_{land, to river} to the river in Fig. 8a and on z_{watertable} in Fig. 8b. The curve_fit fitting method captures the rising limb of the hydrograph better during the period 15 March–15 May because the contribution of Q_{land, to river} is higher compared to the linalg.lstsq fitting method, while the Q_{subsurface, to river} contribution is similar (Fig. 8a). With the curve_fit fitting method, the catchment is wetter; z_{watertable} has a shallower pattern compared to the linalg.lstsq fitting method (Fig. 8b), causing higher Q_{land, to river} values. During the falling limb of the hydrograph, the curve_fit fitting method results in less overestimation caused by a lower Q_{subsurface, to river} contribution, while the Q_{land, to river} contribution is similar (Fig. 8a).
4.2.3 Europe – Umeälven – snow and reservoirs
The Ume River (Umeälven) is one of the largest rivers in Sweden and has been extensively cultivated for hydroelectric power. From a hydrometeorological aspect, snowfall and snowmelt play an important role here. Furthermore, accounting for these hydropower stations in the hydrological model is done on the basis of the information in the GRanD, HydroLAKES and GWSO datasets. This information is rather uncertain as hydropower operators have their own release curves or optimization schemes. Here, we simulated discharge for the period 1979 to 2019 at a daily time step, with landcoverderived model parameters based on CLC 2018 (European Environment Agency, 2018). The values derived for the vertical saturated hydraulic conductivity are quite large, on the order of a few metres per day, and hence no anisotropy factor for horizontal saturated hydraulic conductivity was applied (f_{Kh0} = 1). For this model case, the snow (including the snow avalanche routine) and glacier models are enabled and lakes and reservoirs are included.
Figure 9 presents KGE scores for Swedish Meteorological and Hydrological Institute (SMHI) stations simulated with EHYPE (Swedish Meteorological and Hydrological Institute, 2022) and wflow_sbm. Observed discharges for the stations were obtained from Swedish Meteorological and Hydrological Institute (2024). Because of differences in forcing, simulation periods and the use of KGE methods, a quantitative comparison of model performance scores is not possible. However, it is obvious that for locations (20010 and 20047) largely influenced by reservoirs and lakes, and thus reservoir operations, both models show poor performance. For a more upstream location like station 1630, EHYPE performance is good (KGE > 0.8) and wflow_sbm performance is satisfactory (KGE = 0.66). Wflow_sbm shows lower KGE values for locations 2237 and 2238, while for locations 1733 and 1734 wflow_sbm shows a better performance. Wflow_sbm performance for the unregulated Vindel River (locations 1630, 2237 and 2238), a tributary to the Ume River, could probably be improved on the basis of EHYPE performance. This could be done for example through a further analysis of the impact of different snow model parameters like s_{fall, T threshold} and s_{ddf} on model performance, lake model settings of Lake Storovan, and possibly forcing datasets other than ERA5.
Figure 10 shows simulated and observed discharge for the SMHI stations within the Umeälven catchment for the year 1993. The performance for the most downstream stations 1733 and 1734 is satisfactory. The underestimation of discharge peaks for station 1733 during July and August is mostly caused by underestimating discharge peaks at station 20010 (downstream of Lake Storuman) during the same period. The actual release scheme of Lake Storuman is very likely not captured well enough by wflow_sbm. Wflow_sbm overestimates baseflow and underestimates peak flows for station 2238 (Fig. 10), downstream of Lake Storovan. This may be caused for example by the lake model settings of Lake Storovan or upstream model parameters related to lateral subsurface flow (see also Sect. 4.2.2 and 4.2.1) or snow dynamics, and it has an effect on the underestimation of peak flows at the downstream station 1734 by wflow_sbm. Overall, we show that with the a priori parameter estimation, the performance for the Umeälven catchment is (qualitatively) similar to the EHYPE performance for this catchment.
4.2.4 Europe – Moselle – soil moisture
Besides comparing simulated discharge with observed discharge, it can be useful to compare model state variables like soil moisture or snow water equivalent to actual observations or satellitebased datasets for the purposes of calibration and validation of the hydrological model. Here we perform simulation for the period 1979 to 2019 at a 6hourly time step, with landcoverderived model parameters that are based on CLC 2018 (European Environment Agency, 2018). We use an f_{Kh0} value of 250 based on previous hydrological modelling work of the Rhine Basin in Imhoff et al. (2020). For this model case, the snow (including the snow avalanche routine) model is enabled and lakes and reservoirs are included. The simulated soil moisture dynamics of the first soil layer (0–10 cm) are compared to the SMAP Enhanced L3 Radiometer dataset (O'Neill et al., 2021), averaged over the catchment, for the period 2015 to 2019. Figure 11 shows that wflow_sbm captures the average soil moisture dynamics quite well, with a KGE score of 0.83. Some of the lower and higher soil moisture values of the SMAP Enhanced L3 Radiometer dataset are not captured by wflow_sbm. This could be caused by using the default first soil layer thickness of 10 cm here, while the SMAP Enhanced L3 Radiometer dataset represents the top 5 cm of the soil column. Additionally, difference between wflow_sbm and the SMAP Enhanced L3 Radiometer in the saturated and residual water contents used can also play a role; for example the average saturated and residual water content for the wflow_sbm model is 0.44 and 0.17, respectively, while the SMAP soil moisture product shows values outside of this range (Fig. 11). Figure 12 shows simulated and observed daily discharge for GRDC station 6336050 (Cochem, near the outlet of the Moselle River into the Rhine River), for the period 2007–2008, with a KGE score of 0.74. The KGE score for GRDC station 6336050 for the complete simulation period is 0.71, and thus the overall performance of simulated soil moisture dynamics and discharge at the catchment scale is good.
4.2.5 Europe – Rhine – actual evapotranspiration and snow storage
As a spatially distributed hydrological model, wflow_sbm can easily make direct use of spatial datasets for model calibration, evaluation and data assimilation. Because of the increasing availability of satellitebased earth observations, also at finer temporal and spatial resolutions, hydrological modelling studies increasingly make use of these datasets for calibration, evaluation and dataassimilation purposes (e.g. López López et al., 2016; Demirel et al., 2018; Dembélé et al., 2020). Here, we demonstrate the ability of wflow_sbm to represent the spatial distribution of actual evapotranspiration E_{a} and snow storage S_{snow}, with the a priori estimation of the model parameters and only changing the f_{Kh0} parameter, for the Rhine catchment. Simulation is for the period 2014 to 2019 at a daily time step, with landcoverderived model parameters that are based on CLC 2018 (European Environment Agency, 2018). We use a regionally optimized f_{Kh0} map, initially based on hydrological modelling work of the Rhine Basin in Imhoff et al. (2020) with further improvements as part of hydrological modelling work for Rijkswaterstaat (part of the Dutch Ministry of Infrastructure and Water Management). For this model case, the snow (including the snow avalanche routine) and glacier models are enabled and lakes and reservoirs are included. The KGE score for GRDC station 6435060 (Lobith) at the Dutch–German border for the complete simulation period is 0.85. For the actual evapotranspiration comparison, we use The Global Land Evaporation Amsterdam Model (GLEAM) v3.7b daily actual evapotranspiration data with a spatial resolution of ∼ 28 km (Martens et al., 2017; Miralles et al., 2011) for the period 2015 to 2019. For the forcing variable precipitation, GLEAM and wflow_sbm both use the global precipitation product MSWEP V2.8. Wflow_sbm E_{a} is upscaled to the GLEAM resolution using average resampling. Figure 13 shows the average annual E_{a} for the period 2015–2019 for wflow_sbm and GLEAM. The spatial variability is quite similar (γ=1.07) and the spatial correlation is moderate (r_{s}=0.70), but the matching of the spatial location of grid cells (α=0.19) is not satisfactory, leading to an E_{SP} of 0.13. GLEAM and wflow_sbm both show a higher region of E_{a} in northern Switzerland, which is characterized by lower elevations and flatter terrain than the southern part of Switzerland. For other parts in the Rhine catchment GLEAM shows spatial clusters with higher E_{a} values that are not matched by wflow_sbm. While an E_{SP} value of 0.13 may seem low, E_{SP} can be considered a tough criterion. Koch et al. (2018) identified the SPAtial EFficiency metric (SPAEF) as such, a spatial efficiency metric equivalent to E_{SP}. Additionally, according to Koch et al. (2018) more detailed investigation into the relationship between spatial variability and spatial efficiency metrics such as SPAEF and E_{SP} is required to be able to put these metrics into context, for example for the comparison of different models or catchments. The temporal KGE scores for each grid cell (not shown) show a satisfactory to good performance; the median KGE is 0.79, with a range between 0.54 and 0.91, indicating wflow_sbm can represent E_{a} from GLEAM. Comparing catchmentaverage daily E_{a} from GLEAM and wflow_sbm (Fig. 14) shows a good agreement (KGE = 0.87), with generally an underestimation of E_{a} by wflow_sbm during the months July and August, and GLEAM showing more variability during autumn and winter. Overall, these results indicate that wflow_sbm can represent daily E_{a} variability from GLEAM, although the matching of the spatial location of grid cells is not well represented. This might be caused by, for example, the difference in resolution between GLEAM and wflow_sbm (especially in regions with complex terrain), the potentialevaporation method used by each model (GLEAM uses the Priestley–Taylor equation), irrigation (not (yet) included in wflow_sbm, while in GLEAM this is partly corrected for by the assimilation of satellite soil moisture; Miralles et al., 2011), and differences in snow evaporation approaches between wflow_sbm (only evaporation of intercepted snow) and GLEAM (snow evaporation based on a modified Priestley–Taylor equation).
The CSNOW project provides Sentinel1 (S1) snow depth observations over the Northern Hemisphere mountains at a spatial resolution of ∼ 1 km (Lievens et al., 2019, 2022). For the European Alps the S1 snow depth dataset of Lievens et al. (2019) was improved (Lievens et al., 2022) over the period 2017–2019, and we use this dataset at ∼ 1 km resolution for the comparison with wflow_sbm S_{snow}. Because CSNOW provides snow depth observations and S_{snow} represents snow water equivalent, a direct comparison of spatial patterns is not feasible with the spatial efficiency metric E_{SP}. Therefore, we determine the temporal correlation for each grid cell (Fig. 15) and the spatiotemporal correlation between S1 snow depth retrievals and S_{snow}, estimated with the Spearman rankorder correlation coefficient r_{s}. The S1 snow depth retrievals are matched to the wflow_sbm grid using nearestneighbour resampling, excluding wet snow conditions detected by S1. The spatiotemporal correlation between S1 snow depth retrievals and S_{snow} is 0.88. Figure 15 shows generally a strong to a very strong positive correlation (median of r_{s}=0.87), except in the valleys. Lievens et al. (2022) indicated that the main uncertainties in S1 snow depth retrievals are mainly caused by wet, shallow and occasional snow cover and forest cover. Regional model simulations of snow depth to assess the S1 snow depth retrievals showed a strong positive temporal correlation if the maximum snow depth reached above ∼ 1 m at an elevation above ∼ 1000 m or with a forest cover fraction below ∼ 80 % (Lievens et al., 2022). A similar pattern is revealed by the temporal correlation between the S1 snow depth retrievals and wflow_sbm S_{snow} in Fig. 15.
4.2.6 Africa – Ouémé River – groundwater loss
The Ouémé mesoscale site (Benin), Africa, is part of the AMMACATCH observation network covering a 14 000 km^{2} basin in Sudanian climate on a crystalline basement. It is an interesting test case for wflow_sbm and the automated model setup including a priori estimation of the model parameters. Various studies using a variety of hydrological model concepts, including similar model concepts to wflow_sbm, have been conducted (see Cornelissen et al., 2013, and references therein) for this area. Séguis et al. (2011) reported that around 15 % of water is being lost to the groundwater which is disconnected from the river system. Here, we run the wflow_sbm model for the period 1981 to 2019 at a daily time step both without and with groundwater loss (L_{max}=0 mm d^{−1} vs. L_{max}=0.6 mm d^{−1}; 15 % of ∼ 4 mm d^{−1} of annual average daily rainfall) and analyse the model results for a variety of locations. L_{max} represents a maximum groundwater loss value, and the actual groundwater loss (computed by wflow_sbm) is controlled by the vertical hydraulic conductivity at the soil bottom and the saturated store (see Eq. 85) and may vary spatially and in time. Landcoverderived model parameters are based on VITO v2.0.2 (Buchhorn et al., 2019). For this model case, the snow model is disabled and lakes and reservoirs are not included. Figure 16 shows the KGE scores of discharge for the stations within the Ouémé mesoscale site without and with groundwater loss. Generally the performance of wflow_sbm increases with groundwater loss with a median increase of 0.86 for all stations. Figure 17 shows simulated discharge for the station TEBOU for 2010 with and without groundwater loss. The simulation with groundwater loss clearly shows less overestimation of discharge during the rising limb and peaks of the hydrograph, also reflected in the higher KGE score for the simulation with groundwater loss.
While the simulation with groundwater loss shows generally a better performance, we expect further improvement is possible by checking the effect of different groundwater loss values. For example, during the start of the rising limb of the hydrograph, the simulation underestimates the discharge for 2010 (Fig. 16) and other years (not shown). Applying a lower uniform groundwater loss value (we use the upper range of 15 % reported by Séguis et al., 2011) or applying spatially distributed groundwater loss values based on discharge measurements and a water balance approach for the upstream area (subcatchment) of each station could further improve simulation results.
We presented the wflow_sbm hydrological model as part of the Wflow.jl (v0.7.3) opensource modelling framework for distributed hydrological modelling in Julia, a continuation of the wflow development in the PCRaster Python framework. Wflow_sbm has been applied in various catchments around the world with satisfactory to good performance. With wflow_sbm we aim to strike a balance between lowresolution, lowcomplexity and highresolution, highcomplexity hydrological models. Most wflow_sbm parameters are based on physical characteristics, and at the same time wflow_sbm has a runtime performance well suited for largescale highresolution modelling. We demonstrated some examples of wflow_sbm applications with Wflow.jl (v0.7.3), using HydroMTWflow to set up the wflow_sbm model, including forcing, in an automated way. The wflow_sbm applications illustrate that the a priori model parameter estimates in combination with a manual adjustment of the f_{Kh0} model parameter result in generally satisfactory to good performance for discharge and catchmentaverage soil moisture (Moselle catchment), with, for the Umeälven catchment, similar performance for discharge (qualitatively) to EHYPE. For the Rhine catchment we demonstrated the ability of wflow_sbm to represent the spatial distribution of actual evapotranspiration E_{a} from GLEAM and snow storage S_{snow} for the Alps by comparing to CSNOW S1 snow depth observations. Wflow_sbm can represent daily E_{a} variability from GLEAM, although the spatial location matching is not satisfactory, which could be due to, amongst other things, the difference in resolution between GLEAM and wflow_sbm and different or missing process representations such as snow evaporation and irrigation. For the comparison of S_{snow} and S1 snow depth retrievals, a good performance indicated by the Spearman rankorder correlation coefficient (temporal (median of r_{s}=0.87) and spatiotemporal (r_{s}=0.88)) is obtained. The Ouémé River case illustrates the use of the model parameter L_{max} (maximum leakage) based on data from the literature, resulting in an overall significant increase in model performance. Including the process of leakage to deeper groundwater results in loss of water outside of the model domain, and we recommend including this process only if the scientific literature or geological data indicate that leakage to deeper groundwater is of importance. With the a priori parameter estimation, a working wflow_sbm model is set up quickly and incorrect process representations become apparent. The results, for example for the Ouémé River, indicate that local information and literature studies can help in improving process representation, and if they cannot, this opens the way for a better focus on the missing process representation. This is something that is lacking when a hydrological model is directly calibrated for a given catchment.
While (pedo)transfer functions are available for most of the sensitive wflow_sbm model parameters, this is not the case for the f_{Kh0} model parameter. An interesting approach, as part of future work focusing on model data, could be to develop a transfer function for this model parameter by estimating the transfer function with function space optimization (FSO), a method presented by Feigl et al. (2020). Relevant hydrological processes such as glacier and snow processes, evapotranspiration processes, unsaturated zone dynamics, (shallow) groundwater, and surface flow routing including lakes and reservoirs are part of wflow_sbm. Floodplain dynamics (backwater effects and floodplain storage) are not part of the kinematicwave routing in Wflow.jl, and this may be problematic for accurately simulating discharge and water depths when backwater effects and floodplain storage cannot be ignored (e.g. Zhao et al., 2017). Additionally, the kinematicwave approach is mostly applicable when slopes are steep and less reliable for lowgradient rivers. A recent wflow_sbm development, which is part of v0.7.3, is the improvement of the routing scheme for riverfloodplain dynamics and the improvement of discharge and water depth estimates for lowgradient rivers. The improved routing scheme includes the following options: (1) a 1D local inertial solution for river channel flow and a 2D local inertial solution for floodplain and overland flow, similar to Neal et al. (2012), and (2) a 1D local inertial solution for river channel flow with optional 1D floodplain schematization. Future possible developments related to the improved local inertial routing scheme are (1) to improve the multithreaded performance for the 2D local inertial solution, (2) vectorbased routing (e.g. Mizukami et al., 2021) allowing for more flexible channel routing configurations that are less computationally intensive and (3) to combine different routing solutions (e.g. kinematic wave and local inertial) at the submodeldomain scale (e.g. local inertial for the floodplain). For water resource modelling studies, wflow_sbm is often linked to a networkbased water allocation model (e.g. Meijer et al., 2021). The development of a water demand module (irrigation, livestock, industrial and domestic) and water allocation module is foreseen to fully exploit the gridded capabilities of wflow_sbm. The standard soil column of wflow_sbm extends to 2 m below surface level based on SoilGrids data, and although the soil column depth can be increased, the process modelled by wflow_sbm consists of shallow lateral subsurface flow, with an exponential decline in K_{v0} with depth, which may not be appropriate for simulating deep groundwater. While for many applications deep groundwater processes can be ignored, for the coupling with a groundwater model like MODFLOW or for the extraction of groundwater as part of the foreseen water demand and allocation developments, implementation of a deep groundwater concept is important. Finally, speedup of the wflow code is ongoing work: recently multithreading (single node) was added to the wflow code, and further developments may include distributed computing (using for example Julia's implementation Distributed.jl or the Julia interface to message passing interface (MPI) MPI.jl) and graphics processing unit (GPU) acceleration. In view of these future developments and the current status of the Wflow.jl framework, we have developed the wflow_sbm model, which is applicable worldwide and serves as an important tool to provide relevant information for operational and water resource planning challenges.
Note that variables are defined for model time step t. Some variables use the same Wflow.jl names to improve code readability, as they are handled in separate structs. ^{*}q_av
represents the average surface flow during model time step t.
Wflow.jl is opensource and distributed under the terms of the MIT License. The code and documentation are provided through the following GitHub repository: https://github.com/Deltares/Wflow.jl (last access: 2 April 2024). Wflow v0.7.3 is available through https://github.com/Deltares/Wflow.jl/releases/tag/v0.7.3 (last access: 2 April 2024), Zenodo, with the associated DOI https://doi.org/10.5281/zenodo.10495638 (van Verseveld et al., 2024), and is available as a Julia package. The wflow_sbm model cases presented in this paper are available at Zenodo with the associated DOI https://doi.org/10.5281/zenodo.10370017 (van Verseveld et al., 2023). Development and maintenance of Wflow.jl are conducted by Deltares, and we welcome contributions from external parties.
WJvV and MV wrote the majority of the source code, and JB, HB, DE, LB and MH contributed to the software development and documentation. WJvV wrote the original draft of the paper. AHW, ROI and WJvV performed the setup, analysed the model simulations and wrote Sect. 4.2. JB wrote Sect. 3 and produced Fig. 2. LB produced Fig. 1 (adapted from van Verseveld et al., 2024). All coauthors contributed to the review and editing of the original draft paper.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This work has received funding from the EU Horizon 2020 Programme for Research and Innovation, under grant agreement no 776613 (EUCP: European Climate Prediction system; https://www.eucpproject.eu, last access: 2 April 2024) and from the European High Performance Computing Joint Undertaking (JU) under grant agreement no. 955648 (ACROSS; https://www.acrossproject.eu/, last access: 2 April 2024). The research was cofunded by the Deltares Strategic Research Program “Natural Hazards and RealTime Information”.
This research has been supported by Horizon 2020 (EUCP (grant no. 776613), ACROSS (grant no. 955648)) and the Deltares Strategic Research Program “Natural Hazards and RealTime Information”.
This paper was edited by Mauro Cacace and reviewed by two anonymous referees.
Aerts, J. P. M., Hut, R. W., van de Giesen, N. C., Drost, N., van Verseveld, W. J., Weerts, A. H., and Hazenberg, P.: Largesample assessment of varying spatial resolution on the streamflow estimates of the wflow_sbm hydrological model, Hydrol. Earth Syst. Sci., 26, 4407–4430, https://doi.org/10.5194/hess2644072022, 2022. a
Alfieri, L., Burek, P., Dutra, E., Krzeminski, B., Muraro, D., Thielen, J., and Pappenberger, F.: GloFAS – global ensemble streamflow forecasting and flood early warning, Hydrol. Earth Syst. Sci., 17, 1161–1175, https://doi.org/10.5194/hess1711612013, 2013. a
Arino O., Ramos, J., Kalogirou, V., Defourny, P., and Achard., F.: GlobCover 2009, ESA Living Planet Symposium, 27 June–2 July 2010, Bergen, Norway, https://epic.awi.de/id/eprint/31046/1/Arino_et_al_GlobCover2009a.pdf (last access: 2 April 2024), 2010. a, b
Beck, H. E., Wood, E. F., Pan, M., Fisher, C. K., Miralles, D. M., van Dijk, A. I. J. M., McVicar, T. R., and Adler, R. F.: MSWEP V2 global 3hourly 0.1° precipitation: Methodology and quantitative assessment, B. Am. Meteorol. Soc., 100, 473–500, 2019. a
Bell, V. A., Kay, A. L., Jones, R. G., and Moore, R. J.: Development of a high resolution gridbased river flow model for use with regional climate model output, Hydrol. Earth Syst. Sci., 11, 532–549, https://doi.org/10.5194/hess115322007, 2007. a, b
Benning, R. G.: Towards a new lumped parameterization at catchment scale, PhD thesis, Wageningen University, Wageningen, The Netherlands, https://edepot.wur.nl/216531 (last access: 2 April 2024), 1995. a, b
Bergström, S.: The HBV model – its structure and applications, SMHI Reports Hydrology, Norrköping, Sweden, RH 4, 1992. a, b, c
Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B.: Julia: A Fresh Approach to Numerical Computing, SIAM Rev., 59, 65–98, https://doi.org/10.1137/141000671, 2017. a
Beven, K.: Prophecy, reality, and uncertainty in distributed hydrological modeling, Adv. Water Resour., 16, 41–51, 1993. a
Beven, K.: A manifesto for the equifinality thesis, J. Hydrol., 320, 18–36, 2006. a
Bierkens M. F. P., Bell V. A., Burek, P., Chaney, N., Condon, L. E., David, C. H., de Roo, A., Döll P., Drost, N., Famiglietti, J. S., Flörke, M., Gochis, D. J., Houser, P., Hut, R., Keune, J., Kollet, S., Maxwell, R. M., Reager, J. T., Samaniego, L., Sudicky, E., Sutanudjaja, E. H., van de Giesen, N., Winsemius, H., and Wood, E. F.: Hyperresolution global hydrological modelling: What is next? “Everywhere and locally relevant”, Hydrol. Process., 29, 310–320, 2014. a
Brakensiek, D. L., Rawls, W. J., and Stephenson, G. R.: Modifying scs hydrologic soil groups and curve numbers for rangeland soils, ASAE Paper no. PNR84203, St. Joseph, Michigan, USA, 1984. a
Brooks, R. H. and Corey, A. T.: Hydraulic properties of porous media, Hydrology Papers 3, Colorado State University, Fort Collins, 27 pp., 1964. a, b, c
Brunner, P. and Simmons, C. T.: HydroGeoSphere: A Fully Integrated, Physically Based Hydrological Model, Groundwater, 50, 170–176, 2012. a
Buchhorn, M., Smets, B., Bertels, L., Lesiv, M., Tsendbazar, N.E., Herold, M., and Fritz, S.: Copernicus Global Land Service: Land Cover 100m, epoch 2015, Globe (Version V2.0.2), Zenodo [data set], https://doi.org/10.5281/zenodo.3243509, 2019. a, b, c, d
Chow, V. T., Maidment, D. R., and Mays, L. W.: Applied Hydrology, McGrawHill, New York, 572 pp., ISBN 0070108102, 1988. a, b
Clark, M. P., Schaefli, B., Schymanski, S. J., Samaniego, L., Luce, C. H., Jackson, B. M., Freer, J. E., Arnold, J. R., Moore, R. D., Istanbulluoglu, E., and Ceola, S.: Improving the theoretical underpinnings of processbased hydrologic models, Water Resour. Res., 52, 2350–2365, https://doi.org/10.1002/2015WR017910, 2016. a
Clark, M. P., Bierkens, M. F. P., Samaniego, L., Woods, R. A., Uijlenhoet, R., Bennett, K. E., Pauwels, V. R. N., Cai, X., Wood, A. W., and PetersLidard, C. D.: The evolution of processbased hydrologic models: historical challenges and the collective quest for physical realism, Hydrol. Earth Syst. Sci., 21, 3427–3440, https://doi.org/10.5194/hess2134272017, 2017. a, b
Cornelissen, T., Diekkrüger, B., and Giertz, S.: A comparison of hydrological models for assessing the impact of land use and climate change on discharge in a tropical catchment, J. Hydrol., 498, 221–236, https://doi.org/10.1016/j.jhydrol.2013.06.016, 2013. a
de Bruin, H. A. R., Trigo, I. F., Bosveld, F. C., and Meirink, J. F.: A Thermodynamically Based Model for Actual Evapotranspiration of an Extensive Grass Field Close to FAO Reference, Suitable for Remote Sensing Application, J. Hydrometeorol., 17, 1373–1382, https://doi.org/10.1175/JHMD150006.1, 2016. a
Dembélé, M., Hrachowitz, M., Savenije, H. H. G., Mariéthoz, G., and Schaefli, B.: Improving the predictive skill of a distributed hydrological model by calibration on spatial patterns with multiple satellite data sets, Water Resour. Res., 56, e2019WR026085, https://doi.org/10.1029/2019WR026085, 2020. a, b
Demirel, M. C., Mai, J., Mendiguren, G., Koch, J., Samaniego, L., and Stisen, S.: Combining satellite data and appropriate objective functions for improved spatial pattern performance of a distributed hydrologic model, Hydrol. Earth Syst. Sci., 22, 1299–1315, https://doi.org/10.5194/hess2212992018, 2018. a
Eilander, D. and Boisgontier, H.: hydroMT (v0.4.5), Zenodo [code], https://doi.org/10.5281/zenodo.6107669, 2022. a
Eilander, D., van Verseveld, W., Yamazaki, D., Weerts, A., Winsemius, H. C., and Ward, P. J.: A hydrography upscaling method for scaleinvariant parametrization of distributed hydrological models, Hydrol. Earth Syst. Sci., 25, 5287–5313, https://doi.org/10.5194/hess2552872021, 2021. a, b, c, d
Eilander, D., Boisgontier, H., van Verseveld, W., Bouaziz, L., and Hegnauer, M.: hydroMTwflow (v0.1.4), Zenodo [code], https://doi.org/10.5281/zenodo.6221375, 2022. a, b, c, d, e
Engman, E.: Roughness coeffcients for routing surface runoff, J. Irrig. Drain. Eng., 112, 39–53, https://doi.org/10.1061/(ASCE)07339437(1986)112:1(39), 1986. a
ESDAC: The European Soil Database distribution version 2.0, European Commission and the European Soil Bureau Network, https://esdac.jrc.ec.europa.eu (last access: 2 April 2024), 2004. a
European Environment Agency: Corine Land Cover (CLC) 2018, Version 2020_20u1, Copernicus Services [data set], https://land.copernicus.eu/paneuropean/corinelandcover/clc2018 (last access: 2 April 2024), 2018. a, b, c, d
Fan, J., McConkey, B., Wang, H., and Janzen, H.: Root distribution by depth for temperate agricultural crops, Field Crops Res., 189, 68–74, https://doi.org/10.1016/j.fcr.2016.02.013, 2016. a
Feddes, R. A., Kowalik, P. J., and Zaradny, H.: Simulation of field water use and crop yield, Simulation Monographs, Pudoc, Wageningen, 189 pp., 1978. a, b, c, d
Feigl, M., Herrnegger, M., Klotz, D., and Schulz, K.: Function space optimization: A symbolic regression method for estimating parameter transfer functions for hydrological models, Water Resour. Res., 56, e2020WR027385, https://doi.org/10.1029/2020WR027385, 2020. a
Fenicia, F., Kavetski, D., and Savenije, H. H.: Elements of a ﬂexible approach for conceptual hydrological modeling: 1. Motivation and theoretical development, Water Resour. Res., 47, W11510, https://doi.org/10.1029/2010WR010174, 2011. a
Funk, C., Peterson, P., Landsfeld, M., Pedreros, D., Verdin, J., Shukla, S., Husak, G., Rowland, J., Harrison, L., Hoell, A., and Michaelsen, J.: The climate hazards infrared precipitation with stations – a new environmental record for monitoring extremes, Sci. Data 2, 150066, https://doi.org/10.1038/sdata.2015.66, 2015. a
Gao, H., Hrachowitz, M., Fenicia, F., Gharari, S., and Savenije, H. H. G.: Testing the realism of a topographydriven model (FLEXTopo) in the nested catchments of the Upper Heihe, China, Hydrol. Earth Syst. Sci., 18, 1895–1915, https://doi.org/10.5194/hess1818952014, 2014. a
Gash, J. H. C.: An analytical model of rainfall interception by forests, Q. J. Roy. Meteor. Soc., 105, 43–55, https://doi.org/10.1002/qj.49710544304, 1979. a, b
Gebremicael, T., Mohamed, Y., and der Zaag, P. V.: Attributing the hydrological impact of different land use types and their longterm dynamics through combining parsimonious hydrological modelling, alteration analysis and PLSR analysis, Sci. Total Environ., 660, 1155–1167, https://doi.org/10.1016/j.scitotenv.2019.01.085, 2019. a
Giardino, A., Schrijvershof, R., Nederhoff, C., De Vroeg, H., Brière, C., Tonnon, P. K., Caires, S., Walstra, D., Sosa, J., Van Verseveld, W., Schellekens, J., and Sloff, C. J.: A quantitative assessment of human interventions and climate change on the West African sediment budget, Ocean Coast. Manag., 156, 249–265, https://doi.org/10.1016/j.ocecoaman.2017.11.008, 2018. a
Gupta, H., Kling, H., Yilmaz, K., and Martinez, G.: Decomposition of the mean squared error and nse performance criteria: Implications for improving hydrological modelling, J. Hydrol., 377, 80–91, https://doi.org/10.1016/j.jhydrol.2009.08.003, 2009. a, b
Hassaballah, K., Mohamed, Y., Uhlenbrook, S., and Biro, K.: Analysis of streamflow response to land use and land cover changes using satellite data and hydrological modelling: case study of Dinder and Rahad tributaries of the Blue Nile (Ethiopia–Sudan), Hydrol. Earth Syst. Sci., 21, 5217–5242, https://doi.org/10.5194/hess2152172017, 2017. a
Hengl, T., Mendes de Jesus, J., Heuvelink, G. B. M., Ruiperez Gonzalez, M., Kilibarda, M., Blagotić, A., Shangguan, W., Wright, M. N., Geng, X., BauerMarschallinger, B., Guevara, M. A., Vargas, R., MacMillan, R. A., Batjes, N. H., Leenaars, J. G. B., Ribeiro, E., Wheeler, I., Mantel, S., and Kempen, B.: SoilGrids250m: Global gridded soil information based on machine learning, PLoS ONE, 12, e0169748, https://doi.org/10.1371/journal.pone.0169748, 2017. a, b, c, d
Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., MuñozSabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flem ming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.N.: The ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999–2049, https://doi.org/10.1002/qj.3803, 2020. a
Horn, B. K. P.: Hill Shading and the Reflectance Map, Proceedings of IEEE, 69, 14–47, https://doi.org/10.1109/PROC.1981.11918, 1981. a
Hrachowitz, M. and Clark, M. P.: HESS Opinions: The complementary merits of competing modelling philosophies in hydrology, Hydrol. Earth Syst. Sci., 21, 3953–3973, https://doi.org/10.5194/hess2139532017, 2017. a
Imhoff, R. O., van Verseveld, W. J., Osnabrugge, B., and Weerts, A. H.: Scaling PointScale (Pedo)transfer Functions to Seamless LargeDomain Parameter Estimates for HighResolution Distributed Hydrologic Modeling: An Example for the Rhine River, Water Resour. Res., 56, e2019WR026807, https://doi.org/10.1029/2019WR026807, 2020. a, b, c, d, e, f, g, h, i
Jägermeyr, J., Pastor, A., Biemans, H., and Gerten, D.:Reconciling irrigated food production with environmental flows for Sustainable Development Goals implementation, Nat. Commun., 8, 15900, https://doi.org/10.1038/ncomms15900, 2017. a
Jakeman, A. J. and Hornberger, G. M.: How Much Complexity Is Warranted in a RainfallRunoff Model?, Water Resour. Res., 29, 2637–2649, 1993. a
Karssenberg, D., Schmitz, O., Salamon, P., de Jong, K., and Bierkens, M. F. P.: A software framework for construction of processbased stochastic spatiotemporal models and data assimilation, Environ. Modell. Softw., 25, 489–502, https://doi.org/10.1016/j.envsoft.2009.10.004, 2010. a
Kilgore, J. L.: Development and evaluation of a gisbased spatially distributed unit hydrograph model, MSc. thesis, Virginia Tech, http://hdl.handle.net/10919/35777 (last access: 2 April 2024), 1997. a
Kirchner, J. W.: Getting the right answers for the right reasons: Linking measurements, analyses, and models to advance the science of hydrology, Water Resour. Res., 42, W03S04, https://doi.org/10.1029/2005WR004362, 2006. a
Kling, H., Fuchs, M., and Paulin, M.: Runoff conditions in the upper Danube basin under an ensemble of climate change scenarios, J. Hydrol., 424425, 264–277, https://doi.org/10.1016/j.jhydrol.2012.01.011, 2012. a, b
Knoben, W. J. M., Freer, J. E., and Woods, R. A.: Technical note: Inherent benchmark or not? Comparing Nash–Sutcliffe and Kling–Gupta efficiency scores, Hydrol. Earth Syst. Sci., 23, 4323–4331, https://doi.org/10.5194/hess2343232019, 2019. a
Knoben, W. J. M., Clark, M. P., Bales, J., Bennett, A., Gharari, S., Marsh, C. B., Nijssen, B., Pietroniro, A., Spiteri, R. J., Tarboton, D. G., and Wood, A. W.: Community Workflows to Advance Reproducibility in Hydrologic Modeling: Separating modelagnostic and modelspecific configuration steps in applications of largedomain hydrologic models, Earth and Space Science Open Archive, https://doi.org/10.1002/essoar.10509195.1, 2021. a
Koch, J., Demirel, M. C., and Stisen, S.: The SPAtial EFficiency metric (SPAEF): multiplecomponent evaluation of spatial patterns for optimization of hydrological models, Geosci. Model Dev., 11, 1873–1886, https://doi.org/10.5194/gmd1118732018, 2018. a, b
Kollet, S. J. and Maxwell, R. M.: Integrated surface groundwater flow modeling: A free surface overland flow boundary condition in a parallel groundwater flow model, Adv. Water Resour., 29, 945–958, https://doi.org/10.1016/j.advwatres.2005.08.006, 2006. a
LaverdeBarajas, M., Corzo Perez, G. A., Chishtie, F., Poortinga, A., Uijlenhoet, R., and Solomatine, D. P.: Decomposing satellitebased rainfall errors in flood estimation: Hydrological responses using a spatiotemporal objectbased verification method, J. Hydrol., 591, 125554, https://doi.org/10.1016/j.jhydrol.2020.125554, 2020. a
Lehner, B., Reidy Liermann, C., Revenga, C., Vörösmarty, C., Fekete, B., Crouzet, P., Döll, P., Endejan, M., Frenken, K., Magome, J., Nilsson C., Robertson, J. C., Rodel, R., Sindorf, N., and Wisser, D.: HighResolution Mapping of the World's Reservoirs and Dams for Sustainable River Flow Management, Front. Ecol. Environ., 9, 494–502, https://doi.org/10.1890/100125, 2011. a, b
Lievens, H., Demuzere, M., Marshall, H.P., Reichle, R. H., Brucker, L., Brangers, I., de Rosnay, P., Dumont, M., Girotto, M., Immerzeel, W. W., Jonas, T., Kim, E. J., Koch, I., Marty, C., Saloranta, T., Schöber J., and De Lannoy, G. J. M.: Snow depth variability in the Northern Hemisphere mountains observed from space, Nat. Commun., 10, 4629, https://doi.org/10.1038/s4146701912566y, 2019. a, b
Lievens, H., Brangers, I., Marshall, H.P., Jonas, T., Olefs, M., and De Lannoy, G.: Sentinel1 snow depth retrieval at subkilometer resolution over the European Alps, The Cryosphere, 16, 159–177, https://doi.org/10.5194/tc161592022, 2022. a, b, c, d
Lin, P., Pan, M., Allen, G., Frasson, R., Zeng, Z., Yamazaki, D., and Wood, E.: Global estimates of reachlevel bankfull river width leveraging bigdata geospatial analysis, Zenodo [data set], https://doi.org/10.5281/zenodo.3552776, 2019. a, b, c, d
Liu, S.: Estimation of rainfall storage capacity in the canopies of cypress wetlands and slash pine uplands in NorthCentral Florida, J. Hydrol., 207, 32–41, https://doi.org/10.1016/S00221694(98)001152, 1998. a, b, c
Liu, Z., Martina, M. L. V., and Todini, E.: Flood forecasting using a fully distributed model: application of the TOPKAPI model to the Upper Xixian Catchment, Hydrol. Earth Syst. Sci., 9, 347–364, https://doi.org/10.5194/hess93472005, 2005. a
López López, P., Wanders, N., Schellekens, J., Renzullo, L. J., Sutanudjaja, E. H., and Bierkens, M. F. P.: Improved largescale hydrological modelling through the assimilation of streamflow and downscaled satellite soil moisture observations, Hydrol. Earth Syst. Sci., 20, 3059–3076, https://doi.org/10.5194/hess2030592016, 2016. a, b
Martens, B., Miralles, D. G., Lievens, H., van der Schalie, R., de Jeu, R. A. M., FernándezPrieto, D., Beck, H. E., Dorigo, W. A., and Verhoest, N. E. C.: GLEAM v3: satellitebased land evaporation and rootzone soil moisture, Geosci. Model Dev., 10, 1903–1925, https://doi.org/10.5194/gmd1019032017, 2017. a
Maurer, E., Wood, A., Adam, J., Lettenmaier, D., and Nijssen, B.: A longterm hydrologically based dataset of land surface fluxes and states for the conterminous United States, J. Climate, 15, 3237–3251, https://doi.org/10.1175/15200442(2002)015<3237:ALTHBD>2.0.CO;2, 2002. a
Meijer, K., Verschelling, E., van Verseveld, W., Donchyts, G., Schmeier, S., and Kwadijk, J.: Fit for purpose? Rapid development of water allocation models using global data: Application for the Upper Niger Basin, Environ. Modell. Softw., 145, 105168, https://doi.org/10.1016/j.envsoft.2021.105168, 2021. a, b, c
Melsen, L. A., Teuling, A. J., Torfs, P. J. J. F., Uijlenhoet, R., Mizukami, N., and Clark, M. P.: HESS Opinions: The need for processbased evaluation of largedomain hyperresolution models, Hydrol. Earth Syst. Sci., 20, 1069–1079, https://doi.org/10.5194/hess2010692016, 2016. a
Messager, M. L., Lehner, B., Grill, G., Nedeva, I., and Schmitt, O.: Estimating the volume and age of water stored in global lakes using a geostatistical approach, Nat. Commun., 7, 13603, https://doi.org/10.1038/ncomms13603, 2016. a, b, c, d
Miralles, D. G., Holmes, T. R. H., De Jeu, R. A. M., Gash, J. H., Meesters, A. G. C. A., and Dolman, A. J.: Global landsurface evaporation estimated from satellitebased observations, Hydrol. Earth Syst. Sci., 15, 453–469, https://doi.org/10.5194/hess154532011, 2011. a, b
Mizukami, N., Clark, M. P., Gharari, S., Kluzek, E., Pan, M., Lin, P., Beck, H. E., and Yamazaki, D.: A vectorbased river routing model for Earth System Models: Parallelization and global applications, J. Adv. Model. Earth Sy., 13, e2020MS002434, https://doi.org/10.1029/2020MS002434, 2021. a
Myneni, R. B., Knyazikhin, Y., and Park, T.: MCD15A3H MODIS/Terra+Aqua Leaf Area Index/FPAR 4day L4 Global 500m SIN Grid V006, NASA EOSDIS Land Processes Distributed Active Archive Center [data set], https://doi.org/10.5067/MODIS/MCD15A3H.006, 2015. a, b
Neal, J., Schumann, G., and Bates P.: A subgrid channel model for simulating river hydraulics and floodplain inundation over large and data sparse areas, Water Resour. Res., 48, W11506, https://doi.org/10.1029/2012WR012514, 2012. a
O'Neill, P. E., Chan, S., Njoku, E. G., Jackson, T., Bindlish, R., Chaubell, J., and Colliander, A.: SMAP Enhanced L3 Radiometer Global and Polar Grid Daily 9 km EASEGrid Soil Moisture, Version 5, NASA National Snow and Ice Data Center Distributed Active Archive Center, Boulder, Colorado USA, [data set], https://doi.org/10.5067/4DQ54OUIJ9DL, 2021. a
Pekel, J. F., Cottam, A., Gorelick, N., and Belward, A. S.: High resolution mapping of global surface water and its long term changes, Nature, 540, 418–422, https://doi.org/10.1038/nature20584, 2016. a, b
Pfeffer, W. T., Arendt, A. A., Bliss, A., Bolch, T., Cogley, J. G., Gardner, A. S., Hagen, J.O., Hock, R., Kaser, G., Kienholz, C., Miles, E. S., Moholdt, G., Mölg, N., Paul, F., Radic, V., Rastner, P., Raup, B. H., Rich, J., Sharp, M. J., and Consortium, T. R.: The Randolph Glacier Inventory: a globally complete inventory of glaciers, J. Glaciol., 60, 537–552, https://doi.org/10.3189/2014JoG13J176, 2014. a, b
Pitman, J.: Rainfall interception by bracken in open habitats–relations between leaf area, canopy storage and drainage rate, J. Hydrol., 105, 317–334, https://doi.org/10.1016/00221694(89)90111X, 1989. a, b, c
Rawls, W. J. and Brakensiek, D. L.: Estimation of Soil Water Retention and Hydraulic Properties, in: Unsaturated ﬂow in hydrologic modelling – Theory and practice, edited by: MorelSeytoux, H. J., NATO ASI Series 9, 275–300, Dordrecht, Kluwer Academic Publishing, the Netherlands, ISBN 9789401075596, 1989. a
Rusli, S. R., Weerts, A. H., Taufiq, A., and Bense, V. F.: Estimating water balance components and their uncertainty bounds in highly groundwaterdependent and datascarce area: An example for the Upper Citarum basin, Journal of Hydrology: Regional Studies, 37, 100911, https://doi.org/10.1016/j.ejrh.2021.100911, 2021. a, b
Rutter, A., Kershaw, K., Robins, P., and Morton, A.: A predictive model of rainfall interception in forests, 1. Derivation of the model from observations in a plantation of Corsican pine Agric. Meteorol., 9, 367–384, https://doi.org/10.1016/00021571(71)900343, 1971. a, b
Rutter, A., Morton, A., and Robins, P.: A predictive model of rainfall interception in forests. II. Generalization of the model and comparison with observations in some coniferous and hardwood stands J. Appl. Ecol., 12, 367–380, https://doi.org/10.2307/2401739, 1975. a, b
Samaniego, L., Kumar, R., and Attinger, S.: Multiscale parameter regionalization of a gridbased hydrologic model at the mesoscale, Water Resour. Res., 46, W05523, https://doi.org/10.1029/2008WR007327, 2010. a
Samaniego, L., Kumar, R., Thober, S., Rakovec, O., Zink, M., Wanders, N., Eisner, S., Müller Schmied, H., Sutanudjaja, E. H., WarrachSagi, K., and Attinger, S.: Toward seamless hydrologic predictions across spatial scales, Hydrol. Earth Syst. Sci., 21, 4323–4346, https://doi.org/10.5194/hess2143232017, 2017. a
Schellekens, J., van Verseveld, W., Visser, M., Winsemius, H., Bouaziz, L., Euser, T., De Vries, S., Thiange, C., Boisgontier, H., Eilander, D., Tollenaar, D., Weerts, A., Baart, F., Hazenberg, P., Pronk, M., Lutz, A., Ten Velden, C., Benedict, I., and Jansen, M.: openstreams/wflow: Bug fixes and updates for release 2020.1.2, Zenodo [code], https://doi.org/10.5281/zenodo.4291730, 2020. a, b, c
Schenk, H. J. and Jackson, R. B.: The global biogeography of roots. Ecological monographs, 72, 311–328, https://doi.org/10.1890/00129615(2002)072[0311:TGBOR]2.0.CO;2, 2002. a
Seibert, J., Vis, M. J. P., Kohn, I., Weiler, M., and Stahl, K.: Technical note: Representing glacier geometry changes in a semidistributed hydrological model, Hydrol. Earth Syst. Sci., 22, 2211–2224, https://doi.org/10.5194/hess2222112018, 2018. a, b
Séguis, L., Kamagaté, B., Favreau, G., Descloitres, M., Seidel, J.L., Galle, S., Peugeot, C., Gosset, M., Le Barbé, L., Malinur, F., Van Exter, S., Arjounin, M., Boubkraoui, S., and Wubda, M.: Origins of streamflow in a crystalline basement catchment in a subhumid Sudanian zone: The Donga basin (Benin, West Africa): Interannual variability of water budget, J. Hydrol., 402, 1–13, https://doi.org/10.1016/j.jhydrol.2011.01.054, 2011. a, b
Ŝimůnek, J., van Genuchten, M. T., and Ŝejna, M.: Development and applications of the HYDRUS and STANMOD software packages and related codes, Vadose Zone J., 7, 587–600, 2008. a
Sperna Weiland, F. C., Visser, R. D., Greve, P., Bisselink, B., Brunner, L., and Weerts, A. H.: Estimating Regionalized Hydrological Impacts of Climate Change Over Europe by PerformanceBased Weighting of CORDEX Projections, Frontiers of Water, 3, 713537, https://doi.org/10.3389/frwa.2021.713537, 2021. a, b
Swedish Meteorological and Hydrological Institute: SMHI Vattenweb, https://www.smhi.se/data/hydrologi/vattenwebb, last access: 3 April 2024. a
Swedish Meteorological and Hydrological Institute: Model Performance Europe, https://hypeweb.smhi.se/explorewater/modelperformances/modelperformanceeurope/, last access: 2 June 2022. a
Tanaka, T. and Tachikawa, Y.: Testing the applicability of a kinematic wavebased distributed hydrologic model in two climatically contrasting catchments, Hydrolog. Sci. J., 60, 1361–1373, https://doi.org/10.1080/02626667.2014.967693, 2015. a, b
Todini, E. and Ciarapica, L.: Mathematical models of large watershed hydrology, edited by: Singh, V. P. and Frevert, D. K., Water Resources Publications, Littleton, Colorado, 471–506, ISBN 9781887201346, 2002. a, b
Tóth, B., Weynants, M., Nemes, A., Makó, A., Bilas, G., and Tóth, G.: New generation of hydraulic pedotransfer functions for Europe, Eer. J. Soil Sci., 66, 226–238, https://doi.org/10.1111/ejss.12192, 2015. a, b
Trambauer, P., Werner, M., Winsemius, H. C., Maskey, S., Dutra, E., and Uhlenbrook, S.: Hydrological drought forecasting and skill assessment for the Limpopo River basin, southern Africa, Hydrol. Earth Syst. Sci., 19, 1695–1711, https://doi.org/10.5194/hess1916952015, 2015. a
van Beek, L. P. H., Wada, Y., and Bierkens, M. F. P.: Global monthly water stress: 1. Water balance and water availability, Water Resour. Res., 47, W07517, https://doi.org/10.1029/2010WR009791, 2011. a
van Dijk, A. I. J. M. and Bruijnzeel, L. A.: Modelling rainfall interception by vegetation of variable density using an adapted analytical model, Part 2, Model validation for a tropical upland mixed cropping system, J. Hydr., 247, 239–262, 2001. a, b, c
van Looy, K., Bouma, J., Herbst, M., Koestel, J., Minasny, B., Mishra, U., and Vereecken, H.: Pedotransfer functions in Earth system science: Challenges and perspectives, Rev. Geophys., 55, 1199–1256, https://doi.org/10.1002/2017RG000581, 2017. a
van Verseveld, W., Weerts, A., and Imhoff, R.: Wflow_sbm v0.7.3 model cases (0.1), Zenodo [data set], https://doi.org/10.5281/zenodo.10370017, 2023. a
van Verseveld, W., Visser, M., Buitink, J., Bouaziz, L., Boisgontier, H., Bootsma, H., Pronk, M., Hartgring, S., Eilander, D., Weerts, A., Dalmijn, B., Hofer, J., and Hegnauer, M.: Wflow.jl (v0.7.3), Zenodo [code], https://doi.org/10.5281/zenodo.10495638, 2024. a, b, c, d, e, f, g, h
Vertessy, R. and Elsenbeer, H.: Distributed modeling of storm ﬂow generation in an amazonian rain forest catchment: effects of model parameterization, Water Resour. Res., 35, 2173–2187, https://doi.org/10.1029/1999WR900051, 1999. a, b, c, d, e
Wang, X., Huo, Z., Feng, S., Guo, P., and Guan, H.: Estimating groundwater evapotranspiration from irrigated cropland incorporating root zone soil texture and moisture dynamics, J. Hydrol., 543, 501–509, https://doi.org/10.1016/j.jhydrol.2016.10.027, 2016. a
Wannasin, C., Brauer, C., Uijlenhoet, R., van Verseveld, W., and Weerts, A.: Daily flow simulation in thailand part I: testing a distributed hydrological model with seamless parameter maps based on global data, Journal of Hydrology: Regional Studies, 34, 100794, https://doi.org/10.1016/j.ejrh.2021.100794, 2021a. a, b, c
Wannasin, C., Brauer, C., Uijlenhoet, R., van Verseveld, W., and Weerts, A.: Daily flow simulation in thailand part II: unraveling effects of reservoir operation, Journal of Hydrology: Regional Studies, 34, 100792, https://doi.org/10.1016/j.ejrh.2021.100792, 2021b. a, b
Wigmosta, M. S., Lane, L. J., Tagestad, J. D., and Coleman A. M.: Hydrologic and erosion models to assess land use and management practices affecting soil erosion, J. Hydrol. Eng., 14, 27–41, 2009. a
Wood, E. F., Roundy, J. K., Troy, T. J., van Beek, L. P. H., Bierkens, M. F. P., Blyth, E., de Roo, A., Döll, P., Ek, M., Famiglietti, J., Gochis, D., van de Giesen, N., Houser, P., Jaffé, P. R., Kollet, S., Lehner, B., Lettenmaier, D. P., PetersLidard, C., Sivapalan, M., Sheffield, J., Wade, A., and Whitehead, P.: Hyperresolution global land surface modeling: Meeting a grand challenge for monitoring Earth’s terrestrial water, Water Resour. Res., 47, W05301, https://doi.org/10.1029/2010WR010090, 2011. a
Yamazaki, D., Trigg, M. A., and Ikeshima, D.: Development of a global 90 m water body map using multitemporal Landsat images, Remote Sens. Environ., 171, 337–351, https://doi.org/10.1016/j.rse.2015.10.014, 2015. a
Yamazaki, D., Ikeshima, D., Tawatari, R., Yamaguchi, T., O'Loughlin, F., Neal, J. C., Sampson, C. C., Kanae, S., and Bates, P. D.: A high accuracy map of global terrain elevations, Geophys. Res. Lett., 44, 5844–5853, https://doi.org/10.1002/2017GL072874, 2017. a
Yamazaki, D., Ikeshima, D., Sosa, J., Bates, P. D., Allen, G. H., and Pavelsky, T. M.: MERIT Hydro: a high resolution global hydrography map based on latest topography dataset, Water Resour. Res., 55, 5053–5073, https://doi.org/10.1029/2019WR024873, 2019. a, b
Yang, F., Zhang, G., Yin, X., Liu, Z., and Huang, Z.: Study on capillary rise from shallow groundwater and critical water table depth of a salinesodic soil in western Songnen plain of China, Environ. Earth Sci. 64, 2119–2126, https://doi.org/10.1007/s1266501110384, 2011. a
Zammouri, M.: Case study of water table evaporation at Ichkeul Marshes (Tunisia), J. Irrig. Drain. Eng. 127, 265–271, 2001. a
Zhao, F., Veldkamp, T. I. E., Frieler, K., Schewe, J., Ostberg, S., Willner, S., Schauberger, B., Gosling, S. N., Schmied, H. M., Portmann, F. T., Leng, G., Huang, M., Liu, X., Tang, Q., Hanasaki, N., Biemans, H., Gerten, D., Satoh, Y., Pokhrel, Y., and Yamazaki, D.: The critical role of the routing scheme in simulating peak river discharge in global hydrological models, Environ. Res. Lett., 12, 075003, https://doi.org/10.1088/17489326/aa7250, 2017. a