Effects of Input Data Content on the Uncertainty of Simulating Water Resources

: The widely used, partly-deterministic Soil and Water Assessment Tool (SWAT) requires a large amount of spatial input data, such as a digital elevation model (DEM), land use, and soil maps. Modelers make an effort to apply the most speciﬁc data possible for the study area to reﬂect the heterogeneous characteristics of landscapes. Regional data, especially with ﬁne resolution, is often preferred. However, such data is not always available and can be computationally demanding. Despite being coarser, global data are usually free and available to the public. Previous studies revealed the importance for single investigations of different input maps. However, it remains unknown whether higher-resolution data can lead to reliable results. This study investigates how global and regional input datasets affect parameter uncertainty when estimating river discharges. We analyze eight different setups for the SWAT model for a catchment in Luxembourg, combining different land-use, elevation, and soil input data. The Metropolis–Hasting Markov Chain Monte Carlo (MCMC) algorithm is used to infer posterior model parameter uncertainty. We conclude that our higher resolved DEM improves the general model performance in reproducing low ﬂows by 10%. The less detailed soil-map improved the ﬁt of low ﬂows by 25%. In addition, more detailed land-use maps reduce the bias of the model discharge simulations by 50%. Also, despite presenting similar parameter uncertainty (P-factor ranging from 0.34 to 0.41 and R-factor from 0.41 to 0.45) for all setups, the results show a disparate parameter posterior distribution. This indicates that no assessment of all sources of uncertainty simultaneously is compensated by the ﬁtted parameter values. We conclude that our result can give some guidance for future SWAT applications in the selection of the degree of detail for input data.


Introduction
Water availability and quality are themes of considerable concern for humanity. Predictions of the consequences of anthropogenic and natural changes in the environment are necessary to understand and support decisions regarding water-resources management, water pollution, and flood control [1][2][3]. Such predictions are made with hydrological models ranging from lumped to fully-distributed setups [4][5][6][7]. However, using a model to mimic the real world has proven to be challenging, due to In this study, we investigate the effect of different spatial input data (DEM, LULC and soil maps) on model uncertainty. We use two sources of data for the study: (1) Global model input datasets, which consider general information and may not represent the heterogeneity of the study area.
(2) Regional model input datasets, with fine resolution of the catchment characteristics, commonly not available for free and computationally demanding.
The aim of this study is to understand how the use of different sources (global and regional) of input datasets may affect the goodness-of-fit of the model and parameter posterior distribution and uncertainty. We hypothesize that finer and regionally-adapted datasets provide higher information content and that SWAT setups based on this type of data have a reduced range of parameter uncertainty and outperform those using global datasets.

Study Area
The catchment covers an area of 104 km 2 , with one-fourth located in the north-west of Luxembourg and three-fourths in the south-east of Belgium ( Figure 1). The region is part of the mountainous Ardennes region, with altitudes ranging from 324 to 564 masl. The climate is mild temperate with a total mean annual precipitation of 907 mm from 2006 to 2013. The sparsely-populated area is covered by forests (broad-leaved, coniferous, and mixed), a complex agricultural system and arable lands, pastures, water bodies, and urban areas. The LULC distribution depends on the LULC products being used for this assessment (see also Section 2.3 on LULC data inputs). The prevailing soils developed on schist and sandstone are Cambisol, and Podzols in a few places, showing an overall homogenous distribution of soil types.
Water 2018, 10, x FOR PEER REVIEW 3 of 18 In this study, we investigate the effect of different spatial input data (DEM, LULC and soil maps) on model uncertainty. We use two sources of data for the study: (1) Global model input datasets, which consider general information and may not represent the heterogeneity of the study area.
(2) Regional model input datasets, with fine resolution of the catchment characteristics, commonly not available for free and computationally demanding.
The aim of this study is to understand how the use of different sources (global and regional) of input datasets may affect the goodness-of-fit of the model and parameter posterior distribution and uncertainty. We hypothesize that finer and regionally-adapted datasets provide higher information content and that SWAT setups based on this type of data have a reduced range of parameter uncertainty and outperform those using global datasets.

