Behind the scenes of streamflow model performance

Abstract

Streamflow is often the only variable used to evaluate hydrological models. In a previous international comparison study, eight research groups followed an identical protocol to calibrate 12 hydrological models using observed streamflow of catchments within the Meuse basin. In the current study, we quantify the differences in five states and fluxes of these 12 process-based models with similar streamflow performance, in a systematic and comprehensive way. Next, we assess model behavior plausibility by ranking the models for a set of criteria using streamflow and remote-sensing data of evaporation, snow cover, soil moisture and total storage anomalies. We found substantial dissimilarities between models for annual interception and seasonal evaporation rates, the annual number of days with water stored as snow, the mean annual maximum snow storage and the size of the root-zone storage capacity. These differences in internal process representation imply that these models cannot all simultaneously be close to reality. Modeled annual evaporation rates are consistent with Global Land Evaporation Amsterdam Model (GLEAM) estimates. However, there is a large uncertainty in modeled and remote-sensing annual interception. Substantial differences are also found between Moderate Resolution Imaging Spectroradiometer (MODIS) and modeled number of days with snow storage. Models with relatively small root-zone storage capacities and without root water uptake reduction under dry conditions tend to have an empty root-zone storage for several days each summer, while this is not suggested by remote-sensing data of evaporation, soil moisture and vegetation indices. On the other hand, models with relatively large root-zone storage capacities tend to overestimate very dry total storage anomalies of the Gravity Recovery and Climate Experiment (GRACE). None of the models is systematically consistent with the information available from all different (remote-sensing) data sources. Yet we did not reject models given the uncertainties in these data sources and their changing relevance for the system under investigation.

1 Introduction

Hydrological models are valuable tools for short-term forecasting of river flows and long-term predictions for strategic water management planning, but also to develop a better understanding of the complex interactions of water storage and release processes at the catchment scale. In spite of the wide variety of existing hydrological models, they mostly include similar functionalities of storage, transmission and release of water to represent the dominant hydrological processes of a particular river basin , differing mostly only in the detail of their parameterizations .

In all of these models, each individual model component constitutes a separate hypothesis of how water moves through that specific part of the system. Frequently, the individual hypotheses remain untested. Instead, only the model output, i.e., the aggregated response of these multiple hypotheses, is confronted with data of the aggregated response of a catchment to atmospheric forcing. Countless applications of different hydrological models in many different regions across the world over the last decades have shown that these models often provide relatively robust estimates of streamflow dynamics, for both calibration and evaluation periods. However, various combinations of different untested individual hypotheses can and do lead to similar aggregated outputs, i.e., model equifinality .

To be useful for any of the above applications, it is thus of critical importance that not only the aggregated but also the individual behaviors of the respective hypotheses are consistent with their real-world equivalents. Given the complexity and heterogeneity of natural systems together with the general lack of suitable observations, this remains a major challenge in hydrology e.g.,.

Studies have addressed the issue by constraining the parameters of specific models through the use of additional data sources besides streamflow. , and mapped saturated contributing areas during field surveys to constrain model parameters, while patterns of water tables in piezometers were used by , and . Other sources include satellite-based total water storage anomalies e.g.,, evaporation e.g.,, near-surface soil moisture e.g.,, snow cover information e.g.,, or a combination of these variables e.g.,. Reflecting the results of many studies, showed that streamflow is necessary but not sufficient to constrain model components to warrant partitioning of incoming precipitation to storage, evaporation and drainage.

Hydrological simulations are, however, not only affected by model parameter uncertainty, but also by the selection of a model structure and its parameterization (i.e., the choice of equations). Modeling efforts over the last 4 decades have led to a wide variety of hydrological models providing flexibility to test competing modeling philosophies, from spatially lumped model representations of the system to high-resolution small-scale processes numerically integrated to the catchment scale . and compared global hydrological models and found that differences between models are a major source of uncertainty. Nonetheless, model selection is often driven by personal preference and experience of individual modelers rather than detailed model test procedures .

A suite of comparison experiments tested and explored differences between alternative modeling structures and parameterizations . However, these studies mostly restricted themselves to analyses of the models' skills in reproducing streamflow (“aggregated hypothesis”), with little consideration for the model-internal processes (“individual hypotheses”). The Framework for Understanding Structural Errors (FUSE) was one of the first initiatives towards a more comprehensive assessment of model structural errors, with special consideration given to individual hypotheses .

Subsequent efforts towards more rigorous testing of competing model hypotheses, partially based on internal processes, include , who tested multiple models for their ability to reproduce in situ soil moisture observations as part of the Distributed Model Intercomparison Project 2 (DMIP2). They found that only 2 out of 16 models provided reasonable estimates of soil moisture. In a similar effort, and also compared modeled soil moisture to in situ observations of soil moisture for a range of hydrological models in different environments. In contrast, and used groundwater observations to test individual components of their models. There are actually relatively few studies that comprehensively quantified differences in internal process representation by simultaneously analyzing multiple models and multiple state and flux variables.

Here, in this model comparison study, we instead use globally available remote-sensing data to evaluate five different model state and flux variables of 12 process-based models with similar overall streamflow performance, which are calibrated by several research groups following an identical protocol. The calibration on streamflow was conducted in our previous study , in which we compared models using hourly streamflow observations, leaving the modeled response of internal processes unused.

In a direct follow-up to the above study, we here hypothesize that process-based models with similar overall streamflow performance rely on similar representations of their internal states and fluxes. We test our hypothesis by simultaneously quantifying the differences in the magnitudes and dynamics of five internal state and flux variables of 12 models, in a comprehensive and systematic way. Our primary aim is to test whether models calibrated to streamflow with similar high-performance levels in reproducing streamflow follow similar pathways to do so, i.e., represent the system in a similar way. A secondary objective is to evaluate the plausibility of model behavior by introducing a set of “soft” measures based on expert knowledge in combination with remote-sensing data of evaporation, snow cover, soil moisture and total water storage anomalies.

2 Study area

We test our hypothesis using data from three catchments in the Belgian Ardennes; all of them are part of the Meuse River basin in northwestern Europe: the Ourthe upstream of Tabreux (ID1), the nested Ourthe Orientale upstream of Mabompré (ID2) and the Semois upstream of Membre-Pont (ID3), as shown in Fig. . The Ardennes Massif and Plateau are underlain by relatively impermeable metamorphic Cambrian rock and Early Devonian sandstone. The pronounced streamflow seasonality of these catchments is driven by high summer and low winter evaporation (defined here as the sum of all evaporation components including transpiration, soil evaporation, interception, sublimation and open water evaporation when applicable), as precipitation is relatively constant throughout the year. Snow is not a major component of the water balance but occurs almost every year with mean annual number of days with precipitation as snow estimated between 35 and 40 d yr${}^{-\mathrm{1}}$ . Even if mean annual snow storage is relatively small, snow can be important for specific events, for example, in 2011, when rain on snow caused widespread flooding in these catchments.

Figure 1

(a) Location of the study catchments in Belgium, northwestern Europe. (b) Digital elevation model and catchments of the Ourthe upstream of Tabreux (ID1), Ourthe Orientale upstream of Mabompré (ID2) and Semois upstream of Membre-Pont (ID3). Pixel size of GRACE, GLEAM, MODIS and SCATSAR-SWI1km. Colored dots are the streamflow gauging locations, and black dots are the precipitation stations.

[Figure omitted. See PDF]

The rain-fed Ourthe River at Tabreux (ID1) is fast-responding due to shallow soils and steep slopes in the catchment. Agriculture is the main land cover (27 % crops and 21 % pasture), followed by 46 % forestry and 6 % urban cover in an area of 1607 km${}^{\mathrm{2}}$ and an elevation ranging between 107 and 663 m . Mean annual precipitation, potential evaporation and streamflow are 979, 730 and 433 mm yr${}^{-\mathrm{1}}$, respectively, for the period 2001–2017.

Nested within the Ourthe catchment (ID1), the Ourthe Orientale upstream of Mabompré (ID2) is characterized by a narrow elevation range from 294 to 662 m, with 65 % of the catchment falling within a 100 m elevation band, making this catchment suitable for analyzing snow processes modeled by lumped models. The Ourthe Orientale upstream of Mabompré has an area of 317 km${}^{\mathrm{2}}$, which corresponds to 20 % of the Ourthe area upstream of Tabreux and has similar land cover fractions. Mean annual precipitation, potential evaporation and streamflow for the period 2001–2017 are also relatively similar, at 1052, 720 and 462 mm yr${}^{-\mathrm{1}}$, respectively.

Forest is the main land cover in the Semois upstream of Membre-Pont (ID3), at 56 %, followed by agriculture (18 % pasture and 21 % crop) and 5 % urban cover. The Semois upstream of Membre-Pont is 24 % smaller than the Ourthe upstream of Tabreux, at 1226 km${}^{\mathrm{2}}$, and elevation ranges between 176 and 569 m. Mean annual precipitation, potential evaporation and streamflow are, respectively, 38 %, 4 % and 46 % higher in the Semois at Membre-Pont, at 1352, 759 and 634 mm yr${}^{-\mathrm{1}}$.

3 Data

3.1 Hydrological and meteorological data

Hourly precipitation gauge data are provided by the Service Public de Wallonie and are spatially interpolated using Thiessen polygons. Daily minimum and maximum temperatures are retrieved from the 0.25${}^{\circ }$ resolution gridded E-OBS dataset and disaggregated to hourly values by linear interpolation using the timing of daily minimum and maximum radiation at Maastricht . Daily potential evaporation is calculated from daily minimum and maximum temperatures using the Hargreaves formula and is disaggregated to hourly values using a sine function during the day and no evaporation at night. We use the same forcing for 2000–2010 as in the previous comparison study and follow the same approach to extend the meteorological dataset for the period 2011–2017. Uncertainty in meteorological data is not explicitly considered, but our primary aim is to compare the models forced with identical data. Observed hourly streamflow data for the Ourthe at Tabreux, Ourthe Orientale at Mabompré and Semois at Membre-Pont are provided by the Service Public de Wallonie for the period 2000–2017.

