Effects of uncertainty in soil properties on simulated hydrological states and fluxes at different spatio-temporal scales

Soil properties show high heterogeneity at different spatial scales and their correct characterization remains a crucial challenge over large areas. The aim of the study is to quantify the impact of different types of uncertainties that arise from the unresolved soil spatial variability on simulated hydrological states and fluxes. Three perturbation methods are presented for the characterization of uncertainties in soil properties. The methods are applied on the soil map of the upper Neckar catchment (Germany), as an example. The uncertainties are propagated through the distributed mesoscale hydrological model (mHM) to assess the impact on the simulated states and fluxes. The model outputs are analysed by aggregating the results at different spatial and temporal scales. These results show that the impact of the different uncertainties introduced in the original soil map is equivalent when the simulated model outputs are analysed at the model grid resolution (i.e. 500 m). However, several differences are identified by aggregating states and fluxes at different spatial scales (by subcatchments of different sizes or coarsening the grid resolution). Streamflow is only sensitive to the perturbation of long spatial structures while distributed states and fluxes (e.g. soil moisture and groundwater recharge) are only sensitive to the local noise introduced to the original soil properties. A clear identification of the temporal and spatial scale for which finer-resolution soil information is (or is not) relevant is unlikely to be universal. However, the comparison of the impacts on the different hydrological components can be used to prioritize the model improvements in specific applications, either by collecting new measurements or by calibration and data assimilation approaches. In conclusion, the study underlines the importance of a correct characterization of uncertainty in soil properties. With that, soil maps with additional information regarding the unresolved soil spatial variability would provide strong support to hydrological modelling applications.