Study Area
The catchment covers an area of 104 km 2 , with one-fourth located in the north-west of Luxembourg and three-fourths in the south-east of Belgium ( Figure 1). The region is part of the mountainous Ardennes region, with altitudes ranging from 324 to 564 masl. The climate is mild temperate with a total mean annual precipitation of 907 mm from 2006 to 2013. The sparselypopulated area is covered by forests (broad-leaved, coniferous, and mixed), a complex agricultural system and arable lands, pastures, water bodies, and urban areas. The LULC distribution depends on the LULC products being used for this assessment (see also Section 2.3 on LULC data inputs). The prevailing soils developed on schist and sandstone are Cambisol, and Podzols in a few places, showing an overall homogenous distribution of soil types. The study area location with its stream network. Meteorological conditions are measured in a daily resolution next to Schimpach, Luxembourg (50°00′36′′ E, 5°50′56′′ N), while the discharge of the Wiltz River is measured at the outlet close to the city of Winseler (49°58′04′′ E, 5°53′23′′ N).

Soil and Water Assessment Tool (SWAT) Model
We used SWAT to simulate hydrological fluxes in the study catchment. SWAT is a semidistributed, partly physically-based, partly empirical model that requires a diversity of specific information about weather, topography, soil properties, land-management practices, and vegetation [26]. The topographical information is used to delineate the watershed and estimate the surface area and the slope, as well as the stream network and its characteristics, such as length and width. To increase reliability of the latter, SWAT allows the use of a real river network data as burn-in. Additionally, the SWAT user has the option of defining the upstream drainage area, which is

Soil and Water Assessment Tool (SWAT) Model
We used SWAT to simulate hydrological fluxes in the study catchment. SWAT is a semi-distributed, partly physically-based, partly empirical model that requires a diversity of specific information about weather, topography, soil properties, land-management practices, and vegetation [26]. The topographical information is used to delineate the watershed and estimate the surface area and the slope, as well as the stream network and its characteristics, such as length and width. To increase reliability of the latter, SWAT allows the use of a real river network data as burn-in. Additionally, the SWAT user has the option of defining the upstream drainage area, which is required to define the beginning of a stream, providing the opportunity to define how detailed the drainage network is. We set this value to 4.74 km 2 , as recommend by the ArcSWAT interface. Soil properties and LULC information are used to estimate the land phase of the hydrologic cycle in the model. The land phase controls the amount of water, sediment, and nutrient loadings that flow to the main channel within each sub-basin. A SWAT watershed is partitioned into several sub-basins and then divided into hydrologic response units (HRUs). HRUs are land areas comprised of a unique combination of LULC, land management, soil, and slope. Unlike the sub-basins that are spatially related to one another, there is no information sharing between the HRUs in SWAT, i.e., the different HRUs are modeled separately and are then summed up to determine the total fluxes for the sub-basin in which they are located. We used the ArcSWAT 2012 version, with daily precipitation and minimum and maximum temperatures as forcing data (provided by the Agricultural Technical Services Administration of Luxembourg (ASTA)), and selected the Hargreaves method for estimating daily evapotranspiration rates [45].

Data Input and Model Setup
The model setup complexity is directly determined by the number of HRUs and sub-basins, which is indirectly defined by the heterogeneity of the input data maps. We analyzed eight different SWAT model setups (Table 1), using input data with various spatial resolution and information content combinations. For each data category, we selected a large scale, global map and a more regional data map ( Figure 2). The data used in the different setups are described in the following section.
Water 2018, 10, x FOR PEER REVIEW 4 of 18 required to define the beginning of a stream, providing the opportunity to define how detailed the drainage network is. We set this value to 4.74 km 2 , as recommend by the ArcSWAT interface. Soil properties and LULC information are used to estimate the land phase of the hydrologic cycle in the model. The land phase controls the amount of water, sediment, and nutrient loadings that flow to the main channel within each sub-basin. A SWAT watershed is partitioned into several sub-basins and then divided into hydrologic response units (HRUs). HRUs are land areas comprised of a unique combination of LULC, land management, soil, and slope. Unlike the sub-basins that are spatially related to one another, there is no information sharing between the HRUs in SWAT, i.e., the different HRUs are modeled separately and are then summed up to determine the total fluxes for the sub-basin in which they are located. We used the ArcSWAT 2012 version, with daily precipitation and minimum and maximum temperatures as forcing data (provided by the Agricultural Technical Services Administration of Luxembourg (ASTA)), and selected the Hargreaves method for estimating daily evapotranspiration rates [45].

