Unraveling the time-dependent relevance of input model uncertainties for a lumped hydrologic model of a pre-alpine karst system

Uncertainties in hydrologic model outputs can arise for many reasons such as structural, parametric and input uncertainty. Identification of the sources of uncertainties and the quantification of their impacts on model results are important to appropriately reproduce hydrodynamic processes in karst aquifers and to support decision-making. The present study investigates the time-dependent relevance of model input uncertainties, defined as the conceptual uncertainties affecting the representation and parameterization of processes relevant for groundwater recharge, i.e. interception, evapotranspiration and snow dynamic, on the lumped karst model LuKARS. A total of nine different models are applied, three to compute interception (DVWK, Gash and Liu), three to compute evapotranspiration (Thornthwaite, Hamon and Oudin) and three to compute snow processes (Martinec, Girons Lopez and Magnusson). All the input model combinations are tested for the case study of the Kerschbaum spring in Austria. The model parameters are kept constant for all combinations. While parametric uncertainties computed for the same model in previous studies do not show pronounced temporal variations, the results of the present work show that input uncertainties are seasonally varying. Moreover, the input uncertainties of evapotranspiration and snowmelt are higher than the interception uncertainties. The results show that the importance of a specific process for groundwater recharge can be estimated from the respective input uncertainties. These findings have practical implications as they can guide researchers to obtain relevant field data to improve the representation of different processes in lumped parameter models and to support model calibration.