3.2 Remote-sensing data

3.2.1 GLEAM evaporation

The Global Land Evaporation Amsterdam Model GLEAM, provides daily estimates of land evaporation by maximizing the information recovery on evaporation contained in climate and environmental satellite observations. The equation is used to calculate potential evaporation for bare soil, short canopy and tall canopy land fractions. Actual evaporation is the sum of interception and potential evaporation reduced by a stress factor. This evaporative stress factor is based on microwave observations of vegetation optical depth and estimates of root-zone soil moisture in a multi-layer water-balance model. Interception evaporation is estimated separately using a Gash analytical model and only depends on precipitation and vegetation characteristics. GLEAM v3.3a relies on reanalysis net radiation and air temperature from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 data, satellite and gauge-based precipitation, satellite-based vegetation optical depth, soil moisture and snow water equivalent. The data are available at 0.25${}^{\circ }$ resolution (Fig. b) and account for subgrid heterogeneity by considering three land cover types. We spatially average GLEAM interception and total actual evaporation estimates over the Ourthe catchment upstream of Tabreux for the period 2001–2017.

3.2.2 MODIS snow cover

The Moderate Resolution Imaging Spectroradiometer (MODIS) AQUA (MYD10A1, version 6) and TERRA (MOD10A1, version 6) satellites provide daily maps of the areal fraction of snow cover per 500 m $\times$ 500 m cell (Fig. b) based on the Normalized Difference Snow Index . For each day, AQUA and TERRA observations are merged into a single observation by taking the mean fraction of snow cover per day. The percentage of cells with a fractional snow cover larger than zero and fraction of cells without missing data (i.e., due to cloud cover) for the catchment of the Ourthe Orientale upstream of Mabompré is calculated for each day. For this study, we disregard observations during the summer months (JJA, when temperatures did not drop below 4 ${}^{\circ }$C) and only use daily observations in which at least 40 % of the catchment area has snow cover retrievals not affected by clouds, implying that we have 1463 valid daily observations of mean fractional snow cover. This corresponds to 32 % of all observations of the Ourthe Orientale catchment upstream of Mabompré between 2001 and 2017.

3.2.3 SCATSAR-SWI1km Soil Water Index

SCATSAR-SWI1km is a daily product of soil water content relative to saturation at a 1 km $\times$ 1 km resolution (Fig. b) obtained by fusing spatio-temporally complementary radar sensors . Estimates of the moisture content relative to saturation at various depths in the soil, referred to as the Soil Water Index (SWI), are obtained through temporal filtering of the 25 km METOP ASCAT near-surface soil moisture and 1 km Sentinel-1 near-surface soil moisture . The Soil Water Index features as a single parameter the characteristic time length $T$ . The $T$ value is required to convert near-surface soil moisture observations to estimates of root-zone soil moisture. The $T$ value increases with increasing root-zone storage capacities , resulting in more smoothing and delaying of the near-surface soil moisture signal. The provides the Soil Water Index for $T$ values of 2, 5, 10, 15, 20, 40, 60 and 100 d. Since Sentinel-1 was launched in 2014, the Soil Water Index is available for the period 2015–2017. We calculate the mean soil moisture over all SCATSAR-SWI1km pixels within the Ourthe upstream of Tabreux for the available period.

3.2.4 GRACE total water storage anomalies

The Gravity Recovery and Climate Experiment GRACE, twin satellites, launched in March 2002, measure the Earth's gravity field changes by calculating the changes in the distance between the two satellites as they move one behind the other in the same orbital plane. Monthly total water storage anomalies (in millimeters) relative to the 2004–2009 time-mean baseline are provided at a spatial sampling of 1${}^{\circ }$ (approximately 78 km $\times$ 110 km at the latitude of the study region, Fig. b) by three centers: U. Texas/Center for Space Research (CSR), GeoForschungsZentrum Potsdam (GFZ) and Jet Propulsion Laboratory (JPL). These centers apply different processing strategies which lead to variations in the gravity fields. These gravity fields require smoothing of the noise induced by attenuated short wavelength. The spatial smoothing decreases the already coarse GRACE resolution even further through signal “leakage” of one location to surrounding areas , which increases the uncertainty especially at the relatively small scale of our study catchments. We apply the scaling coefficients provided by NASA to restore some of the signal loss due to processing of GRACE observations . The data of the three processing centers are each spatially averaged over the catchments of the Ourthe upstream of Tabreux and the Semois upstream of Membre-Pont for the period April 2002 to February 2017.

3.3 Data uncertainty

The hydrological evaluation data are all subject to uncertainties . Streamflow is not measured directly but depends on water level measurements and a rating curve. quantify median streamflow uncertainties of $\pm$12 %, $\pm$24 % and $\pm$34 % for average, high and low streamflow conditions, respectively, using a Monte Carlo sampling approach of multiple feasible rating curves for 43 UK catchments. We sample from these uncertainty ranges to transform the streamflow observations (100 realizations). We then quantify signature uncertainty originating from streamflow data uncertainty using the 100 sampled time series for a selection of streamflow signatures (Sect. ). The 5th–95th uncertainty bounds of median annual streamflow, baseflow and flashiness indices result in $\pm$11 %, $\pm$9 % and $\pm$12 %, respectively. These magnitudes are similar to those reported by .

GLEAM evaporation estimates are inferred from models and forcing data which are all affected by uncertainty. Yet uncertainty estimates of GLEAM evaporation are not available. However, GLEAM evaporation was evaluated against FLUXNET data by . For the nearby station of Lonzee in Belgium, they report similar annual rates and a high correlation coefficient of 0.91 between the daily time series. GLEAM mean annual evaporation was compared to the ensemble mean of five evaporation datasets in and shows higher than average values in Europe (of approximately 60 mm yr${}^{-\mathrm{1}}$ or 10 % of mean annual rates for our study area). The partitioning of evaporation into different components (transpiration, interception and soil evaporation) differs substantially between different evaporation datasets, as shown by . GLEAM interception currently only considers tall vegetation and underestimates in situ data and is $\sim$ 50 % lower than estimates from other datasets . These uncertainties underline that GLEAM (and other remote-sensing data) cannot be considered a reliable representation of real-world quantities. However, the comparison of daily dynamics and absolute values of this independent data source with modeled results is still valuable for detecting potential outliers and understanding their behavior. In addition, the different methods used to estimate potential evaporation of GLEAM and our model forcing should not impede us from testing the consistency between the resulting actual evaporation .

The most frequent errors within the MODIS snow cover products are due to cloud–snow discrimination problems. Daily MODIS snow maps have an accuracy of approximately 93 % at the pixel scale, with lower accuracy in forested areas and complex terrain and when snow is thin and ephemeral and higher accuracy in agricultural areas . However, here, MODIS data are used to estimate the number of days with snow at the catchment scale. We expect lower classification errors at the catchment scale as it would require many pixels to be misclassified at the same time. For each day and each pixel of valid MODIS observations, we sample from a binomial distribution with a probability of 93 % that MODIS is correct when the pixel is classified as snow and assume a higher probability of 99 % that MODIS is correct when the pixel is classified as no-snow to prevent overestimation of snow for days without snow . We repeat the experiment 1000 times in a Monte Carlo procedure. This results in less than $\pm$2 % uncertainty in the number of days when MODIS observes snow at the catchment scale.

The soil water content relative to saturation of SCATSAR-SWI1km is estimated from observed radar backscatter through a change detection approach, which interprets changes in backscatter as changes in soil moisture, while other surface properties are assumed static . The degree of saturation of the near surface is given in relative units from 0 % (dry reference) to 100 % (wet reference) and is converted to deeper layers through an exponential filter called the Soil Water Index. The smoothing and delaying effect of the Soil Water Index narrows the range of the near-surface degree of saturation. Therefore, data matching techniques are often used to rescale satellite data to match the variability of modeled or observed data , which suggests the difficulty in meaningfully comparing the range of modeled and remote-sensing estimates of root-zone soil moisture content relative to saturation. However, the dynamics of SCATSAR-SWI1km data have been evaluated against in situ stations of the International Soil Moisture Network, despite commensurability issues of comparing in situ point measurements and areal satellite data. Spearman rank correlation coefficients of 0.56 are reported for $T$ values up to 15 d and 0.43 for a $T$ value of 100 d .

GRACE estimates of total water storage anomalies suffer from signal degradation due to measurement errors and noise. Filtering approaches are applied to reduce these errors but induce leakage of signals from surrounding areas. The uncertainty decreases as the size of the region under consideration increases. However, time series of a single pixel may still be used to compare dynamics and amplitudes of total water storage anomalies despite possible large uncertainty . We estimate an uncertainty in total water storage anomalies of $\sim$ 18 mm in the pixels of our catchments by combining measurement and leakage errors in quadrature, which are both provided for each grid location .

4 Methods

4.1 Models and protocol