Data Input and Model Setup
The model setup complexity is directly determined by the number of HRUs and sub-basins, which is indirectly defined by the heterogeneity of the input data maps. We analyzed eight different SWAT model setups (Table 1), using input data with various spatial resolution and information content combinations. For each data category, we selected a large scale, global map and a more regional data map ( Figure 2). The data used in the different setups are described in the following section.
The pan-European elevation data EU-DEM (D1) is a 3D raster dataset with 30 m resolution and a vertical accuracy of 2.9 m [46,47]. This hybrid product is a weighted averaging of the Shuttle Radar Topography Mission (SRTM) and the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER GDEM). 2.
The regional DEM (D2) was provided by the Luxembourg Institute of Science and Technology (LIST). The D2 data has a spatial resolution of 5 m.
To match the spatial resolution of D1, the D2 DEM was resampled to 30 m using bilinear interpolation. This reinforces our goal of focusing on different information contents rather than the effect of spatial resolution, which has been investigated elsewhere [48]. The use of different DEM sources affects the watershed delineation, resulting in different terrain characteristics and physical structures of the stream network (Table 1), which is expected to influence model performance [32,40,49]. Before discretizing the watershed into eight sub-basins, we used the same upstream drainage threshold for the area required to form the beginning of a stream (4.74 km 2 ), because Wu et al. [50] reported that different thresholds may affect the streamflow. Additionally, the watershed was classified into three land-slope classes: <7%, 7-15% and >15% that are considered moderate, medium and steep slopes, respectively, according to German guidance for soil surveying and mapping [51].
The pan-European CORINE Land Cover 2006 (L1) is a compilation of national LULC datasets. Its production was based on agreed methodology and carried out by the European Environment Agency (EEA) under the framework of the Copernicus program [52]. L1 has a 100 m spatial resolution and represents the LULC status close to the year 2006 with an accuracy above 85%. 2.
The regional LULC map (L2) titled Occupation Biophysique du Sol (2007) was also provided by the LIST.
We used the same LULC categories for both maps (global and regional), which are forest (including broad-leaved, coniferous, and mixed), pasture, agriculture, water, and urban area. Comparing the categories' distribution for both maps, the global map has a substantially larger area of agricultural land than the regional map (43.4% vs. 34.5%) and in turn, smaller shares of pasture (24.9% vs. 32.0%), forest (24.0% vs. 25.2%), surface water bodies (0% vs. 0.01%), and urban areas (7.7% vs. 8.3%). We used the same crop rotation calendar for all setups, taking yearly turns of winter wheat and corn that include fertilization, ploughing, planting, and harvest. The fertilization process comprises yearly amounts of phosphorus for winter wheat and corn (60 kg ha −1 and 48 kg ha −1 , respectively) as well as nitrogen (174 kg ha −1 and 140 kg ha −1 , respectively).

1.
The Harmonized World Soil Database (HWSD, version 1.21) (S1) is a 1 km spatial resolution global soil map produced jointly by the International Institute for Applied Systems Analysis IIASA and the Food and Agriculture Organization of the United Nations (FAO) [53]. The database includes commonly used soil parameters and texture classes. S1 presents one soil class for the Wiltz watershed having two layers, one 0.3 m deep and the other 1 m deep.