Introduction
Hydrologic models serve as important tools for the assessment of dominant hydrodynamic processes in karst systems (Hartmann et al. 2013;Sivelle et al. 2019). In those models, subsurface heterogeneity and the resulting complex hydrodynamic processes typical for karst aquifers are often represented in a simplified way (Fleury et al. 2007;Guinot et al. 2015;Sivelle et al. 2021;Sivelle and Jourde 2020). The assessment of the reliability of a model output is therefore an important step towards an improved description of the karst system (Hartmann et al. 2014a). This assessment is usually done by uncertainty quantification techniques, which investigate the likelihood of a model outcome while considering the unknowns in a hydrologic model (Sarrazin et al. 2018;Teixeira Parente et al. 2019). These unknowns arise from different sources of uncertainties, i.e. structural (Fandel et al. 2020;Henson et al. 2018), parametric (Mazzilli et al. 2012;Moussu et al. 2011) and input uncertainties (Liu et al. 2018;Nerantzaki et al. 2020).
Structural uncertainties evolve from the simplifications required while creating a conceptual model of a real-world system (Gupta and Govindaraju 2019;Rojas et al. 2008). This conceptualization often neglects certain parts of the natural system due to a lack of knowledge, which can lead to an underrepresentation of important hydrodynamic processes (Butts et al. 2004;Lee et al. 2011). Parametric uncertainties arise from the fact that the exact values of model parameters, such as discharge coefficients and storage thresholds, are often not known (Ahmadi et al. 2019;Hu et al. 2019). This is particularly true for lumped conceptual models, whose parameters cannot often be constrained by physical field experiments (Wagener et al. 2003). Hence, for each parameter, a reasonable parameter range needs to be defined in which the true parameter value is located (Seibert 1997). Finally, input uncertainties exist due to missing and/or uncertain input data (Breinl 2016;McMillan et al. 2012) as well as due to simplifications of the processes that finally represent the model input, e.g. the groundwater recharge (Kavetski et al. 2006;Patil et al. 2011;Vrugt et al. 2008).
More recent studies highlighted that groundwater recharge in systems with strong subsurface heterogeneities such as karst systems, exhibits a high sensitivity to changes in climatic forcings (Hartmann et al. 2017). In the specific case of prealpine karst catchments, these forcings controlling groundwater recharge are interception, evapotranspiration and snowmelt processes. Garrigues et al. (2015) and Sarrazin et al. (2018) showed that the sensitivity of groundwater recharge with respect to vegetation related processes, i.e. interception and evapotranspiration, mainly results from the spatial variability of soil properties. Moreover, Ollivier et al. (2021) underline that this sensitivity is further related to often missing information about spatially distributed and vegetation-dependent evapotranspiration dynamics. In cases where snowmelt represents a controlling factor in the water balance of karst areas, Doummar et al. (2018) showed that groundwater recharge estimations are most sensitive to temperature variations. That is mainly due to the importance temperature has for the timing of snow accumulation and melt and the resulting control on a spring's discharge behavior (Liu et al. 2021).
As it is difficult to measure interception, evapotranspiration and snowmelt, additional models are often applied to compute input time series for hydrologic models (Hartmann et al. 2014b;Mazzilli et al. 2012;Ollivier et al. 2020). Interception can be modeled using mechanistic (Gash et al. 1995;Liu 2001) and stochastic modeling approaches (Calder 1996;Hall 2003). Data demanding energy balance methods (Colaizzi et al. 2012;Penman 1948) or simple temperature-based parametrizations (Oudin et al. 2005;Thornthwaite 1948) provide evapotranspiration time series. Snow processes, which were recently shown to play a major role for groundwater recharge in pre-alpine and alpine areas (Jódar et al. 2020;Lucianetti et al. 2020), can be modeled using energy balance methods (Herrero et al. 2009;Marks et al. 1999) or simpler degree-day-factor methods (Girons Lopez et al. 2020;Martinec 1960). Reliable time series of interception, evapotranspiration and snowmelt are a prerequisite for a proper description of the water balance for distributed, semidistributed and conceptual models.
An example for a conceptual model applied to a pre-alpine karst system is the LuKARS model, which was developed by Bittner et al. (2018) for the Kerschbaum springshed in Austria. Given the natural characteristics of the springshed (forested catchment, elevation between 415 and 969 m asl, annual mean temperature of 8°C), evapotranspiration, interception and snowmelt representations are expected to have an important influence on the modeled spring discharge (Bittner et al. 2018). However, since no direct measurements for these input data are available, simple algorithms have to be applied for computing the input time series for the LuKARS model (Bittner et al. 2020a).
While previous studies investigated the parametric uncertainties of LuKARS for the Kerschbaum springshed (Teixeira Parente et al. 2019), the presented article aims to investigate how much the input uncertainties affect model predictions. The hypotheses, which the authors want to test in this study, are, first, that the input uncertainties can vary seasonally, then, that it is possible to derive the specific importance of a single process, e.g. snowmelt, for groundwater recharge from its related uncertainties in the model output. This is of particular importance, as it serves as a practical example which is beneficial to guide researchers and decision-makers in favoring field experiments and data collection differently during different seasons to improve the output of a karst aquifer model. To study the uncertainty propagating to the spring discharge, three different modeling approaches are applied for each of the considered hydrological processes, i.e. interception, evapotranspiration and snow processes. In particular, the methods of DVWK (1996), Gash et al. (1995) and Liu (2001) are applied to compute interception, the methods of Hamon (1961), Oudin et al. (2005) and Thornthwaite (1948) to calculate evapotranspiration and the methods of Girons Lopez et al. (2020), Magnusson et al. (2014) and Martinec (1960) to model snowmelt and accumulation. The selection of lumped approaches is driven by data availability in the study area and in particular by the lack of radiation data. Then, all possible model combinations are run varying the parameters of the input models and by using the sampling algorithm of the Fourier Amplitude Sensitivity Test (FAST; Pianosi et al. 2015). Finally, the study investigates the impact of the input model parametrization on the LuKARS model output and how this changes over time. To conclude, the input uncertainties are compared with the parametric uncertainties, which were computed in an earlier study (Teixeira Parente et al. 2019).