Eight research groups (Wageningen University, Université de Lorraine, Leuven University, Delft University of Technology, Deltares, Irstea (now INRAE), Eawag and Flanders Hydraulics Research) participated in the comparison experiment and applied one or several hydrological models (Fig. ). The models include WALRUS Wageningen Lowland Runoff Simulator,, PRESAGES PREvision et Simulation pour l'Annonce et la Gestion des Etiages Sévères,, VHM Veralgemeend conceptueel Hydrologisch Model,, FLEX-Topo, which was still under development when it was calibrated for our previous study , a distributed version of the HBV model Hydrologiska Byråns Vattenbalansavdelning,, the SUPERFLEX M2 to M5 models , dS2 distributed simple dynamical systems,, GR4H Génie Rural à 4 paramètres Horaire, combined with the CemaNeige snow module and NAM NedborAfstrommings Model,. The main differences and similarities between the models in terms of snow processes, root-zone storage, total storage and evaporation processes are summarized in Tables –.

Figure 2

Simplified schematic overview of 12 model structures adapted from with the aim of highlighting similarities and differences between the models. Solid arrows indicate fluxes between stores, while dashed arrows indicate the influence of a state on a flux. Colored arrows indicate incoming or outgoing fluxes, whereas black arrows indicate internal fluxes. The narrow blue rectangle in GR4H indicates the presence of an interception module without interception storage capacity (Table ). Storages with a color gradient indicate the combination of several components in one reservoir. FLEX-Topo consists of three hydrological response units connected through a shared slow reservoir, and wflow_hbv is a distributed model.

[Figure omitted. See PDF]

Table 1

Description of symbols used for fluxes, storages and parameters in Tables  and .

Symbol	Unit	Description
Fluxes
$E_{\mathrm{P}}$	mm h${}^{-\mathrm{1}}$	Potential evaporation
$E_{\mathrm{I}}$	mm h${}^{-\mathrm{1}}$	Interception evaporation
$E_{\mathrm{R}}$	mm h${}^{-\mathrm{1}}$	Transpiration and soil evaporation
$E_{\mathrm{W}}$	mm h${}^{-\mathrm{1}}$	Sublimation
$E_{\mathrm{A}}$	mm h${}^{-\mathrm{1}}$	Total actual evaporation (sum of soil evaporation, transpiration, (separate) interception and, if applicable, sublimation)
$P$	mm h${}^{-\mathrm{1}}$	Precipitation
$P_{\mathrm{R}}$	mm h${}^{-\mathrm{1}}$	Precipitation entering the root-zone storage (after snow and/or interception if present or fraction/total precipitation)
$Q$	mm h${}^{-\mathrm{1}}$	Streamflow
$Q_{\mathrm{R}}$	mm h${}^{-\mathrm{1}}$	Flux from root zone to fast and/or slow runoff storage
$Q_{\mathrm{P}}$	mm h${}^{-\mathrm{1}}$	Percolation flux from root-zone storage to slow runoff storage
$Q_{\mathrm{C}}$	mm h${}^{-\mathrm{1}}$	Capillary flux from slow runoff storage to root-zone storage
$Q_{\mathrm{G}}$	mm h${}^{-\mathrm{1}}$	Seepage (up/down)/extraction
Storages
$S_{\mathrm{T}}$	mm	Total storage
$S_{\mathrm{W}}$	mm	Snow storage
$S_{\mathrm{I}}$	mm	Interception storage
$S_{\mathrm{R}}$	mm	Root-zone storage
${\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}}$	–	Relative root-zone storage ($S_{\mathrm{R}}/S_{\mathrm{R},\max }$)
$S_{\mathrm{D}}$	mm	Storage deficit
$S_{\mathrm{VQ}}$	mm	Very quick runoff storage
$S_{\mathrm{F}}$	mm	Fast runoff storage
$S_{\mathrm{S}}$	mm	Slow runoff storage
$S_{\mathrm{SW}}$	mm	Surface water storage
Parameters
$C_{\mathrm{E}}$	–	Correction factor for $E_{\mathrm{P}}$
$I_{\max }$	mm	Maximum interception capacity
$S_{\mathrm{R},\max }$	mm	Maximum root-zone storage capacity
$S_{\mathrm{thresh}}$	mm	Threshold of root-zone storage above which $E_{\mathrm{R}}=E_{\mathrm{P}}$
$L_{\mathrm{P}}$	–	Threshold of relative root-zone storage above which $E_{\mathrm{R}}=E_{\mathrm{P}}$
$C_{\mathrm{cst}}$	–	Constant water stress coefficient to estimate $E_{\mathrm{R}}$
$a$, $b$, $S_{\mathrm{0}}$	–	Parameters describing the shape of the streamflow sensitivity
$a_{\mathrm{S}}$	–	Fraction of land surface covered by surface water
$a_{\mathrm{G}}$	–	Fraction of land surface not covered by surface water

In our previous study , we defined a modeling protocol to limit the degrees of freedom in the modeling decisions of the individual participants , allowing us to meaningfully compare the model results. The protocol involved forcing the models with the same input data and calibrating them for the same time period, using the same objective functions. However, participants were free to choose a parameter search method, as we considered it to be part of the modelers' experience with the model, even if this would make comparison less straightforward. The models were previously calibrated using streamflow of the Ourthe at Tabreux (ID1) for the period 2004 to 2007, using 2003 as a warm-up year . The Nash–Sutcliffe efficiencies of the streamflow and the logarithms of the streamflow were simultaneously used as objective functions to select an ensemble of 20 feasible parameter sets to account for parameter uncertainty and ensure a balance between the models' ability to reproduce both high and low flows. The temporal and spatial transferability of the models was tested by evaluating the models in a pre- and post-calibration period and by applying the calibrated model parameter sets to nested and neighboring catchments, including catchment ID2 and ID3 . Results thereof are presented in .

In the current study, we run the calibrated models for an additional period from 2011 to 2017 for the Ourthe at Tabreux (ID1), the Ourthe Orientale at Mabompré (ID2) and the Semois at Membre-Pont (ID3). Note that the calibration (2004–2007) and post-calibration (2008–2017) periods have relatively similar hydro-climatic characteristics in terms of streamflow and overall water balance (Figs. S1 and S2 of the Supplement). The modeling groups have provided simulation results for each catchment in terms of streamflow, groundwater losses/gains, interception evaporation, root-zone evaporation (transpiration and soil evaporation), total actual evaporation, snow storage, root-zone storage and total storage as a sum of all model storage volumes (Table ) at an hourly time step for the total period 2001–2017. We compare these modeled states and fluxes and evaluate them against their remote-sensing equivalents as further explained in Sects.  and . Our results are mainly shown for the Ourthe at Tabreux (ID1), as this catchment was used for calibration of the models. However, the snow analysis is performed for catchment ID2 due to the narrower elevation range. Catchments ID1 and ID3 are used to analyze the spatial variability of total storage anomalies.

Table 2

Number of calibrated model parameters, spatial distribution, and model performance calculated for the period 2008–2017 with the Euclidean distance where a value of 0 would indicate a perfect model. Main characteristics describing snow storage, root-zone storage and total storage per model. Notations are defined in Table .

	GR4H	M5	NAM	wflow_hbv	dS2	M4	M3	M2	PRESAGES	FLEX-Topo	VHM	WALRUS
Number of calibrated parameters	4	9	12	9	4	7	6	5	6	20	12	3
Lumped (L)/semi-distributed (S)/distributed (D)	L	L	L	D	L	L	L	L	L	S	L	L
Euclidean distance $\sqrt{(\mathrm{1}-E_{\mathrm{NS},\mathrm{Q}}{)}^{\mathrm{2}}+(\mathrm{1}-E_{\mathrm{NS},\mathrm{logQ}}{)}^{\mathrm{2}}}$	0.17	0.18	0.18	0.20	0.21	0.23	0.23	0.24	0.24	0.26	0.26	0.34
Snow storage $S_{\mathrm{W}}$ (compared to MODIS snow cover)
Snow module	✓	✓	✓	✓	–	✓	–	–	–	✓	–	✓
Degree-hour method	✓	✓	✓	✓	–	✓	–	–	–	✓	–	✓
Elevation zones	✓	–	–	✓	–	–	–	–	–	✓	–	–
Temperature interval for rainfall and snow	✓	–	–	✓	–	–	–	–	–	✓	–	✓
Melt factor constant in time	–	✓	–	✓	–	✓	–	–	–	✓	–	✓
Melt factor $\sim$ snow storage	✓	–	✓	–	–	–	–	–	–	–	–	–
Refreezing of liquid water	–	–	✓	✓	–	–	–	–	–	–	–	–
Sublimation	–	–	–	–	–	–	–	–	–	✓	–	–
Calibration snow parameters	–	✓	✓	–	–	✓	–	–	–	✓	–	–
Root-zone storage $S_{\mathrm{R}}$ (compared to SCATSAR-SWI1km Soil Water Index)
Separate root-zone module with capacity $S_{\mathrm{R},\max }$	✓	✓	✓	✓	–	✓	✓	✓	✓	✓	✓	–
$\frac{\mathrm{d}{S}_{\mathrm{R}}}{\mathrm{d}t}=P_{\mathrm{R}}-E_{\mathrm{R}}$	–	–	–	–	–	–	–	–	–	–	✓	–
$\frac{\mathrm{d}{S}_{\mathrm{R}}}{\mathrm{d}t}=P_{\mathrm{R}}-E_{\mathrm{R}}+Q_{\mathrm{C}}$	–	–	✓	–	–	–	–	–	–	–	–	–
$\frac{\mathrm{d}{S}_{\mathrm{R}}}{\mathrm{d}t}=P_{\mathrm{R}}-E_{\mathrm{R}}-Q_{\mathrm{R}}$	✓	✓	–	✓	–	✓	✓	✓	✓	–	–	–
$\frac{\mathrm{d}{S}_{\mathrm{R}}}{\mathrm{d}t}=P_{\mathrm{R}}-E_{\mathrm{R}}-Q_{\mathrm{R}}--Q_{\mathrm{P}}+Q_{\mathrm{C}}$	–	–	–	–	–	–	–	–	–	✓	–	–
Total storage $S_{\mathrm{T}}$ (anomalies are compared to GRACE total storage anomalies)
$S_{\mathrm{T}}=-S_{\mathrm{D}}\cdot a_{\mathrm{G}}+S_{\mathrm{F}}\cdot a_{\mathrm{G}}+S_{\mathrm{SW}}\cdot a_{\mathrm{S}}$	–	–	–	–	–	–	–	–	–	–	–	✓
$S_{\mathrm{T}}\left(Q\right)=\frac{\mathrm{1}}{a}\frac{\mathrm{1}}{\mathrm{1}-b}\cdot Q^{\mathrm{1}-b}+S_{\mathrm{0}}$	–	–	–	–	✓	–	–	–	–	–	–	–
$S_{\mathrm{T}}=S_{\mathrm{R}}+S_{\mathrm{F}}$	–	–	–	–	–	–	✓	✓	–	–	–	–
$S_{\mathrm{T}}=S_{\mathrm{W}}+S_{\mathrm{R}}+S_{\mathrm{F}}$	–	–	–	–	–	✓	–	–	–	–	–	–
$S_{\mathrm{T}}=S_{\mathrm{W}}+S_{\mathrm{R}}+S_{\mathrm{F}}+S_{\mathrm{S}}$	–	✓	–	–	–	–	–	–	–	–	–	–
$S_{\mathrm{T}}=S_{\mathrm{R}}+S_{\mathrm{VQ}}+S_{\mathrm{F}}+S_{\mathrm{S}}$	–	–	–	–	–	–	–	–	✓	–	✓	–
$S_{\mathrm{T}}=S_{\mathrm{W}}+S_{\mathrm{R}}+S_{\mathrm{VQ}}+S_{\mathrm{F}}+S_{\mathrm{S}}$	✓	–	–	–	–	–	–	–	–	–	–	–
$S_{\mathrm{T}}=S_{\mathrm{I}}+S_{\mathrm{W}}+S_{\mathrm{R}}+S_{\mathrm{F}}+S_{\mathrm{S}}$	–	–	–	✓	–	–	–	–	–	–	–	–
$S_{\mathrm{T}}=S_{\mathrm{I}}+S_{\mathrm{W}}+S_{\mathrm{R}}+S_{\mathrm{VQ}}+S_{\mathrm{F}}+S_{\mathrm{S}}$	–	–	✓	–	–	–	–	–	–	✓	–	–