2.
SoilGrids (S2) is a 250 m spatial resolution global soil database developed jointly by International Soil Reference and Information Centre (ISRIC-World Soil Information) and other collaborators [54]. S2 has six soil layers with depths of 0.05 m, 0.15 m, 0.30 m, 0.60 m, 1.00 m, and 2.00 m. It provides standard soil properties (soil texture, bulk density, soil organic carbon content, etc.) per each grid cell and layer. To adjust the database to the necessity of SWAT input format, we categorized the cell values according to the soil texture, separating them into three classes: (i) loamy sands and silty-loamy sands with a high percentage of sand and a low percentage of clay (≤17%); (ii) silty loams with a high percentage of silt ≥50%; (iii) and sandy loams and slightly clayey loams with clay >7% and silt <50%.
After defining soil texture values for both maps (global and regional), we calculated SWAT soil physical parameters using the pedotransfer function (PTF) developed by Saxton and Rawls [55]. The soil albedo was calculated using the Baumer equation that is based on the percentage of organic material content [56].
The number of HRUs created by the SWAT model setup is based on the quantity and spatial distribution of LULC, soil, and slope classes. Interpreting the number of HRUs as the setup complexity, the simplest setups use S1 (D2S1L1 and D1S1L1) and have 89 and 90 HRUs, respectively, while the most complex setups use S2 and have 200 HRUs and 199 HRUs (D2S2L2 and D1S2L2, respectively) ( Table 1). When defining the HRUs, SWAT allows the use of a threshold for soil, LULC, and slope to reduce the complexity of the model setup. However, we are not using any threshold, because we want to avoid loss of both information and diversity representation.

Calibration, Parameters, and Uncertainty Analysis
We considered daily discharge measured at the main outlet of the catchment for model validation, For the model calibration and parameter uncertainty analysis, we used the widely accepted Metropolis-Hasting Markov Chain Monte Carlo (MCMC) algorithm [57] as implemented in the SPOTPY Python package [58]. The SPOTPY package is an open-source tool that simplifies the link between models and Bayesian, optimization, and sensitivity analysis methods. Once the user connects the model to SPOTPY, testing the effects of different analysis strategies and analyzing the model performance is straightforward. The package allows the choice among 12 algorithms for calibration, 9 pre-built parameter distribution functions, 16 objective functions, and 21 likelihood functions.
SWAT was calibrated for daily discharge from the main watershed outlet point. The selection of the six model parameters to be calibrated is based on expert-knowledge. We excluded several parameters that are widely used in SWAT calibration [59][60][61], but are closely related to one of the inputs we were concerned with analyzing (DEM, soil, and LULC maps). Examples include the curve number, saturated hydraulic conductivity, and available water capacity of the soil layer. We focused our analysis on six parameters related to snow, groundwater, and main channel processes, which are described in Table 2, including their respective prior uncertainty range.
We assumed a non-informative uniform prior distribution for the parameters, which were sampled 200,000 times in five parallel chains, maximizing a logarithmic likelihood function: where y is observed data, y is simulated data, θ is parameter set, σ is the measurement error of observed data, t is the time step, and n is total number of observations [62]. We used the multivariate approach proposed by Brooks and Gelman [63] to monitor the chain's degree of convergence (R). Convergence is achieved when the parameter variance within a chain and the variance between chains become indistinguishable, generating aR value close to one. LargeR values indicate a notable difference between chains. The posterior distribution was defined taking in consideration a warm-up period of one-fourth of the simulation length and the convergence criteriaR < 1.2 in order to generate the posterior parameter distribution.
We assessed model performance calculating the Nash-Sutcliffe Efficiency (NSE), the logarithmic NSE (log NSE), and the percent bias (PBIAS) for the posterior distribution. NSE is defined as the magnitude of residual variance normalized by the observed data variance [64]. Negative values indicate that the mean observed data is a better predictor than the model simulation, whereas values between 0.5 and 1 are considered a good model fit [65]. The squared residuals of the NSE calculation overemphasize high values [10,66]; therefore, we also considered the logarithmic NSE [67], which is more sensitive to low flows than NSE [66,68]. This metric can be interpreted similarly to NSE, where a good performance is achieved when the log NSE is greater than 0.5. PBIAS is determined by the difference between observed and simulated data, normalized by the observed discharge [69]. Model accuracy is characterized by low values, i.e., the metric has zero set as the optimum value, and values between −3 and 3 are considered a good performance. Moreover, a positive PBIAS value represents a model toward underestimation and a negative value a model toward overestimation [69].
Part of model uncertainty was estimated using the measured P-and R-factors [70]. The P-factor is the percentage of data bracket by the 95% prediction uncertainty (95PPU) which is calculated at the 2.5 and 97.5 percentiles of the simulated data. The R-factor is the ratio of the average distance between the upper and lower 95PPU and the standard deviation of the measured data. Ideally, the P-factor tends to 1 and the R-factor is close to 0.