The study area
The Kerschbaum springshed is located close to the city of Waidhofen a.d. Ybbs (Fig. 1a), about 100 km west of the city of Vienna (Austria; Fig. 1b). The study site covers an area of 2.5 km 2 . This pre-alpine recharge area forms part of the eastern foothills of the Northern Calcareous Alps and is dominated by a lithologic sequence of dolomitic basement rocks (Fig. 1c). The study area shows karst features such as springs, dry valleys and caves. Due to the absence of significant sinkholes, the groundwater recharge can be assumed barely influenced by pointinfiltration processes. Moreover, according to the study of Narany et al. (2019), the Kerschbaum springshed is characterized by a deep karstified groundwater system with a wellconnected network of fractures and conduits.
The land cover is dominated by beech forests. Its spring provides a mean discharge of 34 L/s to the regional water supply and shows a quick reaction time to precipitation and snowmelt events of 1 day . Figure 2 shows the available precipitation, temperature, snow depth and discharge time series for the period from 1 January 2006 to 31 December 2007. These time series were measured at the weather station whose location is shown in Fig. 1a. For more information about the study area, the interested reader could refer to the publication of Narany et al. (2019).

Methodology
The paper briefly describes the lumped karst hydrological model LuKARS of the Kerschbaum springshed. For more information about the model the reader could refer to the publication of Bittner et al. (2020b). The paper then focuses on the description of both commonly applied and recently proposed parameterizations for interception, evapotranspiration and snow processes. Since precipitation and air temperature are the only meteorological parameters measured in the study area, temperature-based methods are used to compute evapotranspiration and snowmelt. Finally, the method used to quantify the uncertainty of all investigated model combinations is described.

The LuKARS model
The LuKARS model is a lumped parameter model developed by Bittner et al. (2018) that considers the dominant hydrotopes in a recharge area as distinct response units. Hydrotopes are defined as landscape units with similar soil and land use characteristics (Arnold et al. 1998). Each hydrotope is characterized by a specific retention capacity. As Fig. 3a shows, shallow and coarse-textured soils lead to low soil storage and high quickflow intensity. In contrast, thick and fine-textured soils lead to high soil storage and low quick flow intensity. The conceptual model considers the hydrotopes to represent the vadose zone (soil-epikarst-infiltration zone) and to be directly connected to the saturated zone, which consists of a single linear storage recharged by each hydrotope independently. The duality of flow behavior is implemented by considering for each hydrotope both the fast flow component through the conduits and the slow diffusive discharge in the matrix. Each hydrotope simulates three different types of flow, i.e. the quickflow (Q hyd ), the matrix infiltration (Q is ) that feeds the baseflow storage (B), and the secondary spring discharge (Q sec ; Fig. 3b). Q hyd represents the discharge that is directly moved to the outlet of the catchment through preferential flow paths such as subsurface conduits, and factors that are  (GBA 2021) including the elevations of the Kerschbaum spring and the summit of the Glashüttenberg responsible for the fast reaction of the spring discharge to rainfall and snowmelt events. Q hyd is implemented considering the hysteretic behavior of the soil-epikarst system that starts after a constant hydrotope specific storage value (E max ) was exceeded and stops after a lower constant threshold (E min ) was reached. Q is is the water that infiltrates into the lower reservoir B (Fig. 3) and, thus, represents the groundwater recharge. Q sec is the flow that discharges outside the investigated recharge area and is activated only when the threshold for secondary spring discharge (E sec ) was exceeded. Q is , Q sec and Q b (the baseflow) are implemented using linear transfer functions. Finally, Q tot is the discharge at the spring. The mathematical equations and a graphical user interface for the model are provided in Bittner et al. (2018) and (2020a), respectively.