4.2 Model evaluation: water balance

All the models are evaluated in terms of the long-term water balance, which indicates the partitioning between drainage and evaporative fluxes and allows us to assess long-term conservation of water and energy. We compare mean annual streamflow with observations and mean annual actual evaporation and interception evaporation with GLEAM estimates for the Ourthe at Tabreux during the evaluation period 2008–2017. A detailed description of streamflow performance for specific events (low and high flows, snowmelt events, transition from dry to wet periods) has been detailed in the previous paper . In the current study, differences in streamflow dynamics are briefly summarized by assessing observed and modeled baseflow indices $I_{\mathrm{baseflow}}$, and flashiness indices $I_{\mathrm{flashiness}}$,, as these are representative of the partitioning of drainage into fast and slow responses. Seasonal dynamics of actual evaporation over potential evaporation and runoff coefficients during winter (October–March) and summer (April–September) are compared between the models.

4.3 Model evaluation: internal states

We compare modeled snow storage, root-zone soil moisture and total storage between the models and with remote-sensing estimates of MODIS snow cover, SCATSAR-SWI1km Soil Water Index and GRACE total storage anomalies, respectively, as shown in Tables –.

Table 3

Main characteristics describing evaporation processes per model (where ✓${}^{\mathrm{1}}$ indicates $L_{\mathrm{P}}=\mathrm{1}$ and ✓${}^{\mathrm{2}}$ indicates $E_{\mathrm{I}}=\mathrm{0}$). Notations are defined in Table .

	GR4H	M5	NAM	wflow_hbv	dS2	M4	M3	M2	PRESAGES	FLEX-Topo	VHM	WALRUS
Correction factor for potential evaporation	–	✓	–	–	–	✓	✓	✓	–	–	–	–
Interception evaporation $E_{\mathrm{I}}$	✓	–	✓	✓	–	–	–	–	–	✓	–	–
Maximum interception storage $I_{\max }$	–	–	✓	✓	–	–	–	–	–	✓	–	–
$I_{\max }\sim \mathrm{1.1}-\mathrm{3.4}$ mm	–	–	–	–	–	–	–	–	–	✓	–	–
$I_{\max }\sim \mathrm{1.4}-\mathrm{2.9}$ mm	–	–	–	✓	–	–	–	–	–	–	–	–
$I_{\max }\sim \mathrm{5.3}--\mathrm{6.9}$ mm	–	–	✓	–	–	–	–	–	–	–	–	–
$E_{\mathrm{I}}=\left\{\begin{array}{ll}{E}_{\mathrm{P}},& \text{if}{S}_{\mathrm{I}}>\mathrm{0}.\\ \mathrm{0},& \text{otherwise}.\end{array}\right.$	–	–	✓	✓	–	–	–	–	–	✓	–	–
$E_{\mathrm{I}}=\left\{\begin{array}{ll}{E}_{\mathrm{P}},& \text{if}P>E_{\mathrm{P}}.\\ P,& \text{otherwise}.\end{array}\right.$	✓	–	–	–	–	–	–	–	–	–	–	–
Transpiration and soil evaporation $E_{\mathrm{R}}$	✓	✓	✓	✓	✓	✓	✓	✓	✓	✓	✓	✓
$E_{\mathrm{R}}=E_{\mathrm{P}}\cdot C_{\mathrm{cst}}$	–	–	–	–	✓	–	–	–	–	–	–	–
$E_{\mathrm{R}}=E_{\mathrm{P}}\cdot \frac{{\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}}\cdot (\mathrm{2}-{\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}})}{\mathrm{1}+E_{\mathrm{P}}/S_{\mathrm{R},\max }\cdot (\mathrm{2}-{\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}})}$	✓	–	–	–	–	–	–	–	✓	–	–	–
$E_{\mathrm{R}}=E_{\mathrm{P}}\cdot C_{\mathrm{E}}\cdot \frac{{\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}}\cdot (\mathrm{1}+m_{\mathrm{1}})}{{\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}}+m_{\mathrm{1}}}$, with $m_{\mathrm{1}}={\mathrm{10}}^{-\mathrm{2}}$	–	✓	–	–	–	✓	✓	✓	–	–	–	–
$E_{\mathrm{R}}=\left\{\begin{array}{ll}(E_{\mathrm{P}}-E_{\mathrm{I}})\cdot \frac{{\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}}}{L_P},& \text{if}{\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}}<L_P.\\ E_{\mathrm{P}}-E_{\mathrm{I}},& \text{otherwise}.\end{array}\right.$	–	–	✓${}^{\mathrm{1}}$	✓	–	–	–	–	–	✓	✓${}^{\mathrm{2}}$	–
$E_{\mathrm{R}}=E_{\mathrm{P}}\cdot \mathrm{f}\left(S_{\mathrm{d}}\right)$	–	–	–	–	–	–	–	–	–	–	–	✓
Total actual evaporation $E_{\mathrm{A}}$	✓	✓	✓	✓	✓	✓	✓	✓	✓	✓	✓	✓
$E_{\mathrm{A}}=E_{\mathrm{R}}$	–	✓	–	–	✓	✓	✓	✓	✓	–	✓	✓
$E_{\mathrm{A}}=E_{\mathrm{R}}+E_{\mathrm{I}}$	✓	–	✓	✓	–	–	–	–	–	–	–	–
$E_{\mathrm{A}}=E_{\mathrm{R}}+E_{\mathrm{I}}+E_{\mathrm{W}}$	–	–	–	–	–	–	–	–	–	✓	–	–

4.3.1 Snow days

As most models used in this study are lumped, it is not possible to spatially evaluate modeled snow cover versus MODIS snow cover. However, we can classify each day in a binary way according to the occurrence of snow, based on a threshold for the percentage of cells in the catchment where snow cover is detected. MODIS snow cover observations are classified as days with and without snow using thresholds of both 10 % and 15 % of snow-covered cells in the catchment to be counted as a day with snow, in a sensitivity analysis. For each model, snow days are distinguished from non-snow days whenever the water stored as snow is above 0.05 mm to account for numerical rounding. For each model (and each retained parameter set), we then compare whether modeled snow coincides with “truly” observed snow by MODIS, for each day with a valid MODIS observation. We create a confusion matrix with counts of true positives when observations and model results agree on the presence of snow (hits), false positives when the model indicates the presence of snow but this is not observed by MODIS (false alarms), false negatives when the model misses the presence of snow observed by MODIS (miss) and true negatives when observations and model results agree on the absence of snow (correct rejections). This allows us to identify the trade-off between, on the one hand, the miss rate between model and remote-sensing observation, as the ratio of misses over actual positives (number of days when snow is observed by MODIS) and, on the other hand, the false discovery rate as the ratio of false alarms over predicted positives (number of days when snow is modeled). We also compare annual maximum snow storage and number of days with snow between the seven models with a snow module (GR4H, M5, NAM, wflow_hbv, M4, FLEX-Topo, WALRUS). The snow analysis is performed in the catchment of the Ourthe Orientale upstream of Mabompré as it features the narrowest elevation range among the study catchments (i.e., 294–662 m a.s.l. versus 108–662 m for the Ourthe upstream of Tabreux) and thus plausibly permits a lumped representation of the snow component.