Posterior Model Performance
The model was calibrated using six parameters described in Table 2. We focused on the parameters that are not directly affected by soil, LULC, or topographic information to facilitate the evaluation of the effects of the various spatial products. Our Bayesian calibration resulted in the formation of two groups driven by the different information content of the soil maps, S1 and S2 (Figure 3).
The NSE values range from 0.76 to 0.82 (Figure 3b), indicating a good performance for all model setups. The small variances among NSEs indicate that high flows are not very sensitive to the differences between the various investigated spatial input data. However, a more-detailed soil map (S2) results in higher NSEs and generates a better representation of high flows/events. The main difference between soil maps S1 and S2 is the soil depth, which is important information for a hydrological model to estimate storage volume and retention time of the water [71]. Gatzke et al. [72] and Ficklin et al. [73] agreed that variations in soil depths may be one of the main causes of differences in the hydrologic output for the SWAT model.
Log NSE values are substantially lower than NSE, varying between 0.28 and 0.59 (Figure 3a), and the log NSE value was greater than 0.5 for only one setup (D2S1L1). This may represent a structural model problem in representing low flows and, as suggested by Geza and McGray [74], it could indicate a bad prediction of base flow since that is dominant during low flow conditions. Furthermore, Gassman et al. [75] suggested that weak results can be attributed in part to inadequate representation of rainfall inputs. As we only had one rainfall gauging station available for our simulation, the model simulates a rainfall-runoff reaction with almost every rainfall input, contrasting the observations. This situation assumes uniformity of precipitation throughout the watershed, neglecting its spatial variability. When available, the use of multiple stations or radar could improve the model prediction [76][77][78]. Another possible explanation for the lack of fitting for low flows is the choice of the likelihood function that optimizes high flows more efficiently [68]. Taking a closer look at the log NSE distribution, setups using the regional DEM D2 present a slightly better performance than those using global DEM D1 (Figure 3a), suggesting a sensitivity to small topographic changes. However, these differences are not clear for high flows (Figure 3b), where D1S1L1 and D2S1L1 present reverse behavior for D1 versus D2 for the NSE. Overall, we conclude Taking a closer look at the log NSE distribution, setups using the regional DEM D2 present a slightly better performance than those using global DEM D1 (Figure 3a), suggesting a sensitivity to small topographic changes. However, these differences are not clear for high flows (Figure 3b), where D1S1L1 and D2S1L1 present reverse behavior for D1 versus D2 for the NSE. Overall, we conclude that the regional D2 product is superior in reproducing streamflow. Many studies agree on the streamflow sensitivity to changes of DEM. Focusing on resolution, Cotter et al. [28] and Chaubey et al. [29], demonstrated that different DEMs result in different watershed areas, stream networks, and slopes and that these characteristics interfere with streamflow prediction. Dixon and Earls [32] showed that a 90 m DEM and a resampled 90 m DEM from a 30 m DEM contributed to different accuracies when predicting streamflow, highlighting the differences in information content when dealing with DEMs with the same resolution. Li and Wong [79] concluded that lowering DEM resolution does not necessarily lead to a poor model performance with respect to the quality of the data source. However, they expressed some concern about assuming that derived lower-resolution data are always superior to another lower resolution from a different source.
Looking at the distribution of PBIAS values (Figure 3c), one can see that the more-detailed LULC map (L2) results in PBIAS values closer to zero, i.e., a slightly better representation of the overall water balance. Earlier findings confirmed that streamflow is not extremely sensitive to LULC source or resolution. Huang et al. [80] found slight differences in the simulated monthly and daily streamflow when using three different LULC maps, in time, and in number of categories. After comparing seven resolutions of the same LULC map, Cotter et al. [28] showed that streamflow is not significantly affected by LULC resolution. By contrast, Tang et al. [81] investigated afforestation measures, which reduced modeled streamflow with SWAT. However, since our LULC maps mainly differed in pasture and agricultural categories we argue that these two LULC options do not effect evapotranspiration decisively, and thus the water balance of our catchment.
When combined with global soil map S1, setups using L2 presented a better model performance for high flows and a worse model performance for low flows than those using L1. However, the opposite is true when the setups are combined with the regional soil map S2, where L1 performs better than L2 for high flows and worse than L2 for low flows. PBIAS values also indicate that SWAT underestimates streamflow when using S1 and overestimates it when using S2. Despite indicating similar model performance for different soil maps, previous studies also indicated changes in the flow prediction. Cotter et al. [28] and Geza and McCray [74] suggested that decreases in soil map resolution may result in slower predicted flow. Wahren et al. [43] showed that soil maps with different data content may provide a different representation of the peak flow after a dry period.