Interception
The approach applied by Bittner et al. (2018) was based on the percentages for interception of beech forest stands proposed by DVWK (1996). This study further considers the methods proposed by Gash et al. (1995) and Liu (2001). DVWK (1996) suggests that 11% of precipitation is intercepted from beeches in the winter season (21 December, d w ), whereas 17% of precipitation is intercepted in the summer season (21 June, d s ). A linear interpolation is applied between these values following Eqs. (1) and (2), which compute daily time series of interception I (mm/day) for the time between 21 December and 21 June and the time if I < 5 mm=day; The approach by Gash et al. (1995) is based on the calculation of a gross rainfall that is needed to saturate the canopy, i.e. P′ g (mm). P′ g is calculated using Eq. (3).
where C m is the stand storage capacity (mm), ER (−) is the ratio between the mean evaporative rate E and the mean rainfall rate of the event for saturated canopy conditions R. The parameter p represents the free throughfall coefficient (−). For a daily time step, if the precipitation P (mm/day) is larger than P′ g (mm/day), I (mm/day) can be calculated with Eq. (4).
The method proposed by Liu (2001) requires the definition of the same parameters as in the method of Gash et al. (1995). However, instead of differentiating between the cases in which precipitation P is greater or smaller than the gross rainfall that is needed to saturate the canopy (P′ g ), the exponential function in Eq. (6) is defined to compute daily interception amounts (I).
Liu (2001) investigated the parameter sensitivities of Eqs. (4) and (6) and showed that an overestimation of either C m or ER results in an overestimation of interception, whereas large Fig. 3 The conceptual modeling approach of LuKARS. a Conceptual representation of the four implemented hydrotopes. Hyd 1 indicates the dolomite quarries with no groundwater recharge and the dominance of surface runoff (SF). Hyd 2 and Hyd 4 represent coarse-textured and fine-textured soils, respectively. b The bucket-type model implementation of dominant hydrotopes. Q sec is the secondary spring discharge, Q hyd the quickflow, Q is the matrix infiltration feeding the baseflow storage B, Q b the baseflow and Q tot the discharge at the spring. E sec , E max and E min are thresholds storage values regulating the activation of the discharge components values of p cause an underestimation of interception. Moreover, the parameter sensitivities depend on the magnitude of a considered rainfall event and the type of canopy. For example, ER is most sensitive in areas dominated by intense rainfall events, whereas ER and C m are most sensitive in areas characterized by small rainfall events. The parameter p is sensitive when the models are applied to areas with small rainfall events and open canopies.