4.3.2 Root-zone soil moisture

We compare the range of the relative root-zone soil moisture ${\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}}$ (${\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}}=S_{\mathrm{R}}/S_{\mathrm{R},\max }$, Table ) between the models for the period in which SCATSAR-SWI1km is available (2015–2017). Time series of catchment-scale root-zone soil moisture are available for all the models except WALRUS and dS2 as these models have a combined soil reservoir (Fig. ). The dS2 model only relies on the sensitivity of streamflow to changes in total storage. In WALRUS, the state of the soil reservoir (which includes the root zone) is expressed as a storage deficit and is therefore not bound by an upper limit (Table ). Root-zone storage capacities ($S_{\mathrm{R},\max }$, mm) are available as a calibration parameter for all the other models. We relate the range in relative root-zone soil moisture to the maximum root-zone storage capacity $S_{\mathrm{R},\max }$, because we expect models with small root-zone storage capacities $S_{\mathrm{R},\max }$ to entirely utilize the available storage through complete drying and saturation.

We then compare the similarity of the dynamics of modeled time series of the relative root-zone soil moisture with the remotely sensed SCATSAR-SWI1km Soil Water Index for several values of the characteristic time length parameter ($T$ in days). The $T$ value has previously been positively correlated with root-zone storage capacity, assuming a high temporal variability of root-zone soil moisture and therefore a low $T$ value for small root-zone storage capacities $S_{\mathrm{R},\max }$ . For each model and feasible realization, we identify the $T$ value that yields the highest Spearman rank correlation between modeled root-zone soil moisture and Soil Water Index. We then relate the optimal $T$ value to the root-zone storage capacity $S_{\mathrm{R},\max }$. This analysis enables us to identify potential differences in the representation and the dynamics of root-zone storage between the models.

4.3.3 Total storage anomalies

For each model, we calculate time series of total storage (Table ) and mean monthly total storage anomalies relative to the 2004–2009 time-mean baseline for comparison with GRACE estimates for the Ourthe upstream of Tabreux (ID1) and the Semois upstream of Membre-Pont (ID3). Both catchments coincide with two neighboring GRACE cells, allowing us to test how well the models reproduce the observed spatial variability. We further relate the modeled range of total storage (maximum minus minimum total storage over the time series) to Spearman rank correlation coefficients between modeled and GRACE estimates of total storage anomalies.

4.4 Interactions between storage and fluxes during dry periods

The impact of a relatively large ($>$ 200 mm) versus relatively small ($<$ 150 mm) root-zone storage capacity on actual evaporation, streamflow and total storage is assessed during a dry period in September 2016 by selecting two representative models with high streamflow model performance (GR4H and M5). The plausibility of the hydrological response of these two model representations is evaluated against remote-sensing estimates of root-zone soil moisture and actual evaporation.

4.5 Plausibility of process representations

The models are subsequently ranked and evaluated in terms of their consistency with observed streamflow, remote-sensing data and expert knowledge with due consideration of the uncertainty in the evaluation data, as detailed in Sect. . We summarize our main findings by evaluating the models in terms of their deviations around median annual streamflow, flashiness and baseflow indices, median annual actual evaporation and interception compared to GLEAM estimates, the number of days with snow over valid MODIS observations, the number of days per year with empty root-zone storage and the very dry total storage anomalies compared to GRACE estimates.

5 Results and discussion

5.1 Water balance

5.1.1 Streamflow

All the models show high Nash–Sutcliffe efficiencies of the streamflow and the logarithm of the streamflow ($E_{\mathrm{NS},\mathrm{Q}}$ and $E_{\mathrm{NS},\mathrm{logQ}}$) with median values of above 0.7 for the post-calibration evaluation period 2008–2017 (Fig. a and Table  for the calculation of the Euclidean distances). The interannual variability of streamflow agrees strongly with observations for each model in the period 2008–2017 (Fig. b). The difference between modeled and observed median streamflow varies between $-\mathrm{5.6}$ % and 5.6 %, and the difference in total range varies between $-\mathrm{9.6}$ % and 20 %. This is in line with our results in the previous paper, in which we also showed that all the models perform well in terms of commonly used metrics . However, there are differences in the partitioning of fast and slow runoff, as shown by the flashiness and baseflow indices ($I_{\mathrm{flashiness}}$ and $I_{\mathrm{baseflow}}$) in Fig. c. The largest underestimation of the flashiness index occurs for M2 and dS2 ($\sim$ 20 %), while FLEX-Topo shows the highest overestimation (26 %). FLEX-Topo and WALRUS underestimate the baseflow index most (41 % and 70 %, respectively), while GR4H and M5 show the highest overestimation (15 % and 21 %, respectively). There is a strong similarity between modeled and observed hydrographs for one of the best performing models M5, as quantified by its low Euclidean distance (Fig. d and Table ). The GR4H model is the only one which includes deep groundwater losses, but they are very limited and represent only 1.6 % of total modeled streamflow of the Ourthe at Tabreux, or approximately 7 mm yr${}^{-\mathrm{1}}$.

Figure 3

Evaluation of modeled streamflow performance for the Ourthe at Tabreux for the period 2008–2017. (a) Nash–Sutcliffe efficiencies of the streamflow $E_{\mathrm{NS},\mathrm{Q}}$ and the logarithm of the streamflow $E_{\mathrm{NS},\mathrm{logQ}}$ (median, 25th/75th percentiles across parameter sets). (b) Modeled mean annual streamflow for hydrological years between 2008 and 2017 across feasible parameter sets. The models are ranked from the highest to the lowest performance according to the Euclidean distance of streamflow performance (see Table ). The dashed line and grey shaded areas show median, 25th/75th and minimum–maximum range of observed mean annual streamflow. (c) Baseflow index $I_{\mathrm{baseflow}}$ as a function of the flashiness index $I_{\mathrm{flashiness}}$ (median, 25th/75th percentiles across parameter sets). Observed values are shown by the grey dashed lines. (d) Observed and modeled hydrographs of model M5 with high streamflow model performance (low Euclidean distance).

[Figure omitted. See PDF]

5.1.2 Actual evaporation

Modeled median annual actual evaporation $E_{\mathrm{A}}$ (computed as the sum of soil evaporation, transpiration, (separate) interception evaporation and, if applicable, sublimation, Table ) for hydrological years between October 2008 and September 2017 varies between 507 and 707 mm yr${}^{-\mathrm{1}}$ across models, with a median of 522 mm yr${}^{-\mathrm{1}}$, which is approximately 10 % lower than the GLEAM estimate of 578 mm yr${}^{-\mathrm{1}}$, as shown in Fig. a. Annual actual evaporation of the VHM model is very high compared to the other models, with a median of 707 mm yr${}^{-\mathrm{1}}$, and approximates potential evaporation (median of 732 mm yr${}^{-\mathrm{1}}$). Calibration of the VHM model is meant to follow a manual stepwise procedure including the closure of the water balance during the identification of soil moisture processes . However, in the automatic calibration prescribed by the current protocol, this step was not performed, which explains the unusual high actual evaporation in spite of relatively similar annual streamflow compared to the other models, as there is no closure of the water balance (Fig. a).

Figure 4

Evaluation of modeled evaporation for the Ourthe upstream of Tabreux for the period 2008–2017. (a) Modeled mean annual actual evaporation $E_{\mathrm{A}}$ and minimum–maximum range of mean annual interception evaporation $E_{\mathrm{I}}$ for hydrological years between 2008 and 2017 across feasible parameter sets. The dark grey shaded area shows the range of potential evaporation $E_{\mathrm{P}}$ used as input for the models. The light grey shaded area shows GLEAM actual and interception evaporation. (b) Daily actual evaporation from GLEAM and modeled by the M5 model. (c) Summer against winter $E_{\mathrm{A}}$/$E_{\mathrm{P}}$ ratios for each model (median and 25th/75th percentiles across parameter sets). (d) Summer against winter runoff coefficient $Q/P$ for each model (median and 25th/75th percentiles across parameter sets), plotted on the same scale. The dashed grey lines indicate the observed median runoff coefficients in summer and winter.

[Figure omitted. See PDF]

Interception evaporation is included in four models, with GR4H showing the lowest annual interception evaporation of 100 mm yr${}^{-\mathrm{1}}$ (19 % of $E_{\mathrm{A}}$ or 10 % of $P$), FLEX-Topo and wflow_hbv having relatively similar amounts of approximatively 250 mm yr${}^{-\mathrm{1}}$ ($\sim$ 45 % of $E_{\mathrm{A}}$ or 26 % of $P$) and NAM having the highest annual interception evaporation of 340 mm yr${}^{-\mathrm{1}}$ (65 % of $E_{\mathrm{A}}$ or 36 % of $P$), as shown in Fig. a. Differences are related to the presence and maximum size of the interception storage ($I_{\max }$), as shown in Table . GLEAM interception estimates of 189 mm yr${}^{-\mathrm{1}}$ are almost twice as high as GR4H estimates, 25 % lower than FLEX-Topo and wflow_hbv, and 44 % lower than NAM values, suggesting a large uncertainty in the contribution of interception and transpiration to actual evaporation. For comparison, measurements of the fraction of interception evaporation over precipitation in forested areas vary significantly depending on the site location, with estimates of 37 % for a Douglas fir stand in the Netherlands , 27 %, 32 % and 42 % for three coniferous forests of Great Britain and 50 % for a forest in Puerto Rico , and are difficult to extrapolate to other catchments due to the heterogeneity and complexity of natural systems.