Uncertainty Analysis
A reliable prediction model must be accurate, i.e., measure the bias representing how close the simulated values are to the observed data, and precise, i.e., representative of the range of the simulated values.
The simulation model accuracy is estimated in Figure 3 and represents its goodness-of-fit that was discussed in the previous section. Figure 3 also shows the simulation precision pictured as the width of each boxplot, which represents the model parameter uncertainty. All setups present similar uncertainties on the simulated discharge and a higher uncertainty for low flow predictions compared to high flow predictions. Figure 4 represents the accuracy and precision of the simulation, showing the 95PPU in the hydrography, how much this interval bracket observed discharge data (P-factor), and how spread the simulated data are (R-factor). Notice that the P-factor indicates how much of the uncertainty we are capturing. All setups present similar P-and R-factor values, with the P-factor varying from 0.34 to 0.41 and the R-factor varying from 0.41 to 0.45. Confirming the difference shown in Figure 3, the hydrographs slightly differ depending on the use of different soil information, presenting smaller P-factors for S1 and greater for S2. These results agree with Shen et al. [48] and Kumar and Merwade [40] who showed a relative difference on prediction uncertainty, which could be due to the ranges of calibrated parameters. observable in the posterior parameter distribution, which is described in the following section. To solve this issue, Ajami et al. [12] and Yen et al. [15] suggested the use of a framework, IBUNE (Integrated Bayesian Uncertainty Estimator) and IPEAT (Integrated Parameter Estimation and Uncertainty Analysis Tool) respectively, to account for major uncertainties sources, such as, parameter, input data and model structure uncertainty. However, both frameworks only consider precipitation as input data, not accounting for the uncertainty of spatial input data such as DEM, LULC and soil maps.  We use the quantile-quantile plot (QQ plot) adapted from the probabilistic forecast to evaluate the uncertainty analysis [13]. The predictive QQ plot compares the cumulative distribution function (cdf) of the simulated discharge to the cdf of a uniform distribution. Consequently, it is possible to assess whether the hypotheses in the calibration framework are consistent with the observed data by analyzing the relation between the QQ plot and the identity line. Figure 5 shows the QQ plot for all setups, and the curve shape suggests that these predictive uncertainties are being underestimated [13], which means that there are further uncertainties that are not considered in our modelling approach. As Leta et al. [82] showed, once all sources of uncertainty are not considered at the same time the parameters will compensate for this lack of information. This compensation is also observable in the posterior parameter distribution, which is described in the following section. To solve this issue, Ajami et al. [12] and Yen et al. [15] suggested the use of a framework, IBUNE (Integrated Bayesian Uncertainty Estimator) and IPEAT (Integrated Parameter Estimation and Uncertainty Analysis Tool) respectively, to account for major uncertainties sources, such as, parameter, input data and model structure uncertainty. However, both frameworks only consider precipitation as input data, not accounting for the uncertainty of spatial input data such as DEM, LULC and soil maps.  Table 1).

