Interactive comment on “ A comprehensive sensitivity and uncertainty analysis for discharge and nitrate-nitrogen loads involving multiple discrete model inputs under future changing conditions ” by Christoph Schürz et al

In this manuscript by Schürz et al., a detailed analysis of the impact of scenario simulations on hydrological variables is presented. In a sensitivity analysis, the sensitivities of five groups are separated. These groups are three types of scenarios (land use, point source and climate) and two model-specific groups (model set-up, model parameterisation). In the analysis, the impact of the input variables and of the uncertainty of the selection of scenarios or model characteristics is presented.


Introduction
Environmental systems are under constant change. Predicting the development of natural resources in a changing system 10 involves large uncertainties (Milly et al., 2008). Climate change, in concurrence with other dynamic processes such as population growth, land use change or economic development pose challenges to the management of water supply and water quality (Duran-Encalada et al., 2017;Yates et al., 2015). Human disturbances can exacerbate the impacts of climate and amplify consequences to water quality (Jiménez et al., 2014) on one hand. On the other hand, stakeholders in environmental systems have to respond to future changes, for instance adapting farm management practices due to changes in temperatures 15 and precipitation patterns (Schönhart et al., 2018). Ideally, an impact assessment considers all future changes that can affect the development of the environment of interest as well as those future changes that can introduce uncertainties in the simulation of the environmental variables of interest.
Changes in environmental systems are typically represented by discrete scenarios in impact studies. Preferably, the set of scenarios representing a dynamic change covers the full range of trajectories along which the development is plausible 20 (Clark et al., 2016). Scenario development involves different data sources and models, which can introduce and propagate uncertainties. For example, climate scenarios have several sources of uncertainty and may include several socioeconomic scenarios (e.g. the current "Representative Concentration Pathways" (RCP; Moss et al., 2010;van Vuuren et al., 2011)) that drive an array of global climate models (GCM) (Knutti and Sedláček, 2013). However, the GCMs also have inherent uncertainty and they provide the boundary conditions for regional climate models (RCM) (e.g. Jacob et al., 2014). Further, the downscaling 25 (Wilby et al., 1998;Wood et al., 2004) of the RCM simulations and the bias correction Seibert, 2013, 2012) are associated with their own uncertainty and are a standard procedures in climate scenario development. Eventually, it is essential to characterize the uncertainties inherent in all processes that affect the simulation of an environmental variable.
To simulate the development of hydrological variables under changing conditions, the developed scenarios are implemented as boundary conditions in hydrological models that are calibrated for historic observations. Yet, often different model setups 30 and different sets of parameters in a model can perform equally well to reproduce historical observations of the variables of interest. Equifinality is a well-known issue in hydrologic modeling that has been extensively addressed in the literature (e.g. Schulz et al., 1999;Beven, 2006;Beven and Freer, 2001;Beven, 1996), where multiple model structures (e.g. Clark et al., 2008) and model parametrizations (e.g. Schulz et al., 1999) represent observations equally well and thus cannot be rejected (Beven, 2006). An adequate representation of historical data does not necessarily assure that different model setups agree when extrapolating to future conditions (Chiew and Vaze, 2015;Milly et al., 2008). Thus, differences in the model setup are a source of uncertainty in the simulation of an environmental variable under future conditions. Altogether, an impact study comprises an abundance of combinations of trajectories of system changes and model setups 5 to describe an environmental system that ultimately characterize the uncertainties in a simulation. Hence, a comprehensive description of the uncertainties in model simulations is a major challenge of any impact study.
Model sensitivity analysis (SA) can be used to derive the impact of different input variables on hydrological target variables.
SA investigates the response of a modeled variable to the variation of model input variables (Saltelli et al., 2004). For a local sensitivity analysis (LSA) the model inputs are varied around a point (often an 'optimum' point) in the model input space.
Global sensitivity analysis (GSA) assesses the sensitivity of a model output for the entire feasible range of model inputs (Gupta and Razavi, 2017;Pianosi et al., 2016). Compared to LSA, GSA usually requires a larger number of computations. Thus, a substantial part of recent GSA literature focuses on the computational efficiency and the robustness of GSA methods (e.g. Pianosi and Wagener, 2015;Razavi and Gupta, 2016a;Sarrazin et al., 2016;Cuntz et al., 2015;Rakovec et al., 2014), but also on increasing the insight into modeled systems from a certain number of model evaluations (e.g. Borgonovo et al., 2017;Dai 15 et al., 2017;Guse et al., 2016a;Massmann et al., 2014;Razavi and Gupta, 2016a).
The complexity and computational demand of a model determine the feasible number of model evaluations and thereby the applicability of a SA method (Razavi and Gupta, 2015). Large atmospheric model applications, for instance, only allow a LSA with a few model evaluations (Gupta and Razavi, 2017;Pianosi et al., 2016). Environmental model applications are usually less computationally expensive and allow a more extensive GSA, illustrated in many environmental modeling studies 20 (e.g. Guse et al., 2016b;Massmann and Holzmann, 2015;Razavi and Gupta, 2016b;Sarrazin et al., 2016). Most applications utilize GSA to identify influential model parameters and to rank model parameters according to their influence on model outputs. Model parameters are usually continuous model inputs. (Saltelli et al., 2008;Baroni and Tarantola, 2014).
Although it is possible to implement composite model inputs (e.g. climate scenarios that affect several climate variables at 25 the same time, or land use scenarios that can impact the entire model setup) in a GSA and therefore to employ GSA in impact studies, a consideration of discrete and composite model inputs can constrain the applicability of GSA and complicate the implementation (Baroni and Tarantola, 2014). In impact studies, the response of an environmental variable to a (future) change in a model input is usually inferred by implementing a scenario realization of the respective model input in a model setup.
From an SA perspective, this approach is equivalent to a local 'one-at-a-time' (OAT) assessment of the model input sensitivity 30 (Saltelli and Annoni, 2010;Baroni and Tarantola, 2014). A local OAT analysis however presumes linear models and noncorrelated inputs which are hardly true for any environmental model application (Rosolem et al., 2012;Baroni and Tarantola, 2014). Thus, to account for interactions of model inputs and model non-linearities the application of GSA is recommended instead (Saltelli and Annoni, 2010;Saltelli and Tarantola, 2002; Baroni and Tarantola, 2014).
Yet, a few studies implemented discrete and composite model inputs in GSA. With the Generalized Probabilistic Framework, Baroni and Tarantola (2014) rendered a solid basis for the implementation of correlated, non-continuous model inputs in GSA and applied the variance-based SA method of Sobol (1993) to assess the response of soil moisture, evapotranspiration, and soil water fluxes to uncertainties in meteorological input data, crop parameters, soil properties, model structure, and observation data. In a synthetic example, Dai and Ye (2015) performed model and scenario averaging to assess the impact of different 5 model structures and scenarios of precipitation on groundwater flow and reactive transport in the soil. In a more recent study, Dai et al. (2017) employed the method of Sobol to identify the relevant system processes for groundwater flow and reactive transport represented in different model structures. Savage et al. (2016) applied GSA to identify the dominant controls in the calculation of flood inundation, to assess whether a high spatial resolution of the flood inundation model or the model parametrization is dominating the simulation. The mentioned studies illustrate the use of GSA with discrete and composite 10 model inputs. Anderson et al. (2014) and Butler et al. (2014) highlight the importance of assessing the uncertainty of future climate change impacts and the identification of relevant drivers and their interactions for climate policy making.
In this paper we demonstrate the utility of GSA and uncertainty analysis in a comprehensive setting of an environmental model impact study and address the following points: -We apply GSA in two environmental modeling impact studies to identify the dominant sources of uncertainties for the simulation of discharge and nitrate-nitrogen (NO − 3 -N) loads. We analyze the impacts of different spatial aggregations of the model setup and different model parametrizations and assess the effects of changes in the land use, point source emissions, and the future climate.
-We analyze the resulting uncertainties in the simulation of the long-term monthly mean discharge and monthly sums of NO − 3 -N loads, as well as flow duration curves (FDCs) of daily discharge and daily NO − 3 -N loads visually. We present 20 ways to visualize the discrete model inputs that provide further insights into the relationships of uncertainties in the simulations and different properties of the discrete realizations of the model inputs.
-Based on the GSA and the visual analysis of the simulated uncertainties we are able to draw conclusion on the simulation of discharge and NO − 3 -N loads as impacted by the model setup, model parametrization and the future scenarios of land use, point source emissions and climate. These conclusions are of course limited to assumptions made in the model setup 25 and in the development of the scenarios.
The paper is structured in the following way: Section 2 contains an overview of the two investigated catchments, the Soil and Water Assessment Tool (SWAT, Arnold et al., 1998) that we implemented in this study, and the preparation of the model input data that we used in the model setup. In Section 2.4 we describe the setup of the SWAT model with different spatial aggregations and illustrate the pre-processing of the SWAT model setups that was necessary to identify the sensitive SWAT model parameters 30 and to define non-unique parameter sets for all model setups. The scenarios of land use, point source emissions and the climate together with the input data and pre-processing to develop the individual scenarios are specified in Section 2.5. Section 2.6 combines the SWAT model setups, the defined non-unique model parametrizations and the developed scenarios of land use, point source emissions and climate in the GSA and explains the methods we applied to analyze the sources of uncertainties for the simulation of discharge and NO − 3 -N loads. The results of the combined GSA framework and the visual analysis are provided in Section 3. We discuss the findings of the GSA application and the visual analysis of the simulation uncertainties for the two case studies in Section 4 and address the specific assumptions that we made during the model setup and the development of the scenarios. The Schwechat river has its source in the Vienna woods at the northeastern boundary of the Northern Limestone Alps with a maximum altitude of 893 m a.s.l. After a natural flow section in the narrow and dominantly forested valley of the "Helenental" (70% of the total catchment area. See Table A3), the Schwechat drains into the Vienna basin with flat topography 15 and a predominance of agriculture, viniculture and settlement areas. The main agricultural crops are winter wheat and summer wheat. Larger areas in the upper part of the catchment are used as pastures (~10% of the total area). The largest settlement is the city of Baden with a population of approximately 26000 inhabitants, while smaller settlements are scattered over the catchment. All municipal wastewaters are collected in three wastewater treatment plants (WWTP, black triangles in Fig. 1 where the the WWTP Baden is the most relevant one with a capacity of 45000 population equivalents (PE). All WWTPs 20 perform carbon removal, nitrification, denitrification and enhanced phosphorus removal. Due to the close proximity to the city of Vienna population growth is a likely prospect for the settlement areas in the lower part of the catchment. The part of the catchment considered in this study has its outlet next to the city of Traiskirchen at an altitude of 185 m a.s.l. and covers an area of approximately 275 km². The long term mean annual precipitation in the Vienna Basin is around 620 mm/yr and the mean annual temperature is 9.9°C. 25 The Raab river originates at the edge of the southeastern Alps. These are characterized by low mountain ranges with a maximum altitude of 1547 m a.s.l., mostly covered by forests (~42% of the total catchment area. See Table A4). The Raab flows through the southern part of Austria and crosses the boarder to Hungary close to the city of Neumarkt an der Raab at an altitude of 232 m a.s.l. The case study encompasses the Austrian part of the Raab with a catchment area of approximately 998 km². The long-stretched river valley is dominated by agricultural activities (~25 % of the total area), with urban areas in 30 between. The slopes along the Raab are covered with heterogeneous patterns of forests, pasture areas and agricultural land use. The main agricultural crops are corn and oil seed pumpkins, but also wheat and vegetable production are common. While the urban areas are of similar small structure as in the Schwechat catchment, leather industries are present in the catchment that release substantial nutrient inputs into the receiving waters, which has resulted in trans-boundary conflicts (Ruzicka et al., 2009). Municipal wastewaters in the Raab catchment are collected in 12 relevant WWTPs (black triangles in Fig. 1) that all have the same standards for wastewater treatment as in the Schwechat catchment, but have almost three times the total capacity (approximately 150000 PE). Six relevant industrial emitters are located along the main reach of the Raab river (white triangles in Fig. 1) that all perform internal waste water treatment following the respective industry-specific regulations for 5 wastewater treatment (e.g., BGBl. II Nr. 10/1999BGBl. II Nr. 12/1999. The average annual precipitation in the Raab catchment is approximately 800 mm/yr and the long term annual mean temperature is 9.0°C.

The Soil and Water Assessment Tool (SWAT)
The SWAT model (Arnold et al., 1998) is a continuous, process based, semi-distributed eco-hydrological model. In this study we implemented SWAT2012 (Rev.622) to simulate daily time series of discharge and NO − 3 -N loads at the catchment outlets. 10 The models' spatial reference to a catchment is given by a subdivision of the basin into subbasins. Areas containing the same land use, soil type and lying in the same slope range are lumped together in each subbasin to form hydrologic response units (HRUs). All processes on the land phase of each subbasin are calculated at the HRU scale and are further propagated into the water phase of each subbasin. The processes calculated on the land phase include water balance components such as interception, infiltration, shallow and deep percolation, surface runoff, lateral flow, groundwater flow, plant uptake and 15 evapotranspiration, or the pathways of nutrients such as the input through atmospheric deposition, or fertilizer application, the transformation into other forms of a nutrient and the transport through surface runoff, percolation, lateral flow and return flow in the groundwater (Neitsch et al., 2011). In the water phase, the nutrients budgets are calculated. Following the calculation of the water balance and the nutrient budgets, the discharge, the nutrient loads and other substances are routed through the linked subbasins to the defined catchment outlet (Neitsch et al., 2011). The required input data to set up a model with SWAT are a digital elevation model (DEM), a raster land use map including the model parametrization and the performed management operations for each land use, a raster soil map with soil physical and chemical parameters for all soil layers, and meteorological 5 input data.

Model input data and data preparation
A DEM with a 10 m resolution was available for Austria from an airborne laser scan (Geoland.at, 2015). Based on the DEM we defined three slope classes with slopes of 0-3%, 3-8%, and >8% in the HRU definition step.
CORINE land cover (EEA, 2015) served as the base land use map to which more detailed agricultural data was added.

10
CORINE does not classify agricultural land uses into crop types. Therefore, tabular data of agricultural land uses at the municipal level derived from the 2010 Austrian agronomic census (Statistik Austria, 2015b) was superimposed onto CORINE data by randomly distributing crops according to the crops' areal share at the municipal level to CORINE pixels containing agricultural and complex cultivation land use. Typical time windows for planting, fertilizer application, tillage and harvest were derived from field experiment records for the individual crops (Land NÖ, 2015) and written to the HRU management files. The 15 management dates were randomized for all HRUs within the time windows derived for a management operation. Dates with strong rainfall or a high soil moisture potential were not used for scheduling management operations. With 70.0% and 42.3% forest land uses were the most dominant land uses in the Schwechat and the Raab catchments, respectively. The SWAT model setups differentiated between deciduous forests, coniferous forests and mixed forests, derived from CORINE land cover (see Tables A3 and A4). All HRUs with one of the three forest types as land use were parameterized with an initial biomass and an 20 initial leave area index to simulate intact forests in both catchments.
The SoilGrids data base (Hengl et al., 2017) is a consistent global soil information system that provides soil physical and chemical parameters at a 250m grid resolution and seven soil depths. We utilized the available soil parameters from SoilGrids and estimated further required soil parameters with pedo-transferfunctions provided by the R package euptf (Tóth et al., 2015).

Meteorology
INCA (Haiden et al., 2011) Preciptation and temperature data in 1km resolution. and 86.3% of municipal point source emissions in the Schwechat and the Raab catchments respectively. Thus, these data cover a substantial part of the municipal emissions. Additionally, daily and weekly internal monitoring data was available for some large WWTP schemes. In most cases however, only information on NO − 3 -N emissions was provided. A general budgeting of nitrogen emissions however showed, that the substantial share of total nitrogen is emitted in form of NO − 3 -N (87% in the Schwechat catchment and 89% in the Raab catchment). For industrial emitters monthly and annual records from internal 5 and external monitoring agencies were available and only allowed an estimation of industrial emissions with coarse temporal resolution, while covering the annual budgets. Again, mainly data for NO − 3 -N emissions were available. Although, nitrogen is emitted in different forms the available data basis only allowed to consider NO − 3 -N loads contributed by point sources. Table 1 provides an overview of the model input data that was used for the SWAT model setup.
Hourly observations of discharge were available for the period from 2003 to 2015 at two gauges for the Schwechat and To compare the observations with the modeled SWAT outputs of discharge and NO − 3 -N loads, daily NO − 3 -N loads and daily mean discharge were calculated from the observation data.  Graphical GIS user interfaces such as ArcSWAT (Winchell et al., 2015) or QSWAT (Dile et al., 2016)  In total, we set up four SWAT models, two with 3 and two with 14 subbasins for the Schwechat catchment and six models for the Raab catchments with two each of 4, 29, and 54 subbasins. For the full HRU setups we kept the resulting HRUs unmodified.
For the model setups with a reduced number of HRUs we eliminated small HRUs. We determined thresholds for land use, soil, and slope classes to remove HRUs that have an area below these found thresholds. The thresholds were determined using 15 the R package 'topHRU' (Strauch et al., 2016). 'topHRU' enables to find thresholds that minimize the number of HRUs of a SWAT model setup while minimizing the aggregation error (sum of changes in the areas of land uses, soils and slope classes of the reduced set of HRUs compared to the full HRU setup). To maintain a comparability between the reduced HRU setups thresholds were selected that result in an aggregation error of maximum 5% in all reduced HRU model setups. Table 2 gives an overview of the final model setups for both case studies. 20 In a parameter screening, we applied a GSA to the simulations of discharge and NO − 3 -N loads at the catchment outlets of all SWAT model setups to identify influential model parameters. Initially, 42 model parameters were selected that are frequently calibrated in SWAT model setups to simulate discharge and NO − 3 -N loads (see e.g. Arnold et al. (2012) and Abbaspour et al. (2007) for a general overview of relevant model parameters, Mehdi et al. (2018) and Haas et al. (2016) for parameters controlling the water balance and nutrient cycles, or Haas et al. (2015) for a review on the dominant nitrogen parameters). The SWAT model setup initializes the model parameters using values obtained from the SWAT data bases (either standard values or user defined, e.g. by pedotransfer functions). The selected initial ranges to modify the model parameters and the selected types of parameter changes (e.g. replace parameter values globally or modify a spatially distributed parameter field by a fraction of a 5 parameter) reflect typical procedures often found in SWAT model calibration studies. An overview of the model parameters that were identified as influential and that were further used in the model impact study is provided in Table A1.
We employed the STAR VARS approach (Razavi and Gupta, 2016a, b) to screen and rank the model parameters. STAR VARS utilizes variograms along each model input dimension of the input space to infer each model inputs influence on a target variable over different scales (where short lag distances approximate the derivative based method of Morris (Morris, 1991) and 10 long distances the method of Sobol (Sobol, 1993)). The calculation of the variograms is based on the tailored STAR sampling design where "star center" points are randomly sampled in the input space. For each center point cross sections are sampled along the input factor dimensions with an equally spaced interval. For each sampled input combination the model is evaluated and variograms along the response surface are calculated. Razavi and Gupta (2016a) proposed integrated measures of the variograms as measures of sensitivity, where the measures IVARS 10 , IVARS 30 , and IVARS 50 represent the integrals over 10%, 15 30%, and 50% of each input dimension respectively and therefore provide the sensitivity of a target variable to a model input over different scales. A detailed description of the method is provided in Razavi and Gupta (2016a) and the STAR sampling is outlined in Razavi and Gupta (2016b). The method proved to be robust and computationally efficient for high dimensional problems (e.g., Razavi and Gupta, 2016b;Sheikholeslami et al., 2019;. 20 We drew STAR samples (Razavi and Gupta, 2016b) with 50 center points and ten parameter samples per parameter dimension that resulted in 18950 parameter combinations per model setup. The Nash Sutcliffe Efficiency criterion (NSE, Nash and Sutcliffe, 1970), the Kling Gupta Efficiency criterion (KGE), including its three components , and a refined version of the Index of Agreement (Willmott et al., 2012) were used to evaluate the simulated time series of daily mean discharge and daily sums NO − 3 -N loads. Additionally, we applied the ratio of the root mean square error and standard deviation 25 (RSR, (Moriasi et al., 2007)) to evaluate different segments of the FDCs of daily discharge and daily NO − 3 -N load simulations (Pfannerstill et al., 2014;Haas et al., 2016). All calculated criteria were included in the parameter sensitivity analysis as target variables. A model parameter was considered to be sensitive if it showed a relative sensitivity of 10% compared to the most sensitive parameter with respect to a specific objective criterion for at least one of the employed objective criteria.
The performed GSA for the model parameters of the different model setups of the Schwechat catchment and the Raab 30 catchment respectively showed very similar results independent of the number of subbasins and HRUs of the individual model setups (Fig. A1). Therefore, for the impact study the same set of model parameters was considered as influential for all model setups of the Schwechat and the Raab, respectively. In total, 19 parameters for the Schwechat and 16 parameters for the Raab were identified to be influential for the analyzed target variables (Table A1). The majority of parameters were identified as influential parameters in the Schwechat and the Raab case study. The parameters SNO50COV, CANMX, CDN, and SDNCO 35 were only relevant for the model setups in the Schwechat and the parameter OV_N was only influential for in the Raab. For the majority of these parameters it is a matter of the selected threshold that defines a parameter to be influential or not. The most dominant parameters were however identified as highly relevant in both case studies.
To represent the model parametrization as an input in the subsequent sensitivity and uncertainty analysis of the environmental impact study, non-unique parameter sets were identified for the Schwechat and the Raab catchments, respectively. The 5 preceding parameter SA revealed that changes in the model parameter values influenced the simulations similarly independent of the subbasin and HRU configurations in the Schwechat and the Raab catchment, respectively. As a consequence, but also to facilitate the separation of the effects of the model setup and the model parametrization in the analysis, we selected parameter combinations as non-unique ones that result in simulations of daily discharge and NO − 3 -N loads that fulfill certain objective criteria together with all model setups of the Schwechat and the Raab, respectively. For the respective 19 and 16 influential 10 model parameters we randomly sampled 100000 parameter combinations and simulated daily discharge and NO − 3 -N loads with all model setups of the Schwechat and the Raab catchments. We evaluated the simulations with the following criteria to accept a parameter set: KGE > 0.5 for daily discharge at the catchment outlets, KGE > 0.4 for daily NO − 3 -N loads at the gauges with longer continuous records (in both case studies the gauging point within the catchment and not at the catchment outlet), percentage bias (Gupta et al., 1999) < 50% for NO − 3 -N loads, and absolute RSR < 1 for different discharge and NO − 3 -N load 15 (according to Pfannerstill et al., 2014;Haas et al., 2016). In total, we identified 43 and 52 behavioral parameter combinations for the Schwechat and the Raab catchments, respectively. The ability of the selected parameter sets used with the different model setups to reproduce the observed data is illustrated in Fig. A2. The initial and final ranges of parameter changes are shown in Table A2. The 43 and 52 parameter combinations are additionally illustrated in parallel coordinate plots for the Schwechat and the Raab in Fig. A3 to show any clustering of individual parameters and interactions between parameters. 20 The majority of parameters are scattered randomly and do not show any clustering or interaction with other parameters. The parameters RCN and NPERCO in the Schwechat catchment show a clear inverse relationship. This implies that the parameters compensate each other in the behavioral model setups. This finding seems plausible for the Schwechat catchment where the NO − 3 -N transport into the receiving waters is strongly groundwater driven and a surplus of NO − 3 -N input is reduced by a decrease in NO − 3 -N percolation. The parameters SLSOIL, SURLAG, and SOL_AWC show a clear bimodal pattern for the Raab 25 catchment. The bimodal patterns of these parameters are strongly related and a compensation effect between these parameters is visible. Model setups with increased slope values (SLSOIL) and longer lag-times of the surface runoff (SURLAG) together with an increased soil available water content (SOL_AWC) resulted in behavioral model and were able to reproduce historic discharge and NO − 3 -N records, similar to the model setups where such clear relationship is not visible.

Scenario definition 30
The study involves future changes of the land use, point source emissions, and the climate. The uncertainties of these variables are expressed as discrete scenarios.
For the land use change scenarios, two scenario story lines (Rounsevell and Metzger, 2010) were developed for the Schwechat and the Raab catchments. A "business-as-usual" scenario extrapolates the trends that we determined for the dominant crops in the time period 1970(Statistik Austria, 2015b to the future (2071 to 2100), while a second "extensive" scenario assumes an extensification of agricultural activities and other intensive land uses in both catchments (Table A5).
In the Schwechat catchment population growth is the strongest factor for a future change in land use (Statistik Austria, 2015aAustria, , 2016. Hence, a transformation from extensive pasture land (-35%) to urban land use and an increase of dense urban areas describe the "business-as-usual" scenario. The "extensive" scenario assumes no change in population and a shift of half 5 of the wheat producing area to extensive pastures.
Since 1970, the areas for corn production increased by 220% in the Raab catchment, mostly for biogas production and at the expense of sugar beets and cereals (Statistik Austria, 2017). For the "business-as-usual" scenario, an increase in the corn area by a further 100% until the end of the century was assumed, replacing extensive pastures (-75%), sugar beets (-80%), legumes (-70%), and winter wheat (-30%). 10 Groundwater protection measures lead to strict regulations for fertilizer application in the Leibnitzerfeld region adjacent to the Raab catchment (LGBl. Nr. 39/2015. Therefore, the "extensive" scenario assumes an adoption of similar nitrogen regulations in the Raab catchment. Thus, decreasing areas with intensive fertilizer application, such as corn by 50% and transforming these areas to extensive pasture land was carried out in this scenario. Two municipal point source emission scenarios for both case studies (Table A6) and two industrial point source emission 15 scenarios for the Raab catchment (Table A7)  Future climate change was considered with 22 downscaled and bias corrected climate change scenarios (Table A8). Regional climate simulations were obtained from the EU-CORDEX project (Jacob et al., 2014), providing 11 GCM-RCM simulations 30 for the emission scenarios RCP4.5 (Smith and Wigley, 2006;Wise et al., 2009) and RCP8.5 (Riahi et al., 2007) respectively.
In this study we utilized daily precipitation sums and daily minimum and maximum temperatures for the time period 2071 to 2100. The EURO-CORDEX climate simulations are available at a spatial resolution of 12.5 km (EUR-11) (Jacob et al., 2014).
Statistical downscaling (Zorita and Von Storch, 1999) was applied to prepare all climate simulations at a resolution of 1 km.
To correct downscaling errors (e.g. Haslinger et al., 2013;Muerth et al., 2013), bias correction (Teutschbein and Seibert, 2013)  was applied to the climate simulations employing quantile mapping (Hempel et al., 2013). Downscaling and bias correction were performed for the historical period 1971 to 2000, involving the reanalysis datasets SPARTACUS (Hiebl and Frei, 2016) for minimum, mean and maximum temperature and GPARD (Hofstätter et al., 2013) for daily precipitation sums. Table 3 summarizes the land use change, point source emissions, and climate change and the model setups and model 5 parametrizations that were used for the analysis of simulated discharge and NO − 3 -N loads in the Schwechat and the Raab catchments. In total, 7000 combinations of land use, point source emissions, climate, model setups and model parametrizations were drawn for both case studies applying a quasi-random sampling Saltelli and Tarantola (2002). The number of combinations results from previous experiments that applying the SA method of Sobol (results not shown) using the sampling strategy proposed by Saltelli and Tarantola (2002). A base sample size of N b = 1000 was used to meet the suggestions shown in Saltelli 10 et al. (2008). Thus, the total sample size of 7000 is defined as N = N b (k + 2), where k is the number of model inputs (k = 5).

Analysis
Although (Sarrazin et al., 2016) report publications that required substantially larger base sample sizes (e.g. N b = 12000 in Nossent et al. (2011), or N b = 8192 in Tang et al. (2007)) for convergence of the ranking of influential continuous model parameters, a sample size of 7000 includes 46% and 12% of all possible model input combinations in the Schwechat and the Raab case studies, respectively. All sampled combinations were assembled to executable SWAT models. Daily discharge and 15 daily NO − 3 -N loads at the outlets of the Schwechat and the Raab catchments were simulated for the period from 2071 to 2100. The analysis of discharge and NO − 3 -N loads follows two main goals i) to identify the dominant controls on the simulation of discharge and NO − 3 -N loads in the two case studies and ii) to assess how the considered inputs control the simulation of discharge and NO − 3 -N loads.

Global sensitivity analysis
To measure the relative importance of the developed model input scenarios, the model setup and the parametrization on the simulation of daily discharge and daily NO − 3 -N loads, we employed GSA using the PAWN sensitivity index (Pianosi and Wagener, 2015). PAWN employs the empirical cumulative distribution function (CDF) of a target variable to infer the model input influence (Pianosi and Wagener, 2015). PAWN is moment-independent and was found to be a robust measure for sensitivity of 5 non-symmetrically distributed outputs of environmental models (Pianosi and Wagener, 2015;Zadeh et al., 2017).
PAWN expresses the sensitivity of a target variable y to a model input x by computing a distance measure between the unconditional CDF F y (y) (where all model inputs are perturbed) and the conditional CDF F (y|xi) (y) (where the model input of interest is fixed and all others are perturbed). Pianosi and Wagener (2015) proposed is the Kolmogorov-Smirnov test statistics as a distance measure. The distance KS j (x j i ) between the CDFs for the model input x i fixed at a value x i = x j i is defined as: To assess the overall sensitivity considering all fixed values of x i , the values of KS j (x j i ) are summarized for all j sampling points. A summary statistics (Pianosi and Wagener (2015) suggested e.g. median or maximum) is applied to compute the PAWN index T i for the model input x i . The model inputs that are analyzed in this study strongly differ in their numbers of discrete realizations. Further, the distribution of the resulting Kolmogorov Smirnov distances can be highly skewed (e.g. the 15 majority of discrete realizations has a low impact, while a few realizations strongly influence the simulation). Therefore, the significance of an average sensitivity of a target variable y i to a model input x i is questionable. In a setting where the strongest impact of a model input x i on a target variable y i is of major interest the application of a maximum statistics is beneficial.
Hence, the PAWN sensitivity index is defined here as: 20 The values x i = x 1 i , . . . , x j i , . . . , x ni i are the n i discrete realizations of the input x i . The resulting PAWN sensitivity index varies between 0 and 1 where a lower value of T i implies a lower influence of the input x i on the target variable y. (2015) introduced the PAWN sensitivity method using a specifically tailored sampling design to infer the PAWN indices T i for continuous model inputs x i . The proposed sampling scheme suggests to draw N c conditional samples at n randomly sampled points of each influencing variable x i , where x i is fixed at a value x i = x j i while all others are perturbed. 25 Recently, Pianosi and Wagener (2018) extended the applicability of the PAWN sensitivity method to estimate T i from a generic random sample of continuous model inputs. To approximate T i the generic sample N is split into n segments along each model input dimension resulting in conditional samples N c with an approximate size of N/n. We employed the proposed updated sampling strategy and adapted it for the use with discrete model inputs. A sample of the size N was drawn. For each model input combination every model input was sampled randomly from its discrete realizations. To infer KS j (x i ) for all discrete values x j i of a model input x i the sample N was split into subsets for all n i discrete values, resulting in subsets of the size N/n i on average. It is important to consider, that the subset size depends on the number of discrete values n i of a model input x i , while the subsets resulting from the sampling scheme proposed by Pianosi and Wagener (2018) had an average size of N/n for all model inputs x i .

Pianosi and Wagener
To account for the effect of different numbers of discrete realizations of the analyzed inputs, but also to assess whether the 5 number of drawn samples of input combinations (N = 7000) was sufficient to perform a GSA with PAWN, confidence intervals were calculated for the PAWN indices applying bootstrapping (Hinkley, 1988;Efron, 1987) using the R package 'boot' (Canty and Ripley, 2017). To calculate the bootstrap mean and the 95% confidence intervals, 1000 bootstrap replicates were drawn (as demonstrated in Sarrazin et al. (2016)).
Signature measures of discharge and NO − 3 -N loads were used as target variables y. Signature measures are measures that 10 describe specific characteristics of simulated time series (Euser et al., 2013) (in this case of daily mean discharge and daily sums of NO − 3 -N loads). We calculated quantile values (0.01, 0.05, 0.20, 0.70, 0.95, and 0.99) of daily discharge and daily NO − 3 -N loads, long-term mean discharges and long-term mean sums of NO − 3 -N loads on an annual basis and for the meteorological seasons spring, summer, autumn, and winter, and mean NO − 3 -N concentrations for different ranges of discharge quantiles (very high discharge (above 0.95 quantile), high discharge (0.95 to 0.70 quantile), medium discharge (0.70 to 0.20 quantile), 15 low discharge (0.20 to 0.05 quantile), and very low discharge (below 0.05 quantile)).  respectively. The climate scenarios were the most relevant inputs for the simulation of seasonal mean discharges and seasonal sums of NO − 3 -N loads. For the simulation of low discharge quantiles (large daily discharges) climate scenarios showed the highest relevance. For the simulation of low discharges however, the importance of the climate scenarios decreases, while the model parametrization becomes more relevant (from dark green to light green in Fig. 2). The mean PAWN indices of climate scenarios drop from 0.74 to 0.47 in the Schwechat catchment and from 0.82 to 0.51 for the simulation of lower discharges, while the mean PAWN indices for the model parametrization show respective increases from 0.43 to 0.87 and 0.44 to 0.80.

Visual analysis of the simulation uncertainties
In general, the model parametrization was highly influential for all calculated signature measures and is comparable to that of the climate scenarios, with mean PAWN indices ranging between 0.43 to 0.90 in the Schwechat and 0.36 to 0.80 in the 5 Raab (fifth row Fig. 2). Particularly, for the simulation of NO − 3 -N concentrations the model parametrization was the most dominant control of the variable simulated. In contrast to the large impact of the model parametrization, the relevance of the model setup was much lower for the simulation of discharge and NO − 3 -N loads and concentrations. Overall, values of the PAWN index for the choice of the model setup did not exceed 0.37, and were much smaller (two to five times) compared to the model parametrization. The model setups yielded insignificantly low PAWN indices for the majority of signature measures 10 with values below 0.1 in the Raab case study (2.5 % CI almost 0 for many signature measures), indicating that the model setup had a low influence on most of the analyzed processes. Only for high discharges and large NO − 3 -N loads a mean value for the PAWN index above 0.1 is visible.

Analysis of the simulation uncertainties of discharge and NO − 3 -N loads
Using all 7000 combinations of land use, point source emissions, climate, model setups, and model parametrizations, the 15 simulated discharges and NO − 3 -N loads deviated by up to 350% (grey bands in Fig. 3) from the simulations of discharge and NO − 3 -N loads in the reference period 2003 to 2015 (dashed line in Fig. 3). In the Schwechat (left column in Fig. 3) wider uncertainty bands are visible for the spring and early summer months. The results for the Raab catchment (right column) show wider uncertainty bands emerged for summer as well as for winter/early spring. A notable difference between the two case studies is how the simulations of long term monthly discharges and NO − 3 -N loads in the reference period compare to the ranges 20 of future simulations. While the majority of model combinations for the Schwechat simulated larger discharges and NO − 3 -N loads for all months in the future, for the Raab catchment the simulations of discharge and especially NO − 3 -N loads are lower in comparison to the reference period.
The analyses of the uncertainty bands with respect to the implemented land use scenarios and the point source scenarios fully confirm the results from the SA (Fig. 4). The attributed uncertainty bands for the two land use scenarios almost entirely    Table A5. The corresponding population growth scenarios (Pop. in the legend) are listed in Table A6 and the corresponding industrial emission scenarios in the Raab catchment (Ind. in the legend) are listed in Table A7.
With the GSA we identified the climate scenarios to have a great influence on all signature measures of the simulated variables. Attributing the uncertainty bands to the individual GCM-RCM combinations unveils diverse outcomes for the future flow regime, the distribution and amplitude of monthly NO − 3 -N loads, as well as the appearance of high and low discharges and NO − 3 -N loads (Fig. 5). A visual analysis of the separated uncertainty bands identifies that the deviations of the mean annual precipitation of the GCM-RCM combinations have a strong impact on the simulation of discharge and NO − 3 -N loads.

5
In comparison to the reference period (dashed line), wetter future climate scenarios (blue) simulated larger discharge and NO − 3 -N loads, while drier future conditions lead to a drastic reduction in discharge and NO − 3 -N loads. These findings further imply that NO − 3 -N applied in fertilizers will remain in the upper soil layers and be transformed (mineralized or immobilized or denitrified) instead of being transported to the receiving waters. A comparison of the NO − 3 -N budgets of simulations with dry and wet climate scenarios for the Raab shows a difference of up to +27% of NO − 3 -N accumulated in the soil, as well as a 10 decrease of 43% and 38% in NO − 3 -N yield in the fast and slow runoff, respectively. Half of the 22 implemented GCM-RCM combinations simulated an increase of more than 75 mm (dark blue) and for two GCM-RCM combinations, an increase of more than 25 mm (light blue) of precipitation for the Schwechat catchment was simulated. In contrast, for the Raab nine and four GCM-RCM combinations simulated a decrease in precipitation of more than 75 mm (dark red) and 25 mm (light red), respectively. Consequently, a decrease in discharge and NO − 3 -N loads due to a 15 decrease in precipitation is pronounced in the Raab catchment, while the majority of simulations of the Schwechat catchment show an increase in discharge and NO − 3 -N loads. While a grouping of the individual climate scenarios with respect to their temperature deviations shows a more indefinite picture, all climate scenarios simulated an increase in temperature. Nevertheless, the expectation that an increase in annual mean temperature increases evapotranspiration and thus reduces discharge and NO − 3 -N loads is not met in Fig. 6. A clear separation 20 of warmer and cooler climate scenarios, as it is observable for precipitation is not the case with temperature. Consequently, the differences in precipitation predominantly account for the influence of the climate scenarios, rather than the differences in temperature.
Although the influence of the model setups was much lower compared to the influence of the climate scenarios or the model parametrization, the analysis of the uncertainty bands for the different model setups provides interesting insights (Fig. 7). 25 The uncertainty bands do overlap to a great extent, which confirms a low impact of the use of different model setups in the simulation of discharge and NO − 3 -N loads. Noteworthy is, that model setups that use the full set of HRUs agree much stronger in their simulations compared to the model setups where the number of HRUs was reduced. The difference between the full HRU and the reduced HRU model setups is distinct in the Schwechat case study. The uncertainty bands of the two full HRU model setups almost completely overlap, although their numbers of subbasins are different (4 and 14 subbasins). The two    Large values of CNOP_till and small values of SOL_AWC reduce the water retention capacity and increase the amplitude of medium and low discharges (third row in Fig. 8). A similar but inverse behavior is visible with medium NO − 3 -N loads (last row in Fig. 8), where a higher water retention results in an increase of NO − 3 -N loads. For the long-term monthly mean discharges and sums of NO − 3 -N loads two effects are observable in Fig. 8. First, smaller values of CNOP_till and larger values of SOL_AWC decrease the upper boundary of the uncertainty bands. Second, selected model parametrizations with large values 5 of CNOP_till and small values of SOL_AWC cause considerably larger discharges in spring and a strongly reduced runoff in the autumn months in the Schwechat case study.

What can we as modelers learn from such analysis
The illustrated case studies emphasized the necessity to characterize, identify and explicitly communicate the uncertainties 10 in a modeling chain, particularly for future simulations of environmental variables where large uncertainties are inherent in several modeling inputs. While the sensitivity analysis of signature measures related to discharge, NO − 3 -N loads and NO − 3 -N concentrations provided a comprehensive overview of the dominant influencing inputs on specific modeled variables, the analysis of the uncertainty bands for the simulation of the modeled variables provided insights into which properties of the model inputs (e.g. mean annual precipitation or mean air temperature of a climate scenario) control the uncertainties and how 15 these control the simulation. The analyses allow to draw conclusions that are beneficial to consecutive steps of an impact study, for instance to refine the impact study setup and to focus on the most influential components and ultimately to reduce the uncertainties in the modeling simulation chain.
The land use scenarios showed an almost negligible impact on the simulation of discharge and NO − 3 -N loads. The discharge and the NO − 3 -N loads at the catchment are however integrated signals for the entire catchment and changes in land use may 20 have a greater importance for particular points in a catchment. Many case studies have applied the SWAT model to assess the impact of land use change on different variables of the water cycle (Wagner et al., 2017;Mehdi et al., 2015b), water quality Mehdi et al., 2015a;Teshager et al., 2016), or sediment yield (Bieger et al., 2013). Bieger et al. (2013) found very low land use change induced increases in discharge for a catchment in China. Only an assumed strong intensification of the agriculture led to a 4% increase in discharge. At the same time however, a strong increase in sediment yield of up to 450% 25 for the summer months was simulated due to the intensification of agriculture. Guse et al. (2015) also found only small changes in simulated discharge caused by future land use change in a German lowland catchment. In absolute numbers the simulated future NO − 3 -N loads showed small differences between the baseline scenario and the two applied methods of land use change presented by Guse et al. (2015). Yet, the temporal patterns in NO − 3 -N loads caused by the different approaches of changing the land use were the major observable difference. Mehdi et al. (2015b) however found that including agricultural land use 30 change into the impact assessment of a southern German watershed strongly increased the NO − 3 -N and total phosphorus loads. Teshager et al. (2016) support the findings of Mehdi et al. (2015b) and also found that corn intensive scenarios lead to an increase in discharge and significant water quality problems while an extensive scenario where mainly switchgrass is planted lead to water quality improvements under future climate change. Consequently, the low impact of land use change found in the present study seems reasonable with respect to other literature, particularly as no extreme scenarios were implemented. This does however not generally imply a low importance of land use change in environmental impact assessments. Land use change or changes in the management can be the most relevant input, particularly when strong future changes, such as possible bans of emittents are considered (Honti et al., 2017).

5
Industrial emitters were the main cause for the impact of point sources on medium to low NO − 3 -N loads. The future scenarios of the development of industrial emitters were however highly uncertain. The developed scenarios are based on expert knowledge. Yet, there is no reliable basis available on status of the industrial emitters by the end of the century. Therefore, the developed scenarios should be noted as feasible futures, rather than e.g. politically realizable futures (Godet and Roubelat, 1996). To set a feasible range as boundaries for the future development of industrial emitters can lead to an overestimation 10 of their impact in comparison to other influencing variables. Nevertheless, the visualization of the NO − 3 -N FDC of the Raab case study highlights the effect of the industrial emissions for medium and small NO − 3 -N loads. Large NO − 3 -N loads however, are hardly affected by the implemented scenarios, indicating that large NO − 3 -N emissions are mainly driven by agricultural activities.
The selection of climate scenarios had a strong influence on the simulation of discharge and NO − 3 -N loads in both case 15 studies. The analysis of the uncertainties bands identified the differences in precipitation between the GCM-RCM combinations as being the main control, while the differences in air temperature had a low impact on the simulation outcome. This finding stands in contrast to other studies. Milly and Dunne (2011) and Sheffield et al. (2012) for example, identified empirical approaches for the calculation of evapotranspiration as the main source for overestimation of the climate's influence on hydrological processes, particularly when evapotranspiration is a function of air temperature (Clark et al., 2016;Shaw and 20 Riha, 2011;Roderick et al., 2014). In the climate scenarios used in this study, the impact of large differences in mean annual precipitation on the simulated outputs exceeded the impact of the differences in air temperature.
The effect of the model setup, with different watershed subdivisions, on the simulation of discharge or water quality variables has been investigated in various studies (e.g. Jha et al., 2004;Momm et al., 2017;Pignotti et al., 2017). Jha et al. (2004) emphasize the greater impact of changes during the HRU definition over the defined number of subbasins, as a consequent 25 change in the distribution of land use, soil, or topography strongly affect runoff and the nutrient budget in a catchment. The analysis of the uncertainties bands with respect to the different model setups clearly confirmed the study by Jha et al. (2004), especially in the case of the Schwechat. Nevertheless, the impact of the model setup was lower than the effect of the model parametrization by a factor of up to five in the Schwechat study and up to eight in the Raab case study. Yet, the model setup strongly affects the computation time. In the present case, where aggregated discharge and NO − 3 -N loads at the catchment 30 outlets were the variables of interest a strong focus on the model parametrization is of higher priority than the spatial distribution of the model setup. Therefore, to maintain short computation times (and at the same time to maintain the distributions of land use, soil, or topography) a model setup with a low number of subbasins without any reduction of the number of HRUs is beneficial.
The impact of parameter non-uniqueness on the simulation of hydrological and water quality variables has been demonstrated previously (e.g.; Wilby, 2005;Mehdi et al., 2018). The importance of the model parametrization for the simulation of discharge and NO − 3 -N loads was confirmed in the present study as well. Large sensitivities of all signature measures of discharge and NO − 3 -N loads to the different model parametrizations were identified . Although all selected parameter sets represented historical observations of discharge and NO − 3 -N loads with a certain goodness of fit (based on defined objective 5 criteria), the colored grouping of the uncertainty bands illustrated that the selected model parameter sets control the simulation of future discharge and NO − 3 -N loads in different ways. Thus, the large impact of the model parametrization and the distinctive patterns identified in the uncertainty bands suggest a great potential to further refine the model parametrization and consequently reduce simulation uncertainties with a more intensive model calibration. Additional information on the time series of observations can help to constrain the model parameters and adequately describe the relevant processes (e.g. Hrachowitz et al., 10 2014;Pfannerstill et al., 2017).

How to attribute subjectivity inherent in the scenarios
Scenarios always reflect subjective assumptions made by the modeler. Assumptions that are made in the scenario development however, can strongly influence a simulation and thus affects a comparison of different model inputs and their impacts on the simulation. All steps in a scenario development involve subjective assumptions and can lack plausibility (Mahmoud et al., 15 2009;van Vuuren et al., 2012), regardless of whether the process involves expert knowledge, the input of stakeholders in an participatory process, or an exploratory approach that extrapolates trends, these practices potentially introduce uncertainties in the definition of scenarios. Technical aspects such as how the scenario is represented in the model are also strongly biased by the modelers decision and represent an additional source of uncertainty (Mahmoud et al., 2009). The communication of the potential uncertainties inherent in the developed scenarios and the boundaries of the explanatory power of an scenario ensemble 20 is essential for the integrity of any impact study (Mahmoud et al., 2009;Jones et al., 2014).
In the present study, several assumptions were made in the development of scenarios that are highly subjective, such as the extrapolated gradient of future land use changes, the drastic changes in future industrial emissions, and also the selection of objective criteria that define a behavioral SWAT model setup. Scenarios must cover a broad range of possible futures and have to be adequately represented in the model setup. An explicit delineation of the implemented scenarios and their limitations is 25 essential to clearly illustrate the limitations of an impact study's conclusions. An immanent risk in any impact study is that the model representation of a future change, or the uncertainties in a model input fail to reproduce the response of a simulated variable that would have taken place in the real environmental system. Hence, a detailed analysis of the simulation uncertainties perfectly complements a SA to identify possible shortcomings in the study setup. Attributing the uncertainty bands resulting from the simulation of an environmental variable to individual model inputs prove to be a useful visual analysis tool that gives 30 the power to illustrate the uncertainties in a transparent way. Furthermore, the colored differentiation provides a visual guidance to judge the impacts of different implemented scenarios.

Sensitivity analysis or hydrologic storylines
The presented approach implements large samples combining scenarios for different model inputs and different model setups and parametrizations in a GSA to identify the dominant contributors of uncertainties in the simulated outputs. The utilization of SA with large sample sizes however, raises the following issues: i) compared to a standard approach to perform an impact assessment, where a few different future scenarios are implemented into a model, the computational demand of a GSA requiring 5 hundreds or thousands of model executions is larger by several orders of magnitude. Thus, a practical implementation of the presented procedure in impact studies is questionable and a strong cooperation between research and the practitioners is essential. ii) scenarios of different model inputs are often interrelated (Mahmoud et al., 2009). A change in one model input therefore for example expects the change of another model input into one direction and makes a change into another direction unlikely. While the implementation of input dependencies, althouh challenging is feasible for continuous model 10 inputs, for instance by a transformation of the input space (e.g., Tarantola and Mara, 2017;Mara and Tarantola, 2012), or the determination of input distribution functions (Hart and Gremaud, 2018), the dependencies of composite model inputs are usually difficult to express mathematically. To identify the dependencies between composite model inputs, expert knowledge is required to properly constrain the model input combinations and therefore complicates the implementation in approaches, such as the presented one. 15 Clark et al. (2016) therefore suggest to identify consistent hydrologic story lines that result in least severe, most likely, and most severe responses of the modeled system. Such an approach would tremendously reduce the number of necessary model evaluations, but also establish consistency between the considered influencing variables. Nevertheless, the feasible combinations of influencing variables that lead to extreme or likely responses of the modeled system are hardly known a priori. Consequently, a sensitivity analysis with a constrained sampling space, to avoid infeasible combinations of influencing 20 variables might be a pragmatic compromise.

Conclusions
In this study we utilized methods for GSA in environmental impact studies to identify the dominant sources of uncertainties for the simulation of environmental variables under future changing conditions. In two Austrian case studies for the rivers Schwechat and Raab, we simulated the river discharge and the NO − 3 -N loads from the catchments under the condition of future 25 changes in climate, land use, and emissions from urban and industrial point sources implementing different SWAT model setups with various model parametrizations.
Both case studies identified climate change and the model parametrization to be the most important (influential) model inputs for the simulation of discharge and NO − 3 -N loads, based on performing a GSA and on the resulting analysis of signature measures of discharge and NO − 3 -N loads (quantiles of discharge and NO − 3 -N loads, seasonal mean discharge and seasonal 30 sums of NO − 3 -N loads and NO − 3 -N concentrations for discharge quantiles). The impact of the model setup on simulated variables of discharge and NO − 3 -N loads was found to be considerably lower than the impact of the model parametrization for the Schwechat and even more distinct for the Raab. The impact of the implemented scenarios for land use and municipal point source emissions were negligible for all analyzed signature measures. Because of a large leather industry in the Raab catchment, the future development of industrial emission in the Raab catchment was found to be relevant for low NO − 3 -N loads and NO − 3 -N concentrations during low discharge. Accompanying the GSA, a detailed analysis of the simulation uncertainties provided additional insights on how the uncertainties in the model inputs control simulated discharge and NO − 3 -N loads. The visualizations we developed supported the 5 identification of the relevant properties of the model inputs that control the simulation uncertainties and provide insight how individual realizations of a model input can affect the simulations. In the climate simulations, we found the precipitation to dominate the simulation outputs, rather than changes in air temperature. Although the impact of the model setup on the simulation of discharge and NO − 3 -N loads was low, the visual analysis of the uncertainty bands illustrated that the HRU definition is an important step in the model setup. The use of the full set of HRUs was identified as the preferred setup in the two case    Table A2.    Table A8. GCM-RCM combinations implemented in the study with their long-term mean annual precipitation sums and long-term mean annual temperatures for the Schwechat and the Raab. Honti, M., Schuwirth, N., Rieckermann, J., and Stamm, C.: Can integrative catchment management mitigate future water quality issues caused by climate change and socio-economic development?, Hydrology and Earth System Sciences, 21, 1593-1609,