GLEAM estimates of actual evaporation show relatively high evaporation rates in winter and are never reduced to zero in summer, as opposed to modeled M5 estimates, as shown in Fig. b. GLEAM actual evaporation minus the separately calculated interception is 94 % of potential evaporation, implying almost no water-limited conditions, as opposed to our models in which actual evaporation in summer (April–September) is, due to water stress, reduced to approximately 73 % of potential evaporation on average for all the models except VHM (Fig. c). Larger differences between the models occur in the ratio $E_{\mathrm{A}}/E_{\mathrm{P}}$ during winter (October–March), when FLEX-Topo, wflow_hbv and VHM show $E_{\mathrm{A}}/E_{\mathrm{P}}$ ratios close to unity and dS2 the lowest values of $E_{\mathrm{A}}/E_{\mathrm{P}}\sim$ 0.75 as shown in Fig. c. The dS2 model differs from all the other models as it relies on a year-round constant water stress coefficient ($C_{\mathrm{cst}}$), independent of water supply, while the stress coefficient depends on root-zone soil moisture content in all the other models (Table ).

Most models slightly overestimate summer runoff coefficients, with values between 0.22 and 0.26, which are very close to the observed value of 0.22, as shown in Fig. d. During winter, runoff coefficients vary between 0.55 and 0.71, which is close to the observed value of 0.66. This implies a relatively high level of agreement between the models in reproducing the medium- to long-term partitioning of precipitation into evaporation and drainage and thus in approximating at least long-term conservation of energy .

5.2 Internal model states

5.2.1 Snow days

MODIS snow cover is detected over most of the catchment area for some time each year between November 2001 and November 2017, except for the periods of November 2006 to March 2007 and November 2007 to March 2008, when snow is detected in less than half of the catchment cells, as shown in Fig. a. The number and magnitude of modeled snow storage events varies between the models (Fig. b). The modeled number of snow days per hydrological year varies from $\sim$ 28 d for FLEX-Topo, WALRUS and wflow_hbv to $\sim$ 62 d for GR4H and $\sim$ 90 d for NAM, M4 and M5, as shown in Fig. c. The variability in median annual maximum snow storage varies from 3 mm for wflow_hbv and $\sim$ 5–6 mm for FLEX-Topo and WALRUS to $\sim$ 10 mm for GR4H, M4, and M5 and 15 mm for NAM. We further evaluate the plausibility of these modeled snow processes by comparing modeled and observed snow cover for days when a valid MODIS observation is available.

Figure 5

(a) Fraction of cells with a MODIS areal fraction snow cover greater than zero in the Ourthe Orientale upstream of Mabompré for the period 2001–2017. MODIS data are available once every 3 d on average. The dashed lines show the two thresholds of 10 % and 15 % selected to distinguish snow days. (b) Modeled snow storage for two contrasting models M5 and WALRUS for the light orange shaded period. (c) Median annual maximum snow storage as a function of number of days per year with snow. Light (yellowish) colors indicate models with higher performance (lower Euclidean distances). The vertical and horizontal error bars indicate the 25th/75th percentiles over time and feasible parameter sets. (d, e) Two-dimensional representation of the false discovery rate as a function of the miss rate when applying thresholds of (d) 10 % and (e) 15 % of cells within the catchment with snow cover greater than zero. In this representation, the perfect model would be at the origin (0 % misses and 0 % false alarms). The dotted lines show the distance from the origin. The vertical and horizontal error bars indicate the uncertainty within feasible parameter sets.

[Figure omitted. See PDF]

The presence of snow modeled by FLEX-Topo, wflow_hbv and WALRUS coincides for 92 % with the presence of snow observed by MODIS. However, these models fail to model snow for $\sim$ 62 % of days when MODIS reports the presence of snow, implying that these models miss many observed snow days, but when they predict snow, it was also observed (Fig. d).

NAM, M4 and M5, on the other hand, predict the presence of snow which coincides with snow observed by MODIS in $\sim$ 68 % of the positive predictions, implying a relatively high probability of false alarm snow prediction of $\sim$ 32 %. However, they miss only $\sim$ 29 % of actual positive snow days observed by MODIS (Fig. d). This suggests that these models miss fewer observed snow days, but they also overpredict snow day numbers, which could be related to the use of a single temperature threshold to distinguish between snow and rain, as opposed to a temperature interval in the other models (Table ).

GR4H is in between the two previously mentioned model categories, with a snow prediction which coincides with observed snow by MODIS in 79 % of the positive predictions and therefore only 21 % of false alarms. The model misses 42 % of actual positive snow days observed by MODIS. GR4H therefore shows a more balanced trade-off between the number of false alarms and the number of observed snow events missed. This is illustrated in Fig. d.

With an increased threshold to distinguish snow days in MODIS, from 10 % to 15 % of cells in the catchment with a detected snow cover (Fig. d and e, respectively), we decrease the number of observed snow days. For all the models, this leads to an increase in the ratio of false alarms over predicted snow days but also a decrease in the ratio of missed events over actual snow days observed by MODIS. However, as all the models are similarly affected by the change in threshold, our findings on the differences in performance between the models show little sensitivity to this threshold.

Despite the large variability in snow response between the models, snow processes are represented by a degree-hour method in all the models, suggesting a high sensitivity of the snow response to the snow process parameterization (Table ).

5.2.2 Root-zone soil moisture

Vegetation-accessible water volumes that can be stored in the root zone largely control the long-term partitioning of precipitation into evaporation and drainage. Most hydrological models include a representation of this root-zone storage capacity $S_{\mathrm{R},\max }$, which is estimated through calibration (Table ). The calibrated root-zone storage capacities vary between 74 and 277 mm across studied models. The root-zone soil moisture content relative to saturation of models with relatively large root-zone storage capacities (here defined as $S_{\mathrm{R},\max }>$ 200 mm) tends to never fully dry out ($>$ 0.20) and saturate ($<$ 0.94) as opposed to models with lower root-zone storage capacities ($S_{\mathrm{R},\max }<$ 150 mm), in which the storage tends to either dry out completely and/or fully saturate (Fig. a). If the vegetation-accessible water storage dries out, this will lead to water stress and reduced transpiration. On the other hand, if the root-zone storage is saturated, no more water can be stored, resulting in fast drainage. The size of the root-zone storage capacity is therefore a key control of the hydrological response, allowing us to explain some of the observed variability between the models.

Figure 6

(a) Range of relative root-zone soil moisture ${\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}}$ in the Ourthe upstream of Tabreux for the period 2015–2017 as a function of the median root-zone storage capacity ($S_{\mathrm{R},\max }$) across parameter sets. The feasible parameters for NAM are split into two groups due to the large variability of $S_{\mathrm{R},\max }$ (subsets with $S_{\mathrm{R},\max }$ of $\sim$ 130 and $\sim$ 240 mm). (b) Root-zone storage capacity $S_{\mathrm{R},\max }$ as a function of the optimal $T$ value for each model realization. Optimal $T$ values are derived at the highest Spearman rank correlation between the Soil Water Index and modeled root-zone soil moisture.

[Figure omitted. See PDF]

We compare the dynamics of modeled and remote-sensing estimates of root-zone soil moisture by calculating Spearman rank correlations between modeled root-zone soil moisture and remote-sensing estimates of the Soil Water Index for the available $T$ values of 2, 5, 15, 20, 40, 60 and 100 d. As the $T$ value increases, the Soil Water Index is more smoothed and delayed. For each model realization, we identify the $T$ value which yields the highest Spearman rank correlation between Soil Water Index and modeled root-zone soil moisture (Fig. b). The optimal $T$ value increases with the size of the calibrated root-zone storage capacity and varies between 15 and 60 d. A small root-zone storage capacity is indeed likely to fill through precipitation and empty through evaporation and drainage more rapidly than a large water storage capacity, leading to a higher temporal relative soil moisture variability. The mismatch between the relatively high root-zone storage capacities of VHM ($S_{\mathrm{R},\max }\sim$ 200 mm) in relation to the relatively low optimal $T$ values of 20 d is likely related to the unclosed water balance (Sect. ). The similarity between modeled root-zone soil moisture and Soil Water Index with optimal $T$ values is high, as implied by Spearman rank correlations varying between 0.88 and 0.90 across models. However, the disparity in optimal $T$ values between the models underlines the different temporal representations of root-zone soil moisture content across models, implying that all these models cannot simultaneously provide a plausible representation of the catchment-scale vegetation-accessible water content.

5.2.3 Total storage anomalies

Total water storage anomalies obtained from GRACE are compared to the storage as simulated by the models, showing relatively similar seasonal patterns, as illustrated in Fig. a for model M5. GRACE total storage anomalies of the Semois upstream of Membre-Pont and the Ourthe upstream of Tabreux are mainly represented by two neighboring cells (Fig. b), allowing us to test how models represent the observed spatial variability. The range of anomalies in the Semois upstream of Membre-Pont is larger than in the Ourthe upstream of Tabreux, implying 18 %, 3 % and 7 % less summer and 19 %, 19 % and 10 % more winter storage in the Semois upstream of Membre-Pont for each of the three GRACE processing centers (Fig. b). Median precipitation is also 37 % higher in the Semois upstream of Membre-Pont than in the Ourthe upstream of Tabreux during winter months (October–March) but relatively similar during summer months (April–September). This difference in precipitation potentially leads to a wider range of modeled anomalies in the Semois upstream of Membre-Pont than in the Ourthe upstream of Tabreux for all the models, as shown in Fig. c. This implies that all the models reproduce the spatial variability between both catchments observed by GRACE. As the models were calibrated for the Ourthe at Tabreux and parameter sets were transferred to the Semois upstream of Membre-Pont, the forcing data are the main difference for explaining the modeled spatial variability.