Posterior Parameter Distribution
The parameter posterior distribution is shown in Figure 6. We noticed a constrained behavior for all calibrated parameters, highlighting the model sensitivity to these parameters. As with the pattern presented in the model efficiency criteria, the parameter posterior distribution separates into two groups with different parameters modal values, and these groups are driven by the soil information S1 and S2. All setups using S1 have similar parameter posterior distribution and are different from the parameter distributions for the S2 setups. Kumar and Merwade [40] and Bossa et al. [83] agree that despite presenting similar model performance, some parameters compensate for the differences in soil maps.
If the daily mean air temperature is less than the snowfall temperature model parameter (SFTMP), then the precipitation is classified as snow in the SWAT model. The corresponding amount of water is then rooted in a snow pack. The snow pack does not melt until its temperature exceeds the snow melt base temperature model parameter (SMTMP). Thereafter, the snowmelt is added to rainfall runoff and the percolation calculation [71]. SFTMP and SMTMP can vary between −5 and 5°C. The posterior distributions show that The SFTMP modal value varies between 0.1 and 4.5 (Figure 6a), and the modal SMTMP value varies between 0.6 for setups using S2 and 3.3 for setups using S1 (Figure 6b).  Table 1).

Posterior Parameter Distribution
The parameter posterior distribution is shown in Figure 6. We noticed a constrained behavior for all calibrated parameters, highlighting the model sensitivity to these parameters. As with the pattern presented in the model efficiency criteria, the parameter posterior distribution separates into two groups with different parameters modal values, and these groups are driven by the soil information S1 and S2. All setups using S1 have similar parameter posterior distribution and are different from the parameter distributions for the S2 setups. Kumar and Merwade [40] and Bossa et al. [83] agree that despite presenting similar model performance, some parameters compensate for the differences in soil maps.
If the daily mean air temperature is less than the snowfall temperature model parameter (SFTMP), then the precipitation is classified as snow in the SWAT model. The corresponding amount of water is then rooted in a snow pack. The snow pack does not melt until its temperature exceeds the snow melt base temperature model parameter (SMTMP). Thereafter, the snowmelt is added to rainfall runoff and the percolation calculation [71]. SFTMP and SMTMP can vary between −5 and 5 • C. The posterior distributions show that The SFTMP modal value varies between 0.1 and 4.5 (Figure 6a), and the modal SMTMP value varies between 0.6 for setups using S2 and 3.3 for setups using S1 (Figure 6b).
The Manning's value for the main channel (CH_N2) is an empirical coefficient related to roughness, and it can vary between 0.01 and 0.25 mm h −1 . The posterior distributions show a modal CH_N2 value varying between 0.11 and 0.12 for setups using S1 and between 0.08 and 0.10 for setups using S2 (Figure 6c). These values may represent "very weedy reaches, deep pools, or floodways with a heavy stand of timber and underbrush" [71].  Table 2 for acronyms and definition of model parameters). The different model setup combinations refer to digital elevation maps (D), LULC maps (L) and soil maps (S) (see Table 1 for further description). Results are grouped by the use of a global soil map (S1) with a solid line and a regional soil map (S2) as a dashed line.
The base flow recession constant (Alpha_BF) is a direct index of groundwater flow response to changes in recharge [71,84]. Alpha_BF can vary between 0.001 and 0.99. The smaller its value, the slower the modeled land response to recharge is. Setups using S2 present a modal Alpha_BF value varying between 0.35 and 0.39, and setups using S1 present two peaks for the modal value of Alpha_BF, one between 0.22 and 0.25 and another between 0.67 and 0.71 (Figure 6e). Groundwater delay time (GW_delay) is the lag between the time that the water exits the soil profile and enters the shallow aquifer [71]. It can vary between 0 and 31 days, and dry regions with deep water tables may present long time delays. Posterior distributions show that GW_delay has two peaks for the modal value for each soil information: varying between 0.6-0.9 and 10.2-12.9 for S1 setups and between 0.8-1.1 and 5.3-7.0 for S2 setups (Figure 6f). The groundwater parameters are primordial for base flow prediction, a process that is dominant during low flow conditions. Our previous results showed that the model does not perform well for low flow events, which may explain the differences on groundwater parameters' posterior distribution. These parameters are compensating for part of the  Table 2 for acronyms and definition of model parameters). The different model setup combinations refer to digital elevation maps (D), LULC maps (L) and soil maps (S) (see Table 1 for further description). Results are grouped by the use of a global soil map (S1) with a solid line and a regional soil map (S2) as a dashed line.
Effective hydraulic conductivity in the main channel alluvium (CH_K2) describes the infiltration through the channel bottom and the higher its value, the higher the water loss rate from the river channel. CH_K2 can vary between 0.001 and 150 mm h −1 . The posterior distribution shows that the CH_K2 model values vary between 134 and 136.6 for setups using S1 and are lower, between 116.6 and 131.3, for setups using S2 (Figure 6d). CH_K2 can be separated into groups according to the alluvium materials where high CH_K2 values, like those we obtained, correspond to an alluvium material formed by clean gravel and sand [71].
The base flow recession constant (Alpha_BF) is a direct index of groundwater flow response to changes in recharge [71,84]. Alpha_BF can vary between 0.001 and 0.99. The smaller its value, the slower the modeled land response to recharge is. Setups using S2 present a modal Alpha_BF value varying between 0.35 and 0.39, and setups using S1 present two peaks for the modal value of Alpha_BF, one between 0.22 and 0.25 and another between 0.67 and 0.71 (Figure 6e). Groundwater delay time (GW_delay) is the lag between the time that the water exits the soil profile and enters the shallow aquifer [71]. It can vary between 0 and 31 days, and dry regions with deep water tables may present long time delays. Posterior distributions show that GW_delay has two peaks for the modal value for each soil information: varying between 0.6-0.9 and 10.2-12.9 for S1 setups and between 0.8-1.1 and 5.3-7.0 for S2 setups (Figure 6f). The groundwater parameters are primordial for base flow prediction, a process that is dominant during low flow conditions. Our previous results showed that the model does not perform well for low flow events, which may explain the differences on groundwater parameters' posterior distribution. These parameters are compensating for part of the model structure uncertainty.