Introduction
The prediction of mathematical environmental models is affected by uncertainty, which arises from inadequate conceptual and mathematical representations of the processes (uncertainty in model structure), inadequate and insufficient knowledge and characterization of system forcing (uncertainty in boundary conditions) and limitations in the measurements or identification of model parameters (parameter uncertainty) (Beven, 2001(Beven, , 2007Refsgaard et al., 2007;Tartakovsky et al., 2012). The need to quantify the predictive uncertainty has led to the development of probabilistic (stochastic) frameworks in many disciplines of environmental sciences and engineering (Altarejos-García et al., 2012;Baroni and Tarantola, 2014;Di Baldassarre et al., 2010;Dubois and Guyonnet, 2011;Savage et al., 2016;Seiller and Anctil, 2014). Currently rigorous quantification of uncertainty is an integral part of science-based predictions and decision support systems (Beven, 2007;Farmer and Vogel, 2016;Liu and Gupta, 2007;Montanari and Koutsoyiannis, 2012).
In hydrological studies, several sources of uncertainty have been studied ranging from atmospheric forcing (Aguilar et al., 2010;Raleigh et al., 2015;Samain and Pauwels, 2013; G. Baroni et al.: Effects of uncertainty in soil properties Vázquez and Feyen, 2003;Zhu et al., 2013) to geology structures (Comunian et al., 2016;Hansen et al., 2014;He et al., 2015;Zech et al., 2015). Among these, the uncertainty related to the soil properties has been widely analysed. Soil properties show, in fact, high heterogeneity at different spatial scales with a hierarchy of spatial structures (Burrough, 1983;Heuvelink and Webster, 2001;Vogel and Roth, 2003) and complex interactions with environmental conditions (Lin, 2010). Despite international initiatives exist to improve the current status of soil characterization (Chaney et al., 2016;Heuvelink et al., 2016;Pelletier et al., 2016;Shangguan et al., 2014), detailed information of the spatial heterogeneity of the soil properties over large areas remains a crucial challenge. For this reason, an increasing number of hydrological modelling studies aim to integrate the uncertainty in soil properties that arise from the unresolved spatial heterogeneity for a proper quantification of the uncertainty of the model results. Since soil properties play a crucial role in the entire water cycle, this topic crosses research fields from lower atmosphere (De Lannoy et al., 2014;Garrigues et al., 2015;Guillod et al., 2013;Osborne et al., 2004;Yu et al., 2014) and surface water (Anderson et al., 2006;Geza and McCray, 2008;Li et al., 2013;Livneh et al., 2015;Salazar et al., 2008) to water and solute transport to groundwater systems (Besson et al., 2011;Hennings, 2002;Yu et al., 2014).
Despite its relevance, however, relatively simple assumptions are adopted to characterize the uncertainty in soil properties and to understand its effect on the hydrological response. In several studies the uncertainty is characterized based on a relatively small number of scenarios (Baroni et al., 2010;Christiaens and Feyen, 2001;Guber et al., 2009;Herbst et al., 2006;Hohenbrink and Lischeid, 2015;Islam et al., 2006;Mirus, 2015;Moeys et al., 2012) or by simple random noise (i.e. variance) added to the original soil properties (Arnone et al., 2016;Chaney et al., 2015;Deng et al., 2009;Garrigues et al., 2015;Han et al., 2014;Loosvelt et al., 2011). Other studies explicitly integrate the complex heterogeneity of the subsurface and the uncertainty in the soil properties is characterized based on spatial correlated random fields, i.e. specifying variance and correlation length (Binley et al., 1989;Fan et al., 2016;Fiori and Russo, 2007;Merz and Plate, 1997;Meyerhoff and Maxwell, 2011). Moreover, many of the above-mentioned studies focused on the effect of the uncertainty in soil properties on a selected hydrologic variable at specific temporal and spatial scales, e.g. rainfallrunoff events (e.g. Arnone et al., 2016;Fan et al., 2016), simulated evapotranspiration (e.g. Garrigues et al., 2015), soil moisture distributions (e.g. Liao et al., 2014) or groundwater recharge (e.g. Moeys et al., 2012). Simultaneous assessments of different hydrological components of the water balance at different spatial and temporal scales are rare. In addition, due to the different settings used in the studies, it is not possible to draw general conclusions about the role of the uncertainty in soil properties. In some cases the refined spatial information of soil properties does not contribute to a more accurate prediction (e.g. Li et al., 2013). In other studies the results showed to be very sensitive to the soil properties (e.g. Livneh et al., 2015). These controversial results foster the debate on the need (or not) for finer resolution soil maps in the different modelling applications (Baveye, 2002;Baveye and Laba, 2015;Heuvelink and Webster, 2001).
In the present study, we investigate impacts of uncertainty of soil properties on hydrological states and fluxes. Uncertainty in soil properties is characterized by three different methods that are consistent in the added noise (i.e. variance), but they differ in the perturbation of the soil spatial structure, i.e. correlation length. The first two methods were previously used in other studies (e.g. Fan et al., 2016;Han et al., 2014). The third method is developed in the present study to introduce small-scale soil variability while preserving the original spatial patterns. Therefore, we hypothesize that local responses of a hydrological system, such as evapotranspiration and soil moisture, will be strongly impacted by the uncertainty introduced at small spatial scale. However, integrated responses like the streamflow aggregate local responses over large areas. We hypothesize that this integrated response will be less impacted by soil properties uncertainty. The extent of the impact is expected to decrease with increasing the aggregation area and to disappear at a specific domain size. In such a condition, the system is stated to be spatially ergodic as the model output is not any more sensitive to the perturbation, i.e. we have the equivalence between spatial and ensemble statistics (Dagan, 1989;Rubin, 2003).
The paper is structured as follow. First, the perturbation methods used for the characterization of the uncertainty of the soil properties are presented. The specific case study is described presenting the catchment, the data used and the specific settings of the perturbation methods. The hydrological model is then introduced together with the uncertainty analysis conducted for the assessment of the effect of the uncertainty in soil properties on the simulated states and fluxes. The results are discussed in Sect. 3, focusing on the effect of the differences detected at different spatial and temporal scales. Final remarks are presented in the conclusions section.

Soil perturbation methods
In this section, the three statistical methods to characterize the uncertainty in soil properties are presented. A sketch for describing the methods is provided in Fig. 1, where one hypothetical horizontal transect through a soil map with three soil units characterized by different percentages of sand is shown as example.
The first method (hereafter denoted as random error method -RE) is based on the assumption that the nominal value in each soil unit is the only source of uncertainty Figure 1. Soil perturbation methods (random error method, RE; spatially correlated method, SC; and conditional points method, CP). The panel on the left shows the percentage of sand of a hypothetical horizontal transect trough the original soil map as an orange line. Within the transect three different soil units are observed, which leads to three different sand contents. Each row of the right panels depicts the steps for setting the perturbation methods. The blue line depicts one realization of the respective perturbation method. The detailed description of these methods can be found in Sect. 2.1. Abbreviations: var -variance, CL -correlation length. while the spatial patterns (i.e. soil units) are considered to be correct. To fulfil this assumption, a simple Gaussian random noise is defined with zero mean and given variance (Fig. 1, step R1). Random values are sampled from the distribution and added to the nominal value of soil properties of each soil unit ( Fig. 1, step R2). This approach was commonly used in several studies with the focus of understanding the effect of the soil properties in forward uncertainty analysis of model response (e.g. Deng et al., 2009) or for creating the forward ensemble in data assimilation tests (e.g. Han et al., 2014).
In the second method (hereafter denoted as spatially correlated method -SC), a similar assumption of additive random values is considered. However, it is also assumed that the uncertainty arises from the presence of smaller soil units that have not been detected in the original soil map (Hennings, 2002). To fulfil this assumption, a spatial structure (i.e. variance and correlation length -CL) is defined (Fig. 1, step S1). Based on that, a spatially correlated random field with zero mean is created (Fig. 1, step S2) and added to the original soil map (Fig. 1, step S3). Random fields are used in this approach to create variability as discussed by Goovaerts (2001) with which simulated short-range components well represent the complexity of the small-scale spatial structure. Readers interested in the details of the generation of random fields are referred to Deutsch and Journel (1998), Goovaerts (1997) and Isaaks and Srivastava (1989).
Finally, in the third approach (hereafter denoted as conditional points method -CP), it is assumed that the nominal value of the original soil units represents some point locations within this unit, but their positions are unknown. The uncertainty arises from the spatial variability within these point locations that is not resolved in the original soil map. To fulfil this assumption, points are randomly distributed over the soil map and the soil properties are associated with each position (Fig. 1, step S1). These values are used to calculate the spatial structure, i.e. the empirical variogram ( Fig. 1, step S2). A variogram model is fitted and a conditional random field is created using the sampled locations as conditional points (Fig. 1, step S3). It has to be noted that the CP method has some similarity with the pilot points approach used for the calibration of hydrogeological models (Carrera et al., 2005). The main difference is the use in this method of new points at each iteration; i.e. the points are located in different positions for each created conditional random field.
It is noteworthy that additional statistical methods for the analysis of soil map are presented in the literature (Goovaerts, 2011;Heuvelink et al., 2016;Kempen et al., 2009;Minasny and McBratney, 2016;Odgers et al., 2014). However, the aim of these methods is to downscale/disaggregate the information available in the original soil map and not to characterize its uncertainty. For this reason, these statistical methods are based on environmental covariates (i.e. environmental variables that co-vary with soil variability) known at higher resolution (i.e. digital elevation model or land use) and they require relative good knowledge of the soil formation and the specific settings to adopt (Kerry et al., 2012;Nauman and Thompson, 2014;Subburayalu et al., 2014;Du et al., 2015). On the contrary, the three methods selected and developed in the present study represent relative simple approaches only based on the information available in the original soil map. They can be applied for the characterization of any type of soil properties (texture, saturated hydraulic conductivity, soil depth etc.) and they reflect different assumptions regarding the uncertainties in the soil properties. For this reason, they can be tuned to characterize uncertainty for soil maps of any scales and they can be easily used with any modelling studies (e.g. sensitivity analysis or data assimilation). Combinations of the methods can also be considered when needed; i.e. soil maps affected by different types of uncertainties.

Study area
The numerical experiments are conducted in the upper Neckar catchment (Fig. 2) that was extensively investigated in previous hydrological studies (Kumar et al., 2010;Samaniego et al., 2010a, b;Wöhling et al., 2013b). This catchment is located in the central uplands of Germany and comprises a catchment area of approximately 4000 km 2 . The catchment has a gradient in elevation from 250 m to 1015 m a.s.l. with a mean elevation of 550 m. The catchment is prevalently characterized by cropped fields and forest but with a remarkably high degree of urbanization (11 %). The long-term mean annual precipitation is around 920 mm yr −1 .
Observed meteorological data, i.e. precipitation as well as minimum, maximum and average daily temperature, were provided by the German Meteorological Service (www. dwd.de/). These observations have been interpolated to a 4 km × 4 km forcing data set for the hydrological model using external drift Kriging. The potential evapotranspiration is estimated using the Hargreaves-Samani method (Hargreaves and Samani, 1985). Data characterizing the land surface are a digital elevation model (Federal Agency for Cartography and Geodesy), a soil map at the scale 1 : 1 000 000 (Federal Institute for Geosciences and Natural Resources -BGR), a hydrogeological map (Federal Institute for Geosciences and Natural Resources -BGR), and land cover information (CORINE, European Environmental Agency -EEA, 2009). The soil map used in the present study (BGR 1 : 1 000 000) contains soil texture (percentage of sand, clay and silt) and bulk density (g cm −3 ) for each soil unit (i.e. polygon of the soil map) and for each soil horizon. For this study, these vertical discretizations are not accounted for and the soil properties of each soil horizon are averaged to the total soil depth of 2 m (Fig. 3). The soil within the catchment is prevalently clay loam but with a relatively high spatial variability represented by 29 soil units (polygons) of different size within the catchment. All these data are discretized to a spatial resolu-tion of 100 m × 100 m. Readers interested in more details on data set and the processing may refer to Kumar et al. (2010), Samaniego et al. (2010b) and Zink et al. (2017). The spatial distributions of cumulative rain, potential evapotranspiration, land use and the mean annual leaf area index are shown in the Supplement (see Fig. S1).

Settings of the soil perturbation methods
In this section, the specific settings of each statistical perturbation method used for the characterization of the soil properties are described. The three methods are used independently to generate three different ensembles to identify the impact of the different uncertainties introduced in the original soil map on simulated states and fluxes.
Considering the random error method (see Fig. 1), a Gaussian random additive noise is used with standard deviation 7 % and 0.07 g cm −3 for soil texture (sand and clay) and bulk density, respectively (Table 1). A correlated sampling design is used to preserve the correlation between the original soil properties (e.g. negative correlation between sand and clay). These noises are selected to perturb the soil properties within the original soil class; i.e. it is assumed that the exact values of the soil properties are unknown but the soil class (e.g. clay loam) is correct. Similar ranges were also applied in other studies (Han et al., 2014;Hennings, 2002).
For the spatially correlated method (see Fig. 1), the parameters for the variogram and co-variogram models are selected to be consistent with the perturbation introduced in the random error method (Table 1). In particular, exponential variogram models are prescribed with the same effective noises used in the random error method (i.e. standard deviation 7 % and 0.07 g cm −3 for texture and bulk density, respectively) and preserving the correlation between the original soil properties. The correlation length of 3 km is selected to represent relative small spatial patterns that were not captured by the original soil map, i.e. patterns smaller than most of the soil units (Fig. 3). The variogram and co-variograms models selected are shown as Supplement (see Fig. S2).
Finally, considering the conditional points method (see Fig. 1), tests are conducted to identify the density of the conditional points within the soil map. One sample at every 3 km × 3 km is found to introduce the same noise prescribed by the other two methods (Table 1). A stratified spatial random sample is used to distribute the points within each soil unit. Based on these, two nested exponential variogram and co-variogram models are fitted to the experimental variograms based on ordinary least-squares residuals (Pebesma, 2004). These variogram models are used to create the conditional random fields. The experimental variograms and the fitted models for one realization (i.e. one random field) are shown, exemplarily, in the Supplement (Fig. S3).
In each method, the perturbed values are forced to a realistic range, i.e. texture values between 0 and 100 % and the sum of textural fractions equal 100 %. Therefore, it has  to be noted that these constrains (i) could modify the Gaussian noise introduced and (ii) could lower the uncertainty in areas of the basin where the actual values are close to the bounds. These constrains did not affect the spatial patterns of the generated soil maps of the present study due to the relative small perturbation introduced and the presence of limited areas with extreme texture conditions. However, atten-tion has to be paid in cases where these features are more relevant.
For each method, an ensemble of 100 realizations is created to characterize the uncertainty in soil properties. The analysis is conducted with the statistical software R 3.2.x (R Core Team, 2013) using add-on packages (Pebesma, 2004). The multi-variate conditional random fields were gen- Table 1. Parameter settings for each perturbation method (random error, spatially correlated and conditional points). Variogram models used for the spatially correlated and conditional points methods are showed in the Supplement (Figs. S2 and S3, respectively).

Perturbation method Parameters
Specific settings

Random Error
Standard deviation 7 % and 0.07 g cm −3 for texture and bulk density, respectively Spatially correlated Variograms and co-variograms models Exponential models (see Fig. S2) Effective variance 50 % 2 and 0.05 g 2 cm −6 for texture and bulk density, respectively. These values are equivalent to the noise (i.e. standard deviation) introduced with the random error method.
Correlation length 3 km

Conditional points
Density of samples 1 sample every 3 km × 3 km

The hydrological model mHM
The effect of the uncertainty in soil properties as characterized by the three perturbation methods on hydrological states and fluxes is analysed using the mesoscale hydrological model (mHM). The mHM (Kumar et al., 2013;Samaniego et al., 2010b) is an open-source, spatially distributed hydrologic model (www.ufz.de/mhm). It considers interception, snow accumulation and melting, soil water retention, evapotranspiration, percolation and runoff generation as main hydrologic processes. The multiscale parameter regionalization (MPR) method embedded in mHM allows for the application of the model at various locations and scales (Kumar et al., 2013;Rakovec et al., 2016). MPR accounts for sub-grid variabilities by estimating model parameters at the scale of the morphological input, e.g. 100 m × 100 m. Subsequently, these parameters are upscaled to the model resolution based on different average rules (harmonic mean, arithmetic mean etc.). For a detailed model description and the MPR scheme, interested readers may refer to Samaniego et al. (2010b) and Kumar et al. (2013). For this study, the soil within mHM is discretized into three layers, the first layer is 5 cm, the second layer is 25 cm and the third has a variable thickness. The depth of latter is based on the information provided by the soil map (2 m). Based on the soil textural properties, mHM estimates effective parameters for porosity, hydraulic conductivity, field capacity and permanent wilting point using a set of pedotransfer functions (Cosby et al., 1984;Twarakavi et al., 2009;van Genuchten, 1980;Zacharias and Wessolek, 2007). The list of the functions is reported in Table S1 in the Supplement.
The model was calibrated and validated in previous studies showing very good capability to match streamflow measurements at catchment of different sizes (Kumar et al., 2010(Kumar et al., , 2013Samaniego et al., 2010b;Wöhling et al., 2013b). The same parameterization is used for the present study. We establish the mHM over the upper Neckar catchment at 500 m spatial resolution resulting in 16 432 grid cells. The model run is conducted at an hourly timescale. All simulations are conducted with a 5-year model spin-up time (1985)(1986)(1987)(1988)(1989) to minimize the effect of inappropriate initial conditions. The implications of uncertain soil properties are evaluated showing the uncertainty in simulated routed streamflow (SF), generated runoff at every grid cell (Q), actual evapotranspiration (AET), volumetric soil moisture (SM) in the upper 30 cm and groundwater recharge (GWR). For each perturbation method 100 simulations were performed that yield a total of 300 simulations. The results obtained during 1 year of forward simulation (1990) are shown, as an example. This year is selected to represent average climate condition of the area (i.e. two rain seasons concentrated in spring and fall and a relatively dry summer season) but with a relatively high variability within the catchment (see Fig. S1).

Uncertainty analysis at different spatio-temporal scales
The uncertainty in simulated states and fluxes is quantified based on the coefficient of variation (CV %) to allow for comparability between the results obtained in the different model outputs. Assuming a generic variable v representing Sect. 3.5 simulated state or fluxes, CV is calculated as follows: where σ is the standard deviation of the variable v at each cell i and time t calculated based on each perturbation method m (i.e. random error method, spatially correlated method or conditional points method) as follows: where N ens is the number of ensemble members (i.e. 100), j one single ensemble member and µ represents the mean of the ensemble at each cell i and time t calculated as follows: The values obtained with the three perturbation methods are compared by aggregating the simulated states and fluxes at different spatial and temporal resolutions. In particular, four analyses are conducted (Table 2). In analysis no. 1, the spatial variability of the uncertainty of the simulated states and fluxes is presented, i.e. depending on the geographical location within the catchment. In this case the average CV calculated for the entire simulation period (i.e. 1 year) in each grid cell is quantified as follows: where T is the number of simulations time steps (i.e. 365 days). This value is used to represent and discuss the av-erage uncertainty obtained in the specific cell i and its spatial variability within the catchment. In analysis no. 2 (Table 2), the daily temporal dynamic of the uncertainty obtained at each grid cell is discussed. For this reason the CV m i,t calculated at the daily time step (Eq. 1) is directly compared for two representative grid cells selected within the catchment (see Fig. 2, point A and B).
The uncertainty on simulated states and fluxes is further compared by aggregating the model outputs at different resolution to identify the effect of the spatial scale on the performance of the model as discussed by Refsgaard et al. (2016). In particular, for use in analysis no. 3 (Table 2), subcatchments of different sizes are defined based on 36 gauging stations located within the catchment (see Fig. 2). The effect of the uncertainty in soil properties to the streamflow routed to the outlet (SF) of each subcatchment is then compared. For the other simulated model outputs (i.e. evapotranspiration, soil moisture and groundwater recharge), the values of each grid cell within the subcatchment are aggregated calculating the average of simulated model output v obtained at the finer resolution as follows: v m,j where N sc is the number of grid cells within the subcatchment sc. The value v m,j sc,t is used in Eqs.
(1)-(4) to calculate and compare the coefficient of variation of the mean simulated states and fluxes for the subcatchments of different sizes.
Finally, in analysis no. 4 (Table 2), the effect of the aggregation of states and fluxes at different resolutions is further analysed based on the approach shown by Hansen et al. (2014) and Rasmussen et al. (2012). In this case, the generic model output v is averaged coarsening the model grid at different resolutions r d (i.e. r 2 = 2 km, r 4 = 4 km, r 8 = 8 km, r 16 = 16 km, r 32 = 32 km). These values are substituted in Eqs. (1)-(3) to calculate the coefficient of variation in each new coarsened grid cell i. In this analysis the average of the CV across the entire domain and over the entire simulation period (i.e. 365 days) is calculated as a summary statistic as follows: where N r represents the number of cell i within the coarsened domain r d .
In addition to the spatial dimension, in this study, the same procedure is also repeated for each spatial aggregation r d considering a time aggregation t d . In particular, all the simulated model outputs v obtained at daily time step are averaged at t 10 = 10 days, t 30 = 30 days, t 60 = 60 days, t 120 = 120 days and t 180 = 180 days. These values are substituted in Eqs. (1)-(4) to calculate the coefficient of variation in each temporal aggregation t d and they are considered to represent the uncertainty in the simulated states and fluxes in case the averaged values are used in the assessment of the performance of the model.
The four analyses described above are conducted based on the results of 100 simulations obtained with the distributed hydrological model for each perturbation methods. A total of 300 simulations, analysed in 12 cases, are discussed in the results section (Table 2).
3 Results and discussion 3.1 Perturbation of the soil properties Three methods are used to perturb the values of the original soil map, i.e. sand (%), clay (%) and bulk density (g cm −3 ). This section discusses the results obtained on clay percentage exemplarily. Similar results are obtained for the other soil properties and for the related soil hydraulic parameters (see Supplement, Figs. S4-S6). For each method, Fig. 4 (top row) shows one realization of the perturbed clay percentage. In addition (Fig. 4, bottom row), one transect along the catchment is selected and the clay percentage of the original soil map, of one realization and of the ensemble spread (95 % confidence interval) are shown. The longitudinal transect was selected to capture the strong variability in the soil units detected along this direction (see Fig. 3).
The random error (RE) method preserves the shapes of the soil units and perturbs just the nominal values. The results therefore show how the contrasts between the soil units are modified and in some cases are exaggerated. For this reason, it is noteworthy to observe that this method could create nonrealistic spatial patterns since soil properties usually show smother changes in space. The results obtained based on the spatially correlated method (SC) show that the shapes of the soil units are still highly identifiable and the sharp changes between the units are still preserved. With this method, however, the random fields superimposed on the original soil map were selected with a correlation length of 3 km (see Sect. 2.3). For this reason, smaller spatial structures than the original soil units are introduced and the sharp changes in the soil properties are not uniformly distributed all over the soil unit. Finally, considering the results obtained with the conditional point (CP) method, the results show that the soil units are visible but the contrasts are completely smoothed eliminating the artefact of the original soil map. However, the spread (grey area) in this transition between soil units (polygons) is wider than the spread detected within each soil unit. The effect is due to the combination of the uncertainty introduced to the nominal value of the soil property and to the exact position of the transition between the soil units.
The spread of the realizations is quantitatively evaluated based on the standard deviation of the ensemble. In particular, Fig. 5a represents the probability distribution of the standard deviation of the clay percentage calculated at ev-  ery grid cell within the catchment (i.e. 16 432 grid cells) for each method. Results obtained based on the three methods, on average, exhibit a high consistency in representing the uncertainty over the catchment (i.e. average standard deviation is for all the methods 7 %). However, some differences are detected in the distributions. The RE method shows a normal distribution with a relatively low variability (i.e. the coefficient of variation of the distribution is 6 %). This is the consequence of the fact that the soil properties within the catchment are perturbed with almost the same magnitude. Similarly, the SC method also shows a normal distribution but with a slightly wider variability (i.e. CV = 8 %). In contrast, the results obtained with the CP method show a very different distribution that is skewed with a tail to high extreme values. These high spreads in the soil realizations are located in the transition between the soil units, in particular if the transition is sharp (see Fig. 4).
Finally, the standard deviation of the values is calculated by aggregating the map for different subcatchments (Fig. 5b) and at different grid resolution (Fig. 5c) based on the analysis described in Sect. 2.5. The spreads of the realizations obtained with the three perturbation methods are of similar magnitude considering the finer resolution (e.g. resolution < 1 km × 1 km). The differences between the perturbation methods become more relevant by aggregating to a coarser resolution. In other words, this means that the introduced uncertainty is in the same order of magnitude, but the spatial patterns are different. In the RE method, the spread is relatively high at all scales and even the mean value of the soil properties of the entire catchment is perturbed (i.e. standard deviation > 0 for the resolution of 60 km × 60 km). The spread obtained with the SC method decreases more rapidly with increasing spatial scale. This is consistent with the correlation length prescribed to the random fields used in this method (i.e. 3 km). However, it is noteworthy to observe how the mean value of the soil properties over the entire catchment is still perturbed also with this method. This is explained by the fact that the random fields superimposed to the original soil map have zero mean over a rectangular domain, but the average can be different when masked to the catchment. The behaviour is exaggerated when a relatively long correlation length in comparison to the size of the domain is used. Finally, the results of the CP method show how just the small scale is perturbed and the spread of the ensemble drops already at the resolution of 5 km to disappear completely when the average over the catchment is considered. This behaviour is consistent with the density of the samples used to constrain the random fields (i.e. one sample every 3 km × 3 km).

Spatial variability of the uncertainty of states and fluxes
In this section, the spatial variability of the uncertainty of the simulated states and fluxes is presented. In this analysis (see Sect. 2.5, Table 2, uncertainty analysis no. 1), the mean coefficient of variation over time (i.e. 1 year) is calculated for each grid cell (i.e. 16 432 grid cells) and the spatial distributions obtained with the three perturbation methods are compared (Fig. 6).
The uncertainties of all hydrological states and fluxes obtained with each perturbation method provide nearly the same magnitude and the same spatial variability, with correlation coefficients calculated between the results obtained by each method higher than 0.8. For this reason, only the spatial distribution of the CVs of the model outputs over the entire catchment obtained with the RE method is shown, as an example (Fig. 6, left). The results obtained with all the three perturbation methods are shown for the transect depicted in Fig. 4 to facilitate the visualization of the relatively small differences (Fig. 6, right).
In general, the results obtained show that, independently from the perturbation method used, the uncertainty in the total runoff (Q) and groundwater recharge (GWR) are highest, with an average CV estimated across the catchment of 11 and 15 %, respectively. Soil moisture (SM) and actual evapotranspiration (AET) appear to be less sensitive to the soil variability with an average CV of 3 and 1 %, respectively (Fig. 6, left). The relatively small differences detected based on the use of different perturbation methods are located in the transition between the soil units (Fig. 6, right) and they are attributed to the higher uncertainty in the soil properties introduced in those areas (see Fig. 5a). Overall, a strong spatial variability in the uncertainty in the model outputs is detected with some differences depending on the considered model output. The uncertainty in runoff is more pronounced in the north-west areas, whereas actual evapotranspiration appears to be more affected in the central-north areas. High uncertainty in simulated soil moisture is distributed across Figure 6. Spatial variability of the uncertainty (CV) in the simulated model outputs (Q is generated runoff; AET is evapotranspiration; SM is soil moisture; GWR is groundwater recharge). In the left column, the results obtained based on the random error method (RE) over the entire catchment are depicted together with the position of the transect (dashed black line) and the two grid cells (blue points). In the right column, the CVs along the transect within the catchment based on the three perturbation methods (random error -RE; spatially correlated -SC; conditional points -CP) are plotted. Vertical dashed grey lines indicate the position of the grid cells A and B within the transect. Please note that all the plots have individual limits for the y axis. the catchment and the uncertainty in simulated groundwater recharge increases close to the catchment outlet (see Fig. 2).
For a further interpretation, the spatial variability of the uncertainty in the simulated model outputs is compared to different boundary conditions and input properties. In particular, the correlation coefficients between the spatial distribution of the CVs of each model outputs and the spatial distribution of clay (%), the mean leaf area index -LAI (m 3 m −3 ) and the annual sum of the potential evapotranspiration -PET (mm) calculated over the simulation period (i.e. 1 year) are calculated. These three factors are selected to represent soil, vegetation and atmospheric conditions. The spatial distributions of these factors are shown in the Supplement (see Fig. S1). Correlation coefficients for each perturbation method are calculated and average and standard deviation based on the three methods are depicted in Fig. 7.  Figure 7. Correlation coefficient calculated between the spatial distributions of the uncertainty (CV) of the simulated model outputs (Q is runoff; AET is actual evapotranspiration; SM is soil moisture; GWR is groundwater recharge) and local environmental conditions (the clay (%) is used to represent the soil; annual mean leaf area index LAI (m 2 m −2 ) is used to represent the vegetation; cumulative potential evapotranspiration PET (mm yr −1 ) is used to represent the atmospheric water demand). The bars represent the mean of the correlation coefficients obtained with the three perturbation methods and the error bars the standard deviation.
The results obtained with the three different methods are consistent between each other also in this comparison (i.e. as represented by the small error bars) showing different correlations for each model output. The uncertainty in the runoff is stronger correlated to the actual value of the soil property. This correlation can be visually identified comparing the spatial variability detected in Fig. 6 (right) and the spatial variability of the soil property shown in Fig. 4 for the same transect. The uncertainty in the actual evapotranspiration is strongly correlated to the atmospheric conditions and, to less extend, to the soil properties. Finally, the uncertainties of the soil moisture and groundwater recharge are correlated to the vegetation characteristics, with a relatively lower effect of soil properties.
To further evaluate the different correlations found for each simulated model output, the correlation matrix between the uncertainty (CV) detected in each model output is calculated (Table 3). On the one hand, the results show that the uncertainties in the fluxes are positively correlated (correlation coefficient > 0.2). This means that when the uncertainty in one specific flux is relatively high, also other fluxes to some degree are uncertain. On the other hand, it is interesting to note that the uncertainty in soil moisture is highly correlated to the groundwater recharge (correlation coefficient = 0.7) while the correlations to the other model outputs are negligible (correlation coefficient < 0.2). This means that the model could have relatively low uncertainty in soil moisture but high uncertainty in evapotranspiration or runoff and vice versa. These results are consistent among all the three Table 3. Correlation matrix of the uncertainty (CV) of the model outputs (Q is generated runoff; AET is actual evapotranspiration; SM is soil moisture; GWR is groundwater recharge) obtained with the three perturbation methods (random error, spatially correlated and conditional points). perturbation methods and they support the use of both states and fluxes for a proper assessment of the performance of hydrological models as it was underlined in several other studies (Ahmadi et al., 2014;Baroni et al., 2010;Conradt et al., 2013;Delsman et al., 2016;McCabe et al., 2005;Rakovec et al., 2016;Silvestro et al., 2015;Wöhling et al., 2013a;Zink et al., 2017).

Temporal variability of the uncertainty of states and fluxes
The daily temporal variability of the uncertainty on the simulated states and fluxes obtained at the model resolution (i.e. 500 m) is presented in this section. In this analysis (see Sect. 2.5, Table 2, analysis no. 2), the coefficient of variation at daily time step for each perturbation method obtained in two grid cells selected within the catchment are compared for an illustrative purpose. The two locations A and B are depicted in Fig. 2. The two grid cells are characterized by (see also Supplement, Fig. S1) a remarkable difference in the precipitation (i.e. almost 1600 and 1000 mm yr −1 , respectively), and by different land use (i.e. crop field and deciduous forest, respectively) but they have almost the same soil properties (i.e. 19 % sand and 59 % clay for grid cell A; 19 % sand and 66 % clay for grid cell B). The grid cells are selected to represent different uncertainties of the model outputs (see Fig. 6).
In particular, grid cell A shows relatively high uncertainty in simulated soil moisture (CV ∼ 4 %) while grid cell B shows relatively low uncertainty (CV ∼ 2 %). The three perturbation methods provide nearly the same results with a correlation coefficient higher than 0.8. For this reason, only the results obtained with the RE method are shown in Fig. 8. The figure also shows the actual values of simulated states and fluxes for comparison (i.e. mean value and 95 % confidence interval of the ensemble simulations obtained with the random error method). The results show how the uncertainty of the total runoff is relatively high during the entire simulation period with a tendency of increasing the uncertainty during a high flow period. The behaviour is particularly evident in the grid cell B (i.e. correlation coefficient between CV and simulated runoff is 0.6). In contrast, the actual evapotranspiration is close to the potential rate for most of the simulation period and, for this reason, it is not sensitive to changes in soil properties. As expected, the uncertainty is only detected during summer time when soil moisture is relatively low and the actual evapotranspiration rate decreases in comparison to the potential evapotranspiration. This result also explains the low correlation detected between the uncertainty in soil moisture and evapotranspiration (see Table 3). The temporal variability obtained for the uncertainty in soil moisture shows a more complex behaviour depending on the grid cell considered. In grid cell A, the CV increases with the increasing of soil moisture while it decreases in grid cell B. The different behaviours are explained comparing the actual soil moisture values. In grid cell A, the soil moisture values are relatively low (0.25 m 3 m −3 ) while, in grid cell B, the values are close to saturation (0.4 m 3 m −3 ). Finally, groundwa-ter recharge also shows a strong temporal dynamic with a tendency of higher uncertainty with increasing groundwater recharge in grid cell A (correlation coefficient = 0.2) while the correlation is negligible in grid cell B (correlation ∼ 0).
Overall, it is noteworthy to observe how the uncertainty in soil moisture is relatively constant in time while the uncertainty in the fluxes shows much stronger temporal variability. This different behaviour can be explained considering two main characteristics. On the one hand, the presence of nonlinear relations between states and fluxes generates threshold behaviour for which the uncertainty in soil moisture could be limited to ranges where the fluxes are not affected. This is for instance the case when the uncertainty in soil moisture is limited to relative wet conditions (i.e. above plant stress) and for this reason it does not affect the evapotranspiration. Similarly, the uncertainty in soil moisture could be limited in relatively dry conditions and the runoff could be not affected. On the other hand, there is a tendency of compensation in the uncertainty in the model outputs for which an overestimation of the actual evapotranspiration could be related to an underestimation of the groundwater recharge (or vice versa). In these conditions the soil moisture could be still well defined without providing any indication of the degradation of the model performance. As a result, the low uncertainty in soil moisture does not represent the overall uncertainty in the model. Overall, this analysis underlines the role of the different hydrological conditions (e.g. dry or wet) for understanding the effect of the uncertainty in soil properties on the model response. Similar conclusions are supported by the use of temporal sensitivity and identifiability analysis to better capture the role of the different uncertainties in the parameters analysed (Ghasemizade et al., 2017;Guse et al., 2016;Pianosi and Wagener, 2016;Wagener et al., 2003).

Spatial uncertainty of states and fluxes at subcatchments
The uncertainties (CV) of simulated states and fluxes are also compared by aggregating the results over subcatchments of different sizes (see Sect. 2.5, Table 2, analysis no. 3). The results obtained with the three perturbation methods are shown against the catchment size in Fig. 9. As presented by Refsgaard et al. (2016), the uncertainty in all the model outputs reduces with increasing catchment area. Assuming an arbitrary threshold (i.e. CV) acceptable for a specific model application, this analysis identifies, on the one hand, the spatial limit of model predictive capability for the specific application. On the other hand, it identifies the resolution above which it might become important to have a better understanding of the soil spatial variability. This resolution is referred to as the Representative Elementary Scale (RES) by Refsgaard et al. (2014) and it provides a clear and simple framework for the assessment of the performance of distributed models. However, it is interesting to note that the three perturbation methods generated very dif- Figure 9. Uncertainty, i.e. coefficient of variation (CV), of hydrological states and fluxes at catchments with different sizes (SF is streamflow, AET is evapotranspiration, SM is soil moisture, GWR is groundwater recharge). Exponential curves are fitted to the data. Please note that all figures have individual limits for the y axis. ferent results and, assuming the same arbitrary threshold for each method, different RESs are identified. The RE method creates higher uncertainty in all the subcatchments and even the mean of states and fluxes over the entire catchment is uncertain (i.e. 60 km × 60 km). The SC method shows a similar pattern but the uncertainty is lower in all the subcatchments. Finally, the uncertainty based on the CP method decreases already at small catchment sizes of e.g. 2 km × 2 km.
The different results in the uncertainty in the model outputs obtained by the use of the different perturbation methods are consistent with the different uncertainties introduced in the soil properties (Fig. 5b). This result supports the conclusion that the different RESs identified are related to the underlying correlation length (CL) scale used in each perturbation method (Refsgaard et al., 2016). The RE method perturbs the value of the entire soil units and it does not generate spatially ergodic soil parameters fields, i.e. aggregated hydrological responses still show a non-vanishing uncertainty at large catchment. The SC method introduced correlation length of 3 km and the effect on the uncertainty in the aggregated model output reached a remarkable reduction (e.g. > 90 % of the uncertainty in all the simulated states and fluxes is reduced) when the entire catchment is considered. Finally, the CP method introduces uncertainty only at small spatial scales while the longer spatial patterns are preserved. For this reason the domain is ergodic already at relatively low catchment size (i.e. 20 km × 20 km).
These characteristic lengths (RES vs. CL) identified by the use of the three different soil perturbation methods are in agreement with previous studies conducted in surface hydrology (Binley et al., 1989;Fan et al., 2016;Herbst et al., 2006;Merz and Plate, 1997) and in stochastic subsurface hy-drology (Dagan, 1989;Fiori and Russo, 2007;Rubin, 2003), where a suitable value for defining ergodic system was found to be ∼ RES/CL > 20. For this reason, it is notable the equivalence of the ergodic concept introduced in subsurface hydrology (Dagan, 1989;Rubin, 2003) and the RES concept in the case the arbitrary threshold (i.e. CV) is set to zero. However, two important characteristics can be further underlined. First, it is notable how catchments with similar size provide different degrees of uncertainties in the model output. This behaviour is in agreement with the results discussed in Sect. 3.2 showing different sensitivity on the soil perturbation depending on the different boundary conditions and model set-up (i.e. depending on the location within the catchment). For this reason, the results support the difficulties to find a universal RES that is not affected by the uncertainty in the soil properties for the entire catchment. Despite the RES concept has some differences with the Representative Elementary Area (REA) concept introduced in past literature (see Refsgaard et al., 2016, for further discussion about the differences), it is noteworthy how this result is in agreement with the difficulties for finding a universal REA discussed also in those studies (e.g. Fan and Bras, 1995;Wood et al., 1988). Second, different sensitivities arise depending on the model output considered. Soil moisture is more sensitive to the perturbation of soil properties since the relative change between the three different methods is the highest among the four hydrological variables under investigation. This behaviour is particularly evident when considering the results obtained with the random error method. In this case, a relatively small perturbation introduced in the mean of the entire catchment (60 km × 60 km) explains already most of the uncertainty in the simulated soil moisture. The uncertainty slightly increases with decreasing catchment size. In comparison, all the fluxes are much less affected by the small perturbations introduced for the entire catchment but they become increasingly pronounced with a decreasing catchment size. For this reason, the RES is also different depending on the model output considered.

Uncertainty of states and fluxes at different spatio-temporal scales
A similar scaling analysis is also conducted averaging states and fluxes by coarsening the grid resolution and by aggregating at different temporal scales (see Sect. 2.5, Table 2, analysis no. 4). The results obtained with the three perturbation methods are presented in Fig. 10. The spatial aggregation of the model output, as represented in the x axis of Fig. 10, shows the same effect obtained by aggregating the model output based on catchment of different sizes (Fig. 9). For this reason, the two analyses (aggregating by catchment vs. coarsening the grid resolution) can be considered equivalent in the identification of the effect of the spatial resolution on the uncertainty in the model outputs. However, the results described in the previous sections Figure 10. Spatio-temporal uncertainty analysis by aggregating the model results at different spatial and temporal resolutions. The three columns refer to the results obtained by (left) random error method -RE, (middle) spatially correlated method -SC and (right) conditional points method -CP. The rows refers to the different model outputs (i.e. Q is runoff, AET is actual evapotranspiration, SM is soil moisture; GWR is groundwater recharge). Note that a smooth approximation is depicted to facilitate the visualization of the actual CVs values.
showed also a strong variability in space and in time. For this reason, the use of the mean coefficient of variation calculated over time and across all the grid cells to represent the model performance (Eq. 6) can be misleading, e.g. underestimating the actual uncertainty in the model output. Instead, the use of the maximum CV calculated over the catchment and over the simulated period could be used to better represent the model performance. In addition, the extension of the analysis to the temporal scale (y axis in Fig. 10) emphasizes the clear tradeoff of the performance of the model between the spatial resolution and the temporal resolution. In particular, assuming an arbitrary threshold (i.e. CV) as a limit of model predictive capability for the specific model application (Refsgaard et al., 2016), the spatial and temporal analysis shows how the simulated states and fluxes should be aggregated in time to maintain acceptable performance when increasing the space resolution.
Overall, also for this spatio-temporal analysis, the results obtained with the three perturbation methods are very different and the RES (here defined as the spatial and tempo-ral scale at which it might become important -or not -to have a better understanding of the soil spatial variability) strongly depends on the perturbation methods used. Since the three perturbation methods reflect different uncertainties introduced in the original soil map, the analysis emphasizes the importance of identifying the correct approach to characterize the uncertainty for each model application and for further model improvements. For the specific case study presented here, it is notable how the streamflow at the catchment outlet (i.e. spatial resolution > 32 km), which was used for calibration of the model in previous studies (Kumar et al., 2013;Samaniego et al., 2010b), is sensitive only to the perturbation of long soil spatial structures introduced with the random error method. For this reason, it could be assumed that the uncertainty in soil properties introduced with the RE method is well compensated by the calibration and the RE method could not represent the actual uncertainty in the specific model application. The same could be considered for the results obtained with the spatially correlated method, as soon as subcatchments of different sizes are used in the calibration. In contrast, this model output is not sensitive to the perturbations introduced at small scale (e.g. conditional points method). On the one hand, this means that small soil variabilities are not relevant when the model application focuses on the streamflow prediction. On the other hand, this results underlines that it is not possible to infer (e.g. calibrate) these small spatial soil patterns based on the streamflow observations. For this reason, the conditional points method appears to be a simple and effective method to preserve the general spatial pattern of the original soil map while introducing uncertainty due to the unresolved spatial heterogeneity within the soil units. This type of uncertainty affects the streamflow only for small subcatchments (size < 1 km × 1 km) while introducing relevant effects on the local hydrological states (i.e. soil moisture) and fluxes (e.g. groundwater recharge). This method therefore could be considered as a valuable choice to account for the uncertainty of soil properties for this type of model applications; i.e. when well calibrated hydrological models based on streamflow measurements are used.
A different behaviour is noted for the distributed hydrological states and fluxes (evapotranspiration and groundwater recharge). These variables represent in fact local conditions (i.e. spatial resolution < 1 km) and they show the same degree of uncertainty independently from the perturbation method used. This means that these localized states and fluxes can be used to infer local properties but it is not possible to use this type of observations to calibrate the values for larger areas. For this reason, the use of, e.g., remote sensing products as total water storage anomalies and evapotranspiration is an effective approach for constraining and improving local model parameterization.

Conclusions
In the present study, uncertainty in soil properties is characterized based on three statistical perturbations methods. This uncertainty is propagated applying the distributed hydrological model mHM. The uncertainty in the simulated states and fluxes are analysed at different spatial and temporal scales. The main conclusions are summarized as follows: 1. The effect of uncertainty in soil properties depends on the hydrological model output. In particular, the uncertainty in the fluxes are relatively positive correlated; i.e. if the uncertainty in one of the simulated flux is high, also the other fluxes show, to some degrees, uncertainties. On the contrary, the uncertainty in the simulated soil moisture shows a more complex relation as its uncertainty does not always represent the overall uncertainty in the simulated fluxes. This behaviour is explained by the non-linear relation between states and fluxes and the occurrence of threshold conditions in the model response. For this reason, these results support the need for more than one variable (e.g. soil moisture and streamflow) for a proper assessment of the overall performance of hydrological models (Rakovec et al., 2016;Zink et al., 2017).
2. The uncertainty in states and fluxes depends on the specific locations and on the boundary conditions. In particular, the uncertainty in the model results shows strong temporal and spatial variability over the catchment with complex interactions to local environmental conditions (i.e. atmosphere, vegetation and soil). These results highlight the role of specific model settings (i.e. parameters and boundary conditions) for a proper characterization of the model response and the difficulty to generalize the result for other applications (i.e. different study areas and weather conditions). Similar conclusions were obtained based on sensitivity analysis conducted using hydrological models in different catchments (e.g. Shin et al., 2013;van Griensven et al., 2006) and they support the use of spatial and temporal diagnostic tools for a better understanding of the input-output space (Ghasemizade et al., 2017;Guse et al., 2016;Pianosi and Wagener, 2016;Wagener et al., 2003).
3. The uncertainty in states and fluxes depends on the spatio-temporal resolution used for the analysis. In particular, the uncertainty in all the model outputs decreases with decreasing spatial and temporal resolution.
Assuming an arbitrary threshold (e.g. CV) acceptable for a specific model application as proposed by Refsgaard et al. (2016), this scaling analysis identifies the Representative Elementary Scale (RES). On the one hand, this scale represents the resolution at which the model produces acceptable limits of predictive capabil-ity. On the other hand, it identifies the resolution above which it might become important to have a better understanding of the soil spatial variability. For this reason, this analysis proves to be a simple and practical approach for the assessment of spatially distributed models. However, in the present study the difficulties to identify a universal RES were identifies since it depends on locations, time and model output. For these reasons, the present study proposes three possible extensions of the RES approach: the use of the maximum CV, the temporal aggregation and the assessment of multi-variables. The first extension should better capture the model performance due to the strong spatial and temporal variability that could be present in the uncertainty within the catchment. The second extension could be used to emphasize the trade-off between temporal and spatial resolution of the model application. Finally, the third extension should provide a better assessment of the overall performance of the model.
4. The assumptions and the methods used for the characterizations of the uncertainty in soil properties plays a crucial role. In particular, the above conclusions are supported by the results obtained with all the three soil perturbation methods used in this study. However, the absolute value of the uncertainty detected in states and fluxes at different spatial and temporal scales strongly depends on the perturbation methods. For this reason, the results underline the importance to properly characterize the specific sources of uncertainty to transform a pure numerical exercise to specific results that are able to better support the model applications. The three methods developed and used in the present study represent three relatively simple approaches that can be considered to account for different types of uncertainty in a soil map. In particular, this study proposes a new perturbation method (here called conditional points method) able to introduce small-scale soil variability while preserving the original spatial patterns. In this context, however, the availability of soil map with additional information regarding not only the actual mean value within the soil units but also information representing the unresolved variability (variance and correlation length of the subdominant soil units) would provide strong support to hydrological modelling applications.
5. Finally, the analysis conducted in the present study identifies important information to be used for possible model improvement, either by collecting additional data regarding the soil properties or for inverse modelling and data assimilation frameworks. In particular, integrated fluxes such as river discharge of large catchments are shown not to be impacted by small-scale soil variabilities (i.e. standard deviation) but only by long spatial structures (i.e. long correlation lengths). For this reason, additional details in the soil map do not improve the model performance on streamflow but rather other sources of uncertainties should be considered for that (e.g. vegetation properties). For the same reason, this integrated observation cannot be used to infer local parameters (i.e. parameter of finer resolutions) but only mean characteristics of the input parameters (e.g. average soil properties over the soil units). On the contrary, local states and fluxes proved to be very sensitive to local variation in the soil properties (i.e. standard deviation). For this reason, a soil map with finer resolution data is found to be an important factor for decreasing the uncertainty in these local model outputs. For the same reason, these simulated outputs can be used to infer local soil parameters in calibration or data assimilation. Despite the transition between these two extreme conditions for which the uncertainty in soil properties is (or is not) important is quite smooth and it depends on the output considered and on the boundary conditions, this analysis provides a strong support to prioritize the model improvements in specific model applications. For this reason, similar studies can be considered for comparing statistical methods to characterize other sources of uncertainty relevant in catchment hydrology (e.g. precipitation, vegetation parameters).
Data availability. The providers of the model input (digital elevation model, atmospheric data, soil properties, etc.) are properly referenced in the paper and in the Acknowledgements. Additional information regarding the data availability can be obtained by contacting the authors.
The Supplement related to this article is available online at doi:10.5194/hess-21-2301-2017-supplement.