Figure 7

(a) Total storage anomalies modeled by M5 and compared to GRACE for the Ourthe upstream of Tabreux. The grey band shows the variability in total storage anomalies of the three processing centers. (b) Range of GRACE total storage anomalies for the three processing centers for the Semois upstream of Membre-Pont compared to the Ourthe upstream of Tabreux for the period 2001–2017. (c) Modeled total storage anomalies for both catchments. (d) Spearman rank correlations between GRACE and modeled total storage anomalies as a function of the range of modeled total storage for the Ourthe upstream of Tabreux.

[Figure omitted. See PDF]

The models are also able to represent the observed temporal dynamics of total storage anomalies, as suggested by Spearman rank correlation coefficients ranging between 0.62 and 0.80 for the Ourthe upstream of Tabreux (Fig. d). There is, however, no relation between the Spearman rank correlations of the anomalies and the total modeled storage range (difference between maximum and minimum values), as shown in Fig. d. PRESAGES, WALRUS, VHM and dS2 have the largest ranges of total modeled storage, varying between 260 and 280 mm, and are also characterized by relatively large root-zone storage capacities (PRESAGES and VHM) or no separate root zone (WALRUS and dS2), while the total storage range of all the other models is between 200 and 220 mm. The similarity in total storage range between most models is likely related to the identical forcing data and the similarity in the long-term partitioning of precipitation into drainage and evaporation (Sect. ). However, the absolute values of total storage during a specific event or the partitioning into internal storage components may vary between the models (Sect. ).

5.3 Interactions between storage and fluxes during dry periods

As previously seen in Fig. a, the relative root-zone soil moisture content of the GR4H model is always above 0.2 for the 3 years for which SCATSAR-SWI1km data are available, as opposed to M5, which fully dries out for some time during the summers of 2015–2017. The Normalized Difference Vegetation Index of MODIS NDVI, also does not show a sharp decrease during these periods (Fig. a, b). Actual evaporation in M5 is strongly reduced during these dry soil moisture periods, unlike GR4H, as shown in Fig. c, d. When zooming into the dry period around September 2016, Fig. e, f show median relative root-zone soil moisture in GR4H of $\sim$ 0.24 versus $\sim$ 0.01 for M5, while SCATSAR-SWI1km has a higher median value of $\sim$ 0.55 (for both optimal $T$ values of 20 and 40 d). The dryness of root-zone soil moisture in M5 leads to median daily evaporation of 0.8 against 1.3 mm d${}^{-\mathrm{1}}$ for GR4H and prolonged periods of almost-zero evaporation in M5 (e.g., 31 August –3 September, 9–15 September and 22–30 September), while this neither occurs in GR4H nor in GLEAM actual evaporation, as shown in Fig. g, h. Despite the high streamflow performance of model M5 (Fig. , Table ), it is unlikely that transpiration is reduced to almost zero for several days in a row each summer in a catchment where approximately half of the area is covered by forests. This is also not supported by the remote-sensing data of soil moisture, NDVI and evaporation. High streamflow performances, therefore, do not warrant the plausibility of internal process representation. Despite the dried-out root-zone storage in M5, there is still water available in the slow storage to sustain a baseflow close to observed values, as shown in Fig. j, l. The streamflow responses of GR4H and M5 are both close to observations (Fig. i, j) in spite of differences in storage and evaporation, suggesting different internal process representations for a similar aggregated streamflow response during a low-flow period.

Figure 8

(a,b) Modeled relative root-zone soil moisture ${\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}}$, SCATSAR-SWI1km Soil Water Index with optimal $T$ value and NDVI for the period 2015–2017 for GR4H (yellow) and M5 (orange), respectively. The error bars and bands show the standard deviation of the remote-sensing data within the catchment area. (c, d) Actual evaporation $E_{\mathrm{A}}$ by GR4H and M5 for the period 2015–2017, showing a large reduction of evaporation during summer for M5, unlike GR4H and GLEAM actual evaporation. (e, f) Zoomed-in modeled ${\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}}$ and SCATSAR-SWI1km root-zone soil moisture for the grey shaded period of September 2016 in (a), (b), (c), and (d). (g, h) Potential, modeled and GLEAM actual evaporation, (i, j) modeled and observed streamflow $Q$, and (k, l) total storage $S_{\mathrm{T}}$ for the September 2016 dry period. The narrow uncertainty band of the GR4H model is related to its converging parameter search method.

[Figure omitted. See PDF]

5.4 Plausibility of process representations

The models are ranked and evaluated for a selection of criteria using observed streamflow, remote-sensing data and expert knowledge (Figs.  and S3 in the Supplement). All the models deviate less than $\pm$6 % from observed median annual streamflow (Fig. a), which is less than the estimated uncertainty of 11 % (Sect. ). In contrast, the modeled flashiness and baseflow indices of most models deviate more than the estimated uncertainty (Fig. b, c). FLEX-Topo is the only model with a clear overestimation of the flashiness index, which relates to the calibration aim of having a flashy model to reproduce small summer peaks .

Figure 9

Ranking and evaluation of model behavior for a selection of criteria based on observed streamflow, remote-sensing data and expert knowledge. The grey shaded areas are soft indications of more plausible behavior based on uncertainty estimates and expert knowledge. Model ranks as a function of the (a) deviation from observed median annual streamflow; (b) deviation from the flashiness index; (c) deviation from the baseflow index; (d) deviation from median annual GLEAM actual evaporation; (e) deviation from median annual GLEAM interception for models with a separate interception module; (f) number of days with snow cover for valid MODIS observations between 2001 and 2017 for models with a snow module; (g) annual number of days when the root-zone storage is dry (filled with less than 1 % of its capacity); (h) deviation from the 1st percentile of GRACE total storage anomalies for the three centers. The error bars show the 25th–75th range across the ensemble of feasible parameter sets. Results are for the Ourthe catchment upstream of Tabreux, except for the snow analysis.

[Figure omitted. See PDF]

Modeled median annual total actual evaporation deviates by approximately $-\mathrm{10}$ % from GLEAM estimates, except for the $+$22 % overestimation of the VHM model due to the issue of the unclosed water balance, as shown in Fig. d. These results are consistent with the evaluation study of GLEAM compared to other evaporation products , which reports higher than average values for GLEAM in Europe ($\sim$ $+$10 % at our latitude).

Four models explicitly account for interception with a separate module. Median annual interception rates deviate substantially from GLEAM estimates ($-$47 % to $+$80 %) as shown in Fig. e. There is a high uncertainty in the partitioning of evaporation into different components in evaporation products, and GLEAM likely underestimates interception rates . Therefore, we consider a large uncertainty of $+$50 % to evaluate and rank the models. The GR4H interception is lower than GLEAM estimates. However, an interception storage was recently included in an hourly GR model (GR5H) to better represent the interception processes .

All the models substantially underestimate the number of days when snow is observed by MODIS at the catchment scale for all valid MODIS observations (cloud cover $<$ 40 % and excluding summer months), as shown in Fig. f. Yet we estimate a low uncertainty of less than 2 % around this number (Sect. ). The NAM, M4 and M5 models are closest to MODIS estimates, but they are characterized by high false alarm rates (Fig. d), which implies a mismatch in the modeled and observed days with snow for valid MODIS observations. Based on expert knowledge and the trade-off between miss rate and false discovery rate (Fig. d, e), we expect the annual number of days with snow storage to be between 28 and 62 d yr${}^{-\mathrm{1}}$ as modeled by wflow_hbv, WALRUS, FLEX-Topo and GR4H, whereas the $\sim$ 90 d yr${}^{-\mathrm{1}}$ of NAM, M4 and M5 seems too high.

The FLEX-Topo and M2 to M5 models are characterized by an empty root-zone storage for approximately 10 d yr${}^{-\mathrm{1}}$ (${\stackrel{\mathrm{‾}}{S}}_{\mathrm{R}}$ $<$ 1 %) as shown in Fig. g. These models have in common that evaporation from the root zone occurs at a potential rate and is not (or hardly) reduced when soils become dry until the point where the storage is empty. This is the case for models with a very low or no evaporation reduction parameter $L_{\mathrm{P}}$. This behavior is not supported by the remote-sensing data of evaporation, soil moisture and NDVI (Sect. ) nor by theory on root water uptake reduction under dry conditions . The additional slow groundwater reservoir added in model M5 compared to M2–M4 leads to a smaller root-zone storage capacity as the available storage is partitioned into the root-zone storage and the additional groundwater store. The smaller root-zone storage capacity of model M5 exacerbates the number of annual days with empty storage. This highlights the complex interactions in internal dynamics even in parsimonious lumped models with similar mean annual streamflow performance.

Catchments with relatively large root-zone storage capacities underestimate GRACE estimates of very dry storage anomalies most (Figs.  and h). The uncertainty of GRACE is represented by the estimates of the three processing centers and the $\sim$ 18 mm uncertainty estimate mentioned in Sect. . FLEX-Topo has a low root-zone storage capacity and is the only model which overestimates the very dry storage anomalies. Models with root-zone storage capacities of around 110 to 150 mm show the most consistent behavior with GRACE estimates of very dry storage anomalies.

6 Implications

While streamflow alone may be used to evaluate hydrological models, we subsequently use these models to understand internal states and fluxes in current and future conditions or to make operational streamflow predictions (e.g., HBV and GR types of models are used by the Dutch and French forecasting services). Our findings show that similar streamflow responses obtained by models calibrated according to an identical protocol rely on different internal process representations. While not unexpected, it implies that we might get the right answers but for the wrong reasons , as these models cannot at the same time all be right and different from each other .