Conclusions
This study builds on previous studies by others, which provided detailed insights into the use of different spatial input datasets. So far, a detailed intercomparison of DEM, LULC and soil maps was missing in order to address the influence of different spatial input datasets (DEM, LULC and soil maps) on SWAT performance estimating discharge and its uncertainty. Spatial data originating from different sources may provide different data content despite having the same resolution. We analyzed eight different setups using global (D1, L1, and S1) and regional (D2, L2, and S2) DEM, LULC, and soil maps. Similar to previous studies by Chaubey et al. [29] who investigated DEM in different resolution and Yen et al. [35] who studied the effect of different land-use maps, our results indicate that the choice of regional or global information may depend on the focus of the analysis, because SWAT performance varies for high and low flows. SWAT predicted high flows efficiently for all setups, suggesting that the model is minimally affected by spatial input data differences for high flows. However, there is a notable difference among the setups when predicting low flows. When comparing all setups, those using the regional D2 and the global S1 are always more efficient in representing low flows, suggesting that SWAT is either sensitive to small topographic changes but cannot necessarily make use of additional soil information, or that there are other model errors that we have not considered. The model setup D2S1L1 represents the only setup with a log NSE greater than 0.5, indicating that additional data or missing processes are needed to improve the model's capability to simulate low flows in this catchment.
Each parameter presents a similar posterior distribution according to the soil map used, indicating smaller parameter values for the setups using the regional S2 map. Additionally, there is a constrained behavior of all parameter posterior distributions for all setups, highlighting the model sensitivity to these parameters. All setups present similar uncertainty on the output and higher uncertainty for low-flow than for high-flow predictions. Furthermore, the QQ plot results show that all model setups underestimate the model uncertainty, suggesting that additional sources of uncertainties should be considered simultaneously. To proceed further in addressing all sources of modelling uncertainty, one promising way forward would be to combine our methodology with the work presented by [12,15] with their IPEAT and IBUNE frameworks in order to account for meteorological input data and model structure uncertainty.
Author Contributions: C.C. performed the calculations, analyzed the results and wrote the paper; S.J. provided expertise on SWAT, supervised the model setup and provided expertise on the area of study; T.H. developed the parameter estimation tool, supervised the work and contributed to the paper writing; M.B. contributed to the paper writing; L.B. provided the general idea of this paper, supervised the work and contributed to the paper writing.