Evapotranspiration
The Thornthwaite (1948) evapotranspiration model was used to calculate the potential evapotranspiration (ET pot ) in the original LuKARS model of the Kerschbaum spring (Bittner et al. (2018). In this work, the simulation approaches proposed by Hamon (1961) and Oudin et al. (2005) are additionally applied. It is important to note that Bittner et al. (2018) and Teixeira Parente et al. (2019) used ET pot as actual evapotranspiration (ET act ), since the results obtained for the annual ET pot losses were in good agreement with ET act computed in previous studies for the same study area (Markart et al. 2006). Specifically, for the years -2007, Bittner et al. (2018 calculated a total of 45% ET pot , which is in good agreement with the 43% ET act presented in Markart et al. (2006). In the presented study, ET act is hence considered to be equal to the ET pot to be able to compare the different model configurations under the same conditions. Thornthwaite (1948) provides estimates for monthly ET pot and assumes that once the mean monthly temperature becomes larger than 0°C, ET pot becomes 0. In this method, the hours of daylight are assumed to be 12 and each month is 30 days long. Then, ET pot (mm/month) is calculated as shown in Eq. (7), where k TH (mm/month) is a proportionality constant, T mean (°C) the mean monthly temperature and H (°C) is the heat index defined in Eq. (8), with r an exponent given by Eq. (9).
r ¼ 6:75e −7 H 3 þ 7:71 e −5 H 2 þ 1:792 e −2 H þ 0:49239ð9Þ As proposed by Bittner et al. (2018), the monthly ET pot are divided by the number of days and the resulting daily ET pot are to be representative for the 15th day of a month. In order to obtain daily ET pot (mm/day) from this method, a linear interpolation is applied between these representative ET pot values. Oudin et al. (2005) derived an empirical equation to estimate daily ET pot (mm/day) as input for lumped rainfall-runoff models. They tested various ET pot modeling approaches for numerous catchments in France, Australia and the United States. Their goal was to identify those atmospheric variables which provide the best streamflow predictions when being used as input for ET pot models. The equation they derived is shown in Eq. (10), where H O is the extraterrestrial solar radiation [MJ/(m 2 /day)], k OU [m 3 /kg/(1,000 MJ 2°C )] is a proportionality constant, T mean is the mean daily temp (°C), and 0.408 is an approximation for the latent heat flux (MJ/kg). It is important to note that ET pot is 0 (mm/day) if T mean ≤ 5 (°C). The third method applied to calculate ET pot was proposed by Hamon (1961), who derived a simple procedure to be used in water balance estimations. The goal was to use readily available data, i.e. daily air temperature (T mean ), for ET pot estimations. The derived methodology is based on the saturated water vapor concentration at T mean , i.e. e 0 (T mean ) (kPa) and is expressed by Eq. (11), where k HA (mm/day) is a proportionality constant and N is the maximum number of daylight hours. e 0 (T mean ) was approximated using the modified Magnus equation proposed by Alduchov and Eskridge (1996) that is shown in Eq. (12).
Snowmelt Bittner et al. (2018) used the method proposed by Martinec (1960) to model snowmelt and storage for the Kerschbaum spring recharge area. This study additionally considers the methods described by Girons Lopez et al. (2020) and Magnusson et al. (2014). All snow routines considered in this framework assume that the energy available for snowmelt is proportional to air temperature. This means that below a certain threshold temperature T T , precipitation falls as snow, whereas rainfall occurs for temperatures above this threshold. The proportionality of snowmelt (M) is controlled by the degree-day factor C 0 [mm/(day°C)] and the daily mean temperature T mean (°C). Moreover, all the considered snow models neglect sublimation processes, which is often the case in degree-day methods. The degree-day method proposed by Martinec (1960) simulates M (mm/day) using Eq. (13).
While Martinec (1960) considers C 0 to be constant, Braun and Renner (1992) argue that this factor should be changing over time, since environmental conditions, e.g. solar inclination and snow albedo, vary seasonally. Girons Lopez et al. (2020) describe a seasonally varying degree-day factor, i.e. C 0,n , based on a sine function. The intensity by which C 0 varies is controlled by an amplitude factor C 0,a [mm/(day°C)]. Then, C 0,n is computed as shown in Eq. (14), where n is the time (day).
Finally, Magnusson et al. (2014) approach the calculation of snowmelt with the exponential function shown in Eq. (15), where M M represents a snowmelt transition (°C). This method allows for melting to occur even below freezing temperature.
According to the investigation of Girons Lopez et al. (2020), who tested a variety of modifications to different temperature-based snow routines, the most sensitive parameters in the models applied in this study are the snowmelt transition M M and the temporally varying degree-day factor C 0, n .

Parameter sampling and investigated model combinations
Appropriate parameter ranges are defined for all unknown parameters in the equations describing interception, evapotranspiration and snowmelt. The parameter ranges selected for the nine considered calibration parameters, i.e. C m , p, ER, k OU , k HA , C 0 , T T , C 0, a and M M , are shown in Table 1. The specific range of values for each parameter is based on previous studies, which are indicated in Table 1.
This study considers the model configuration from Bittner et al. (2018) as the reference model, whose results are used to evaluate the performance of all the considered model combinations. All investigated model combinations are shown in Table 2 The sampling algorithm of the Fourier Amplitude Sensitivity Test (FAST), which was developed by Cukier et al. (1978) and implemented in the SAFE toolbox (Pianosi et al. 2015), is used in this study to obtain for each model combination a set of values that covers the full parameter space of each parameter. To avoid unrealistic water budgets, the model input time  DVWK (1996), Thornthwaite (1948) and Martinec (1960), rather than defining parameter ranges and samples, this study keeps the parameters fixed and equal to those found in Bittner et al. (2018) (Table 1). The model combinations including the total number of investigated samples are summarized in Table 2. For the sake of completeness, Appendix shows the calibrated LuKARS model parameters found in Bittner et al. (2018). The hydrotope specific parameters control the behavior of the soilepikarst system, of the matrix infiltration and of the quickflow through the conduits, while the baseflow storage parameters determine the response of the saturated zone.

Results
The results of this study are presented in the following sections. First, the uncertainties of single processes in relation to the parametric uncertainties of the Kerschbaum spring LuKARS model which were computed in an earlier study (Teixeira Parente et al. 2019). Notice that this comparison is mainly qualitative, since the parameter uncertainty was estimated using only the combination DVWK -Thornthwaite -Martinec with fixed input parameters. Then, section 'Evaluation of all model combinations' focuses on the comparison of the results of all evaluated model combinations and highlights how the input uncertainties change when considering more processes to be unknown, i.e. interception, evapotranspiration and snowmelt. The analyses focus on the comparison of the interquartile range of the LuKARS model outputs as an indicator for model uncertainty. Moreover, the interquartile ranges are normalized with the observed spring    discharge to make the interpretation of input uncertainties independent of high and low flow periods, occurring during snowmelt and snow accumulation periods, respectively. Figure 4 shows the cumulative input values generated with all parameter samples for each applied algorithm. As the analyses focus on the interquartile range of model outputs, also the generated cumulative recharge values of that range are shown. Cumulated sums allow a better visualization through continuous plots. It is observed that the inputs computed with the different algorithms are well distributed around the input time series used in the original Kerschbaum spring LuKARS model. Slight deviations are only visible for ET pot in 2006, where the method of Hamon (1961) partially overshoots the inputs generated with the method of Thornthwaite (1948), and in 2007, where the Thornthwaite (1948) method overshoots the interquartile range of inputs computed with the method of Oudin et al. (2005). These deviations are related to the linear interpolation which is applied to derive daily values from the Fig. 4 The plots show the interquartile range of the cumulative input values for each applied algorithm. For comparison, the black lines highlight the inputs used in the study of Bittner et al. (2018). a The interception inputs, b the potential evapotranspiration inputs and c the potential snowmelt inputs monthly ET pot , obtained with the methodology of Thornthwaite (1948).

Model input uncertainties related to single hydrological processes
In this section, the focus is on the uncertainties related to single hydrological processes, i.e. interception, evapotranspiration and snowmelt, how these uncertainties change in different periods of the year and how they compare to the parametric uncertainties previously computed by Teixeira Parente et al. (2019). The model parameters, which were considered in the evaluation of the parametric uncertainties of the LuKARS model, are highlighted with an asterisk (*) in the Appendix and are the discharge coefficients and exponents, minimum and maximum storage capacities, and activation level of secondary spring discharge for hydrotopes Hyd 2, Hyd 3 and Hyd 4 (Teixeira Parente et al. 2019). Figure 5 shows that the interquartile ranges resulting from uncertainties in interception (model combinations Gash and Liu) are generally smaller than the evapotranspiration and parametric uncertainties. Moreover, the uncertainties related to interception do not show a distinct seasonal variation in 2006; however, a slight seasonal variation with increasing interception over the summer period can be observed in 2007. Moreover, both Gash and Liu model follow the same temporal dynamics and lead to very similar interquartile ranges, showing that the choice of the interception model is not very significant for this case study.
Regarding the interquartile ranges resulting from the use of the Oudin et al. (2005) and Hamon (1961) methods, it is noted that the uncertainties related to evapotranspiration are characterized by a clear seasonal variability. The method of Oudin et al. (2005) brings the largest difference in the normalized interquartile range, which increases over the summer seasons in 2006 from 0.03 (29 March) to 0.06 (2 August) and in 2007 from below 0.04 (21 January) to 0.1 (5 September) and decreases again in the winter seasons. Overall, the normalized interquartile range of evapotranspiration is smaller than the parametric uncertainties. Moreover, most of the time the uncertainties of evapotranspiration are higher than the snowmelt uncertainties, even over the winter period 2006-2007. The uncertainties in snowmelt are higher than the evapotranspiration uncertainties for an extended period only in the early year 2006. Also in this case, as it was observed for interception, the difference in the normalized interquartile range between the Hamon and the Oudin models is rather small, reaching a maximum value of 0.04 on 4 September 2007. In case of snow processes, the choice of the model appears to be more relevant than for evapotranspiration and interception. On the one hand, the mean of the differences in the normalized interquartile range between both ET models, i.e. 0.011, is higher than the mean of the normalized interquartile range differences between the snow models, i.e. 0.008. On the other hand, the maximum difference in the normalized interquartile range is identified between the Girons Lopez and Magnusson models, i.e. is 0.06 on 21 March 2006. This is reasonable, since snow processes do not play a role over the whole time of a year, whereas ET does.

Evaluation of all model combinations
The minimum and maximum percentage discrepancies between the simulated and observed spring discharge for each model combination are shown in Fig. 6. Figure 7 shows the interquartile ranges of all model evaluations, including the results of the parametric uncertainty study performed in Teixeira Parente et al. (2019) and the simulated spring discharge obtained from the original Kerschbaum LuKARS model. Comparing the interquartile ranges of the different model combinations (Fig. 7), it is seen that including more uncertain hydrological processes does not necessarily lead to an increase in the output variance. As an example, Fig. 8a,b shows two different cases, characterized by an increase and a decrease in output variance with increasing process complexity, respectively. Figure 8a compares the model combinations Liu, Liu-Magnusson and Liu-Magnusson-Oudin, therefore introducing progressive uncertainties in interception, snowmelt and evapotranspiration. Here, the normalized interquartile ranges increase with the number of hydrological inputs considered as uncertain. Figure 8b considers the model combinations Magnusson, Magnusson-Oudin and Gash-Magnusson-Oudin. In contrast to the previous case, the model combination considering snow and evapotranspiration as uncertain input (Magnusson-Oudin) shows larger output variability than the model considering uncertainties in all the three processes (Gash-Magnusson-Oudin). Looking at all model combinations in Fig. 7, the highest normalized interquartile ranges are noted for model combinations considering snowmelt and evapotranspiration to be uncertain.

Discussion
When considering the uncertainties related to single processes, no significant seasonal variation is observed in the normalized interquartile range related to uncertainties in interception. In the particular case of a broadleaf forest, in Waidhofen a.d. Ybbs beech forest, a more pronounced seasonal variation should be expected due to the higher interception capacity of the leafs in the summer period. This change in canopy cover is, however, not considered in the modeling approaches of Gash et al. (1995) and Liu (2001). In order to obtain more realistic interception estimates, future works should represent variable canopy cover by considering the gross rainfall needed to saturate the canopy, i.e. P′ g , to be changing over time rather than constant. This could be achieved by a temporally varying stand storage capacity (C m ).
In contrast to uncertainties related to interception, a pronounced seasonal variation characterises the uncertainties related to evapotranspiration. The reason why the uncertainties in evapotranspiration are higher in summer 2007 as compared to summer 2006 can be found in the mean summer temperatures of both years (Fig. 2). The mean temperature between April and September in 2006 was 13.17°C. In comparison, a mean temperature of 14.16°C was observed in the same period in 2007. This difference of 1°C leads to the observed increased uncertainties in ET pot . Thus, the results using temperature-based approaches for computing ET pot show that Fig. 6 The bars show the minimum and maximum percentage discrepancies between each model combination and the observed spring discharge the uncertainties of ET pot increase with the available temperature for evapotranspiration. Given the fact that evapotranspiration related uncertainties can even reach the range of parametric uncertainties, e.g. in summer 2007, an appropriate representation of evapotranspiration is crucial to reasonably calculate the groundwater recharge as input for modeling a karst spring discharge in pre-alpine karst systems. Given the seasonal variation of uncertainties related to evapotranspiration, a reasonable representation of evapotranspiration for computing groundwater recharge is even more important during the summer period. This specific knowledge can guide researchers in gathering better field data when specific discharge conditions, e.g. mean and low flow conditions in summer, should be favored in model calibration. As recent studies highlighted the specific role of snowmelt for groundwater recharge in alpine and pre-alpine catchments, this study further investigates if this specific importance is also reflected in increased uncertainties in modeled spring discharge when snowmelt is   Figure 2b shows that in the winter season 2005-2006 the snow cover stayed for several months. Whereas, no long-lasting and pronounced snow cover was observed in the winter season 2006-2007. The same pattern of uncertainties can also be identified when considering snowmelt to be uncertain in combination with other processes, i.e. interception and evapotranspiration. Figure 9 shows the maximum interquartile range for all model combinations including uncertainties in snow processes, i.e. 18 model combinations, and for those which do not include these uncertainties, i.e. 8 model combinations. In case of uncertain snowmelt, the overall interquartile range significantly changes when snowmelt controls groundwater recharge and, thus, the modeled spring discharge (e.g. from January 2006 to May 2006). Moreover, while only considering snowmelt estimations to be uncertain does not lead to a significant increase in model output uncertainties in the evaporative season (Fig. 5), Fig. 9 shows that snowmelt estimations increase model output uncertainties in cases when snowmelt, interception and evapotranspiration are uncertain, e.g. in summer 2007. This can be explained by the fact that different processes can compensate for over-or underestimated water budgets of other processes. For example, an overestimation of the snowmelt can be compensated by an underestimation of the evapotranspiration. Summarizing, the results show that the higher uncertainties in snowmelt occur when the simulated spring discharge is controlled by snow processes. This highlights that the specific importance of snowmelt for groundwater recharge can be identified in the snowmelt related uncertainties when modeling a karst spring discharge. Hence, an erroneous assessment of snowmelt-related groundwater recharge can negatively affect the simulated spring discharge. On the contrary, the results do not exhibit a clear impact of uncertainties in snowmelt on the spring discharge during the evaporative season in summer, when the snowmelt related uncertainties are not significant (<2% of the observed discharge, Fig. 5).
In the specific case of the Kerschbaum LuKARS model, the baseflow storage of the dolomite-dominated aquifer has a high storage capacity and is not immediately affected by changing hydrometeorological conditions. However, this can be different in more limestone-dominated karst systems and requires further investigations.
Finally, the results of all model combinations highlight that considering more processes to be uncertain does not necessarily lead to an increase of the normalized interquartile range of modeled spring discharge. This is particularly true when considering evapotranspiration and snowmelt to be uncertain compared to considering also interception to be uncertain. These results highlight once more that two uncertain processes can compensate each other, leading to a reduction of the interquartile range of model outputs.

Conclusion
This study investigated how the input uncertainties of a lumped parameter model, i.e. LuKARS, vary temporally when applying the hydrologic model in a pre-alpine recharge area. Therefore, snowmelt, evapotranspiration and interception were computed with three different algorithms each. The resulting groundwater recharge was used as input for LuKARS considering each possible model combination and focusing on the uncertainties of each single process. Moreover, the input uncertainties imposed by each process Fig. 9 Contribution of uncertainties related to snowmelt to the total LuKARS model input uncertainties. The two bands show the maximum interquartile ranges of the investigated model combinations that include uncertainties in snowmelt (grey) and of those which do not (red). An effect of snow process uncertainties can be observed throughout the years with a more pronounced impact during the winter were compared to the parametric uncertainties obtained in previous studies (Teixeira Parente et al. 2019).
No clear tendency towards increasing model output uncertainties can be identified when more hydrological input time series are considered to be uncertain. Indeed, it was found that with increasing number of uncertain input the interquartile range of the LuKARS model outputs can even decrease. This shows that two or more uncertain processes can compensate each other.
The results of this study further show that model input uncertainties show temporal variations depending on how much the groundwater recharge and the modeled spring discharge is controlled by one specific process, e.g. snowmelt and evapotranspiration. Thus, the results highlight that the importance of a specific process for groundwater recharge can be derived from the respective input uncertainties. Further, this research identified that uncertainties in snow processes can even be higher than parametric uncertainties.
This study investigated the time-dependent relevance of the model input uncertainties for pre-alpine conditions typical of Central Europe. A similar approach should be applied to ex-tend these results to karst catchments with different climate conditions and land uses. An intercomparison study could be based on the recently developed WoKaS database (Olarinoye et al. 2020). Moreover, it is of particular interest to apply the presented methodology to karst catchments characterized by different recharge processes. Thus, a comparison between catchments dominated by diffusive recharge, as the Kerschbaum springshed, and those dominated by point infiltration is suggested.
The knowledge gained from investigating temporally varying model input uncertainties can guide researchers and water managers in gaining relevant data needed for improving the reliability of hydrologic model results. In this case, e.g., the uncertainties in snowmelt could be reduced by implementing snow measurement stations in the recharge area. Moreover, the information on the temporal variability of model input uncertainties helps to derive which data are needed to improve the reliability of model output results during different times of year.

Declarations
Conflict of interest On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.