Almost all the models show a similar long-term partitioning of precipitation into drainage and evaporation, as they are forced and constrained by the same data, also leading to relatively similar volumes of total storage. However, the partitioning of total storage into several internal storage components differs between the models, resulting in distinct runoff responses as expressed by the baseflow and flashiness indices.

None of the models is systematically consistent with the information available from streamflow observations, remote-sensing data and expert knowledge. However, some processes either play a limited role in the overall water balance or can be compensated by other processes. Snow occurs every year but is not a major component of the streamflow regime , interception evaporation can be compensated by root-zone evaporation, and very dry periods only occur for several weeks per year when streamflow is already very low. There is also a large uncertainty in each of the data sources, which makes us reluctant to use them to determine hard thresholds to reject models. Instead, we ranked the models for a selection of “soft” criteria and found that NAM, wflow_hbv and PRESAGES are overall most consistent with the evaluation data, with median ranks of 2–3 (Fig. S3 in the Supplement). While an overall ranking may be useful for practitioners, modelers benefit more from the specific ranking for each criterion to detect specific model deficiencies that could be improved in the model structure. An overall ranking is only a mere indication which should be interpreted carefully due to uncertainty in the evaluation data and the applied calibration strategy.

The presence of interception or a slow storage (absent in M2–M4 but added in M5) affects the representation of other internal processes, including transpiration and/or root-zone soil moisture, implying that individual internal model components are altered by the presence/absence of other potentially compensating processes. Adding an additional internal model component changes the internal representation of water storage and fluxes through the system, which should be kept in mind if model parameters were to be fixed in alternative model structures. Furthermore, model improvements through additional process components and/or adapted parameterization should be evaluated in terms not only of the aggregated response, but also the partitioning of fluxes and storages through the system (e.g., does the groundwater component improve the baseflow index at the expense of the availability of root-zone soil moisture during dry periods?). Models should be confronted with expert knowledge, e.g., on the occurrence of days with water stress or snow storage, to assess the plausibility of internal states and fluxes .

Applying these models to a future, more extreme climate in the same region might lead to contrasting insights regarding impacts of climate change, as also shown by studies of , and , in which model structures may lead to different signs of change in mean streamflow. Using one model or the other to assess the effect of rising temperatures on snow could lead to very different timescales of snow storage decline. Vegetation already experiences more intense water stress in some models compared to others, and this would be exacerbated in more extreme drought scenarios . More intense precipitation events could affect interception evaporation and therefore water availability in the root zone differently from one model to another. Beyond model structure, the experience each modeler has with their model and associated calibration procedure to constrain model parameters may also impact the simulation results .

Our findings should, therefore, encourage modelers to use multiple data sources for model calibration and evaluation, as already suggested by many other studies . Remote-sensing estimates of soil moisture, evaporation and total storage anomalies are available at the global scale, and in spite of potential biases with models, the temporal dynamics are useful for constraining our models . Additionally, it seems essential to support decision-makers by studies relying on multi-model and multi-parameter systems, as also suggested by and , to reveal uncertainties inherent to the heterogeneous hydrological world .

This study is the result of a joint research effort of scientists and practitioners gathering each year in Liège at the International Meuse Symposium to exchange interdisciplinary and intersectoral knowledge related to the Meuse basin. Although coordination of large international teams may be challenging, international studies favor close collaboration between scientists and practitioners, who can learn from each other to accelerate modeling advances . Another advantage of comparing modeling results of several research groups is to quickly detect small mistakes in the modeling process, including shifts in the time series or using forcing data of one catchment to model another catchment. While hydrograph characteristics were the main focus of the previous study , we gain distinct insights into the plausibility of model behavior by evaluating additional facets of internal process representation using remote-sensing data.

7 Knowledge gaps and limitations

Many aspects of the hydrological response remain unknown and can hardly be evaluated against observations. While in situ observations of snow, evaporation or soil moisture are rarely available at a sufficient spatio-temporal scale, remote-sensing estimates have the advantage of high spatial resolution and worldwide spatial coverage, though they often rely on models themselves and are affected by high and often unknown uncertainty. Comparing models with these independent observations is valuable for evaluating their consistency and detecting outliers. However, these observations cannot be considered representative of the truth as they rely on many assumptions themselves, hindering “real” hypothesis testing. The ratio of actual over potential evaporation as a result of water stress at the catchment scale, therefore, remains highly uncertain . While areal fractions of snow cover can be estimated by MODIS, the presence of clouds limits the usability of the data, and knowledge of the catchment-scale snow water equivalent is lacking. If remote-sensing estimates of near-surface relative soil moisture are available, root-zone water content remains uncertain, and while GRACE provides estimates of total storage anomalies, we lack knowledge on absolute total water storage. The spatial variability and the temporal dynamics of these remote-sensing products provide useful, additional, independent information to understand the hydrological puzzle but certainly not all the answers to evaluate the states typically included in process-based models. Measurements are, therefore, of crucial importance to increase our understanding of hydrological processes at the catchment scale, which in turn will improve the quality of remote-sensing products and model development .

The evaluation of model behavior is conditional on the calibration procedure, which was freely chosen by the individual contributing institutes. The use of different or more calibration objectives and in-depth uncertainty estimation may have resulted in different conclusions in terms of the plausibility of the behavior of each model.

We performed a thorough analysis of 12 models, five variables and three catchments. We deliberately chose to limit the number of study catchments to balance depth with breadth, allowing us to dive into process-relevant insights.

8 Conclusions

Similar streamflow performance of process-based models, calibrated following an identical protocol, relies on different internal process representations. Most models are relatively similar in terms of the long-term partitioning of precipitation into drainage and evaporation. However, the partitioning between transpiration and interception, snow processes and the representation of root-zone soil moisture varies significantly between the models, suggesting variability of water storage and release through the catchment. The comparison of modeled states and fluxes with remote-sensing estimates of evaporation, root-zone soil moisture and vegetation indices suggests that models with relatively small root-zone storage capacities and without reduction in root water uptake during dry conditions lead to unrealistic drying-out of the root-zone storage and significant reduction of evaporative fluxes each summer. Expert knowledge in combination with remote-sensing data further allows us to “softly” evaluate the plausibility of model behavior by ranking them for a set of criteria. Even if none of the models is systematically consistent with the available data, we did not formally reject specific models due to the uncertainty in the evaluation data and their changing relevance for the studied catchments. The dissimilarity in internal process representations between the models implies that they are not necessarily providing the right answers for the right reasons, as they cannot simultaneously be close to reality and different from each other. While the consequences for streamflow may be limited for the historical data, the differences may exacerbate for more extreme conditions or climate change scenarios. Considering the uncertainty of process representation behind the scenes of streamflow performance and our lack of knowledge and observations on these internal processes, we invite modelers to evaluate their models using multiple variables, encourage more experimental research, and highlight the value of multi-model multi-parameter studies to support decision-making.

Data availability

Streamflow and precipitation data were provided by the Service Public de Wallonie in Belgium (Direction générale opérationnelle de la Mobilité et des Voies hydrauliques, Département des Etudes et de l'Appui à la Gestion, Direction de la Gestion hydrologique intégrée (Bld du Nord 8-5000 Namur, Belgium)). Hourly radiation data were retrieved from the portal of the Royal Netherlands Meteorological Institute (http://www.knmi.nl/nederland-nu/klimatologie/uurgegevens, ). Daily temperature data were retrieved from the E-OBS (version 17.0) OPeNDAP server . Actual evaporation estimates from the Global Land Evaporation Amsterdam Model (GLEAM) are available through the SFTP server of GLEAM https://www.gleam.eu/,. MODIS snow cover fractions are available for download from the Earthdata portal at 10.5067/MODIS/MOD10A1.006 and 10.5067/MODIS/MYD10A1.006 . MODIS vegetation indices data are available for download at 10.5067/MODIS/MOD13A1.006 and 10.5067/MODIS/MYD13A1.006 . The Soil Water Index SCATSAR-SWI1km is available from the Copernicus Global Land Service at https://land.copernicus.eu/global/products/swi . GRACE land data are available at http://grace.jpl.nasa.gov , supported by the NASA MEaSUREs program. The modeled states and fluxes of each model are available online in the 4TU data repository at 10.4121/13221038.v1 .

The supplement related to this article is available online at: https://doi.org/10.5194/hess-25-1069-2021-supplement.

Author contributions

LJEB, GT, FF, MH, LAM, and AHW designed the study. GT, TdBE, JB, CCB, JdN, SM, BG, FF, JN and LJEB ran their models and provided simulated time series of all variables. LJEB conducted all the analyses and wrote the manuscript. All the authors discussed the design and results and contributed to the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Acknowledgements

We thank Deltares and Rijkswaterstaat for the financial support to conduct this analysis. The authors would like to thank the Service Public de Wallonie, Direction générale opérationnelle de la Mobilité et des Voies hydrauliques, Département des Etudes et de l'Appui à la Gestion, Direction de la Gestion hydrologique intégrée (Bld du Nord 8-5000 Namur, Belgium) for providing the precipitation and streamflow data. We are grateful for the E-OBS dataset from EU-FP6 project UERRA (http://www.uerra.eu, last access: 21 February 2021) and the Copernicus Climate Change Service and the data providers in the ECA&D project (https://www.ecad.eu, last access: 21 February 2021). We thank the editor Nadav Peleg, Keith Beven, Shervan Gharari, Uwe Ehret and three anonymous referees for their constructive comments which helped us improve the manuscript.

Review statement

This paper was edited by Nadav Peleg and reviewed by Keith Beven, Uwe Ehret, and three anonymous referees.