Model-aided quantification of dissolved carbon and nitrogen release after windthrow disturbance in an Austrian karst system

Abstract. Karst systems are important for drinking water supply. Future climate projections indicate increasing temperature and a higher frequency of strong weather events. Both will influence the availability and quality of water provided from karst regions. Forest disturbances such as windthrow can disrupt ecosystem cycles and cause pronounced nutrient losses from the ecosystems. In this study, we consider the time period before and after the wind disturbance period (2007/08) to identify impacts on DIN (dissolved inorganic nitrogen) and DOC (dissolved organic carbon) with a process-based flow and solute transport simulation model. When calibrated and validated before the disturbance, the model disregards the forest disturbance and its consequences on DIN and DOC production and leaching. It can therefore be used as a baseline for the undisturbed system and as a tool for the quantification of additional nutrient production. Our results indicate that the forest disturbance by windthrow results in a significant increase of DIN production lasting  ∼  3.7 years and exceeding the pre-disturbance average by 2.7 kg ha−1 a−1 corresponding to an increase of 53 %. There were no significant changes in DOC concentrations. With simulated transit time distributions we show that the impact on DIN travels through the hydrological system within some months. However, a small fraction of the system outflow (  1 year.


Introduction
Karst systems contribute around 50 % to Austria's drinking water supply (COST, 1995).Karst develops due to the dissolvability of carbonate rock (Ford and Williams, 2007) and it results in strong heterogeneity of subsurface flow and storage characteristics (Bakalowicz, 2005).The resulting complex hydrological behaviour requires adapted field investigation techniques (Goldscheider and Drew, 2007).Future climate trajectories indicate increasing temperature (Christensen et al., 2007) and a higher frequency of hydrological extremes (Dai, 2012;Hirabayashi et al., 2013).Both will influence the availability and quality of water provided from karst regions because temperature triggers numerous biogeochemical processes and fast throughflow water has a disproportional effect upon water quality.Furthermore, forest disturbances (windthrows, insect infestations, droughts) pose a threat to water quality through the mobilization of potential pollutants, and these disturbances are likely to increase in the future (Johnson et al., 2010;Seidl et al., 2014).
One way to quantify the impact of changes in climatic boundary conditions on the hydrological cycle is through use of simulation models.Special model structures have to be applied for karst regions to account for their particular hydrological behaviour (Hartmann et al., 2014a).A range of models of varying complexity are available from the literature that deal with the karstic heterogeneity, such as groundwater flow in the rock fracture matrix and dissolution conduits (Jourde et al., 2015;Kordilla et al., 2012), varying recharge areas (Hartmann et al., 2013a;Le Moine et al., 2008) (Rimmer and Salingar, 2006;Tritz et al., 2011).
Nitrate and dissolved organic carbon (DOC) have both been considered in drinking water directives and water preparation processes (Gough et al., 2014;Mikkelson et al., 2013;Tissier et al., 2013;Weishaar et al., 2003).Though nitrate pollution of drinking water is usually attributed to fertilization of crops and grassland, an excess input of atmospheric nitrogen (N) from industry, traffic and agriculture into forests has caused reasonable nitrate losses from forest areas (Butterbach-Bahl et al., 2011;Erisman and Vries, 2000;Gundersen et al., 2006;Kiese et al., 2011).The Northern Limestone Alps area is exposed to particularly high nitrogen deposition (Rogora et al., 2006) and nitrate leaching occurs at increased rates (Jost et al., 2010).In addition to this, forest disturbances such as windthrow and insect outbreaks disrupt the N cycle and cause pronounced nitrate losses from the soils, at least in N-saturated systems that received elevated N deposition due to elevated NOx in the atmosphere (Bernal et al., 2012;Griffin et al., 2011;Huber, 2005).In contrast to N deposition, atmospheric deposition of DOC is low (Lindroos et al., 2008) and thus has not been identified as major driver of DOC leaching from subsoil (Fröberg et al., 2007;Kaiser and Kalbitz, 2012;Verstraeten et al., 2014).Moreover, studies show contrasting results but point to increased DOC (TOC) leaching from soil and catchments after forest disturbances (Huber et al., 2004;Löfgren et al., 2014;Meyer et al., 1983;Mikkelson et al., 2013;Wu et al., 2014).
While many studies identify N and DOC as source of contamination in karst systems (Einsiedl et al., 2005;Jost et al., 2010;Katz et al., 2001Katz et al., , 2004;;Tissier et al., 2013) or provide static vulnerability maps (Andreo et al., 2008;Doerfliger et al., 1999), only very few studies use models to quantify the temporal behaviour of a contamination through the systems (Butscher and Huggenberger, 2008).Some studies use N and DOC to better understand karst processes (Charlier et al., 2012;Mahler and Garner, 2009;Pinault et al., 2001) or for advanced karst model calibration (Hartmann et al., 2013b(Hartmann et al., , 2014b)), but to our knowledge there are no applications of such approaches to quantify the drainage processes of N and DOC, and particularly so after strong impacts on ecosystems (e.g.windthrow) that release reasonable amount of nitrate from the forest soils.
In this study, we consider the time period before and after storm Kyrill (early 2007) and several other storm events (2008) that hit central Europe.The storms, from now on referred to as the wind disturbance period, caused strong damage to the forests in our study area, a dolomite karst system.We apply a new type of semi-distributed model that considers the spatial heterogeneity of the karst system by distribution functions.We aimed at comparing the hydrological and hydrochemical behaviour (DOC, DIN) of the system before and during the wind disturbing period.In particular, we wanted to understand if and how DOC and DIN input to the hydrological system changed by the impact of the storms.Furthermore, we used virtual tracer experiments to create transit time distributions that expressed how the impact of the storms propagated through the variable dynamic flow paths of the karst system.This allowed us to assess the vulnerability of the karst catchment to such impacts.

Study site
The study site, LTER Zöbelboden, is located in the northern part of the national park "Kalkalpen" (Fig. 1).Its altitude ranges from 550 to 956 m a.s.l. and its area is ∼ 5.7 km 2 .Mean monthly temperature varies from −1 • C in January to 15.5 • C in August.The average temperature is 7.2 • C (at 900 m a.s.l.).Annual precipitation ranges from 1500 to 1800 mm and snow accumulates commonly between October and May with an average duration of about 4 months.The mean N deposition in bulk precipitation between 1993 and 2006 was 18.7 kg N ha −1 yr −1 , out of which 15.3 kg N (82 %) was inorganic (approximately half as NO − 3 -N and half as NH + 4 -N; Jost et al., 2010).Due to the dominance of dolomite, the catchment is not as heavily karstified as limestone karst systems, but shows typical karst features such as conduits and sink holes (Jost et al., 2010).The site can be split into steep slopes (30-70 • , 550-850 m a.s.l.) and a plateau (850-950 m a.s.l.), with the plateau covering ∼ 0.6 km 2 .Chromic Cambisols and Hydromorphic Stagnosols with an average thickness of 50 cm and Lithic and Rendzic Leptosols with an average thickness of 12 cm can be found at the plateau and the slopes, respectively (WRB, 2006).Both the plateau and the slopes are mainly covered by forest.Norway spruce (Picea abies (L.) Karst.)interspersed with beech (Fagus sylvatica L.) was planted after a clear cut around the year 1910.The vegetation at the slopes is dominated by semi-natural mixed mountain forest with beech (Fagus sylvatica) as the dominant species, Norway spruce (P.abies), maple (Acer pseudoplatanus), and ash (Fraxinus excelsior).At the slopes no forest management has been conducted since the establishment of the national park.

Available data
A 10-year record of input and output observations was available.Starting from the hydrological year 2002/03, it envelops well the stormy period that began in January 2007.It included daily rainfall measurements and stream discharge measurements from stream sections 1 and 2 (Fig. 1).We obtained the discharge of the entire system with a simple topography-based up-scaling procedure that is described in more detail in Hartmann et al. (2012a).Irregular (weekly to monthly) observations of DOC, DIN and SO 2− 4 concentrations are available for precipitation and at weir 1. DOC, NO − 3 , SO 2− 4 and NH + 4 (since January 2010) samples were filtered (0.4 µm) before the analysis.NH + 4 concentrations were measured by spectrophotometry (Milton Roy Spec- 4 concentrations in runoff are very small or not detectable.Therefore we calculated DIN outputs as NO − 3 -N.Additionally, irregular observations of snow water equivalent at the plateau allowed for independent setup of the snow routines.

Recent disturbances
Kyrill in the year 2007 and some similarly strong storms subsequent to 2008 caused some major windthrows as well as single tree damages.A windthrow disturbance of ∼ 5 ha occurred upstream of weir 1.Though no direct measurements exist as to the total extent of the windthrow area we estimate that 5-10 % of the study site has been subject to windthrow (Kobler et al., 2015).We did not observe a significant change in intra-and interannual variability in DOC concentrations and discharge before and during the wind disturbance period (Fig. 2a, e).Runoff concentrations of DIN showed clear responses to the disturbances.With the first windthrow event it started to increase until 2008/09 and slowly decreased again in 2010/11 (Fig. 2c).Comparison of DOC concentrations with discharge before and during the wind disturbance period revealed a similar pattern.As shown by other studies on DOC mobilization (e.g., Raymond and Saiers, 2010), a positive correlation between concentrations and discharge (on log10 scale) occurred for DOC with concentrations up to 6 mg L −1 during high discharge (similar to Hagedorn et al., 2000).However, there was no obvious difference between the pre-disturbance period (Fig. 2b).

Model hydrodynamics
The semi-distributed simulation model considers the variability in karst system properties by statistical distribution functions spread over Z = 15 model compartments (Fig. 3).That way it simulates a range of variably dynamic pathways through the karst system.The detailed equations of the model hydrodynamics are similar to its previous applications (Hartmann et al., 2013a(Hartmann et al., , c, 2014b)).They are described in the Appendix.Since in our case the model is used to simulate the discharge of the entire system and a weir within the system, some small modifications had to be performed.Preceding studies showed that weir 1 (Fig. 1) receives its discharge partially from the epikarst and partially from the groundwater, reaching it partially as concentrated and partially as diffuse flow (Hartmann et al., 2012a).Consequently we derive its where f Epi is the fraction from the epikarst and (1 − f Epi ) the fraction from the groundwater.f Epi,conc and f GW,conc represent the concentrated flow fractions of the epikarst and groundwater contributions, respectively.Table 1 lists all model parameters, with a short description for each.

Model solute transport
To model the non-conservative transport of DOC and, DIN and SO 2− 4 , we equipped the model with solute transport routines.SO 2− 4 was included as an additional calibration variable because it proved to be important to reduce model equifinality (Beven, 2006) by adding additional information about groundwater dynamics (Hartmann et al., 2013a, b).The inclusion of these three solutes allowed for a more reliable estimation of model parameters (Hartmann et al., 2012b(Hartmann et al., , 2013a) ) and, later on, the evaluation of possible changes in the dynamic of solute concentrations during the stormy period.For most of the model compartments solute transport simulations simply followed the assumption of complete mixing.However, to represent net production and leaching of DOC and DIN in the soil, as well as dissolution of SO 2− 4 in the rock matrix, additional processes were included in the model structure.Similar to preceding studies (Hartmann et al., 2013a(Hartmann et al., , 2014b where a Geo [-] is another variability parameter and G max,SO 4 [mg L −1 ] is the equilibrium concentration of SO 2− 4 in the matrix.DOC is mostly mobilized at the forest floor (Borken et al., 2011).Stored in the soil or diffusively and slowly passing downwards, large parts of the DOC are absorbed or consumed by microorganisms, but when lateral flow and concentrated infiltration increase, net leaching of DOC increases as well.For that reason our DOC transport routine only provides water to the epikarst when it is saturated (Eq.A4) with increasing DOC net production toward the more dynamic model compartments (Fig. 3).Its DOC concentration P DOC,i [mg L −1 ] for each model compartment is found by where a DOC [-] is the DOC variability constant and P DOC [mg L −1 ] is the DOC net production at soil compartment 1.Similar to other studies that assessed N input to a karst system (Pinault et al., 2001), we used a trigonometric series to assess the time variant net production of DIN, P DIN,i [mg L −1 ], to the soil: Here, P DIN is the mean amount of dissolved inorganic N in the soil solution, while A N [mg L −1 ] and S PH,DIN [d] are the amplitude of the seasonal signal and the phase shift of seasonal DIN uptake (immobilization by plants and soil organisms) and release (net DIN in the soil water) cycle, respectively.J D is the Julian day of each calendar year.Due to its seasonal variation, P DIN,i can also be negative, meaning that uptake of DIN takes place.

Model calibration and evaluation
With 14 model parameter that controlled the hydrodynamics and 7 parameters that allow for the non-conservative solute transport, the calibration of the model was a highdimensional problem.For that reason we have chosen the Shuffled Complex Evolution Metropolis (SCEM) algorithm (Vrugt et al., 2003), which proved to be capable of exploring high-dimensional optimization problems (Fenicia et al., 2014;Feyen et al., 2007;Vrugt et al., 2006).As performance measure we used the Kling-Gupta efficiency (KGE; Gupta et al., 2009).For calibration, KGE was weighted equally among all solutes, one-third for the discharge of the entire system and two-thirds for the discharge of weir 1, whose observation precision was regarded to be more reliable than the up-scaled discharge.KGE is defined as where r is the linear correlation coefficient between simulations and observations, and µ S /µ O and σ S /σ O are the means and standard deviations of simulations and observations, respectively.α expresses the variability and β the bias.
To check for the stability of the calibrated parameters, we perform a split-sample test (Klemeš, 1986).Since the predisturbance time series was too short to be split into two equally long periods, we perform a both-sided split-sample test by bootstrapping two independent 4-year time series of observations (first sample: discrete sampling of 50 % of the values of each observed time series; second sample: remaining 50 % of the observations).We calibrate our model with the first sample and evaluate it with the second sample, and vice versa.A parameter set is regarded stable when the calibration with both samples yields similar parameter sets and their KGE concerning discharge and the solutes does not reduce significantly when applying them to the other sample.

Change in hydrochemical behaviour with the stormy period
After the model evaluation, we use the different components of the KGE in Eqs. ( 5) and ( 6) to explore the impacts of the storm disturbance period on the hydrochemical components.By assuming that the model is able to predict to hydrochemical behaviour that prevailed without the impact of the storms adapting the hydrochemical parameters of the model in Eqs. ( 3)-( 4) and analysing the difference between the adapted hydrochemical simulations and the non-adapted simulations, we are able to quantify the change in solute mass balance due to the storm impact.We define the time span for our adaption as the time when the different components of KGE exceed the range of their pre-disturbance variability.
During this time period we compensate for the apparent dewww.biogeosciences.net/13/159/2016/Biogeosciences, 13, 159-174, 2016 viations by adapting the hydrochemical parameters.This is done twice -once by manual adaption and another time using an automatic calibration scheme.Their new values will indicate changes in the seasonality, production or interannual variations.

Transit time distributions
The signal of the storm impact will travel at various velocities and via pathways through the karst system.While fast flow paths and small storages will transport the signal rapidly to the system outlet, slow pathways and large storages will delay and dilute the signal.Transit time distributions indicate how fast surface impacts travel through the hydrological system.We derive transit time distributions from the model by performing a virtual tracer experiment with continuous injection over the entire catchment at the beginning of the impact of the stormy period.When a model compartment reaches 50 % of the tracer concentration is considered median transit time.The thus-derived transit times will elaborate how the hydrological system propagates the signal through the system including all slow and fast pathways as defined by Eqs. ( 12) and ( 18).As for DIN and DOC we assume complete and instantaneous mixing with each model storage (soil, epikarst, and groundwater) at each compartment; the time that we refer to as "mean transit time" of a model compartment is the time the virtual tracer needs to pass through the particular model storage.In combination with the fluxes that are provided from each of the model compartments, it is possible to quantify the fractional contribution of fast and slow flow paths, respectively.We will apply the virtual tracer from the previously assessed beginning of the impact until the end of the time series to assess the transit time distribution.In addition, we apply a second virtual tracer that also lasts only for the disturbance period (as estimated in Sect.3.3) to evaluate the filter and retardation potential of the karst system.

Model performance
Table 1 shows the calibrated parameters for the two samples.They indicate a thick soil and a relatively thin epikarst.The dynamics expressed by the storage constants indicate days and weeks for the conduits (model compartment i = Z) and the epikarst, respectively.The distribution coefficient of the groundwater is larger than the soil/epikarst storage constant.For DOC and DIN there are a natural production rates of 1.6-1.8 and −1.35−0.1 mg L −1 , respectively.The DOC distribution coefficient is between 0.9 and 1.1.The phase shift and amplitude for DIN showed that there is a seasonal variation in DIN net production, with its maximum release at April each year for both of the samples.SO 2− 4 is dominated by the concentration in the precipitation input with some leaching in the soil and sulfides in the dolomite.Its variability constant is quite low (< 0.1).Weighted KGEs, as well as their values for the individual simulation variables, are relatively stable.Overall, calibration on both samples provided similar parameter values.Due to its higher stability concerning the evaluation period, we chose the second sample for further analysis.
The discharge simulations follow adequately the variations in the observations (Fig. 4), although some small events are not reproduced by the model and although the simulations of the weir's discharge tend to underestimate peak flows.No obvious differences can be seen between the pre-disturbance and wind disturbance period.The hydrochemical simulations tend to follow the observations as well (Fig. 5).However, there is sometimes some underestimation of the DOC peaks for the pre-disturbance period.The DIN simulations appear to be more precise during the pre-disturbance period, but there is a systematic underestimation when the disturbance takes place.

Model performance during the wind disturbance period
There is a deviation between pre-disturbance and disturbance period simulated and observed variability and bias for DIN (Fig. 6).A similar tendency can be found for DOC.However, only for DIN are the deviations different to the variations already found during the pre-disturbance period (which is also the calibration/validation period).The variations in DOC appear to be systematic, too, but they fall within its ranges of variability during the pre-disturbance period.

Adaption of N parameters for the wind disturbance period
The very first signs of the impact were found on 1 May 2007, and they lasted to the end of the hydrological year 2010/11.In a first trial (Table 2), the model parameters for the DIN production were adapted manually to compensate for the changes in observed DIN concentrations with a focus on reducing the difference indicated by the bias, β, and variability, α, components of KGE DIN .In a second trial, we use an automatic calibration scheme to achieve the optimum KGE DIN .As indicated by the highest KGE (Table 2), the automatic calibration provided the highest KGE DIN , but this is achieved by improving variability α and correlation r.Almost no improvement is reached for the bias β.Even though resulting in a slightly lower improvement of KGE DIN the manual calibration results in a much more acceptable reduction in the bias (Fig. 6).Its parameter values showed a production rate, P DIN , of DIN almost 2 mg L −1 larger than the pre-disturbance value; an amplitude, A DIN , around 1 mg L −1 smaller; and a phase shift, S PH,DIN , towards a week earlier in the year, resulting in a more acceptable simulation of DIN dynamics during the disturbance period (Fig. 7).

Transit time distributions
The transit time distributions show that the soil and epikarst system reacts quite rapidly to the virtual injection.Fifty percent of the injection concentration is reached within ∼ 60 days (Fig. 8a), while most of groundwater system requires ∼ 100 days to reach 50 % of the injection concentration, with a few flow paths requiring up to 300 days (Fig. 8c).A similar behaviour is found when the impact ends (Fig. 8b, d).This behaviour also shows that some of the slowest flow paths just reach the input concentration before they start to decline again.

Reliability of calibrated parameters and model simulations
Most of the calibrated model parameters are in ranges that are in accordance with other modelling studies or field evidence.General differences between the calibrated parameter values of the both-sided split-sample test may mostly be due to the comparatively low resolution of the hydrochemical variables (SO 4 , DOC and DIN), which actually increased as a result of the bootstrapping procedure.However, the good multi-objective simulation performance of the model, as well as its evaluation by the split-sample test, indicate an overall acceptable performance of the model.At almost 3-8 days the epikarst storage constant is in accordance with field studies on the epikarst storage behaviour that found retention times of some days to a few weeks (Aquilina et al., 2006; www.biogeosciences.net/13/159/2016/Biogeosciences, 13, 159-174, 2016 for comparison the KGE components and their interannual variability are also shown for the pre-storm period and after the correction of the DIN production model parameters during the wind period.Perrin et al., 2003).The soil and the epikarst storage capacity are quite large.These high values may be explained by structural errors of the model that result in unrealistic calibrated parameter values, in particular possible parameter interactions between their storage capacities and storage coefficients.Since the soil and the vegetation controls the fraction of rain that is lost to evapotranspiration, this high calibrated value might be due to tree roots ranging through the soil into the epikarst (Heilman et al., 2012) or rock debris (Hartmann et al., 2012a).Similar to the epikarst storage constant, the conduit storage constant, K C , is, with its value of 1.1 days, in the range of previous modelling studies (Fleury et al., 2007;Hartmann et al., 2013a).The high values of the epikarst variability constant and the groundwater constant indicate a low development of preferential flow paths in the rock, which is typical for dolomite aquifers (Ford and Williams, 2007).A low degree of karstification was already known for our study site (Jost et al., 2010) and the calibrated recharge areas fall well within the ranges found in previous modelling studies (Hartmann et al., 2012a(Hartmann et al., , 2013c)).
The hydrochemical parameters mostly show realistic values.A DOC production parameter, P DOC , of ∼ 1.6-1.8mg L −1 resulted in realistic simulated concentrations at the weir.For DIN production the two calibration samples result in values of −1.4 and 0.1 mg L −1 , going along with amplitudes of 3.4 and 1.8, respectively.Hence, there appears to be some correlation between the production and amplitude parameters, P DIN and A DIN .Negative values indicate that, during some periods of the year, all DIN is consumed by plants or soil organisms and that the production period is shorter but more pronounced due to its larger value of amplitude.However, we expect these differences to be minor since the phase shift S PH,DIN of both calibration samples is almost the same, as well as their annual maximum (P DIN + A DIN ) of 2.01 and 1.95 mg L −1 .This behaviour indicates a maximum of DIN production and leaching at the time of the year when snowmelt reaches its maximum (March to April) and when DIN uptake by plants is still low (Jost et al., 2010).The dissolution equilibrium concentrations of 2.7-3.1 mg L −1 for SO 2− 4 indicate the abundance of the precipitation input, oxidation of sulfides (e.g.pyrite) in the dolomite and traces of evaporates in the small Plattenkalk occurrences (Kralik et al., 2006).

Impact of storms
The deviation between simulated and observed time series (Fig. 5) already indicates that DIN is the only solute that shows a clear impact of the storms.This is further corroborated by considering the individual components of KGE in Fig. 6.It is well known that nitrate leaching to the groundwater increases sharply after tree damage (dieback) in forests where N is not strongly limited (Bernal et al., 2012;Griffin et al., 2011;Huber, 2005).Such disturbances disrupt the N cycle.The loss of tree N uptake favours nitrification of surplus NH + 4 by microorganisms.Moreover, above-(i.e.foliage) and belowground (i.e.fine roots) litter from dead trees enhances the mineralization of organic matter, ammonification and nitrification.Both processes are accelerated by increased soil moisture and soil temperature due to the loss of the forest canopy.Subsequently, leaching of N increases with increased seepage fluxes due to decreased interception and water uptake by trees.Since the simple DIN routine of the model cannot take into account such changes, the underestimated DIN concentrations and their amplitude show the effect of forest disturbance on the leaching of DIN from the studied catchment.There is also an apparently systematic deviation of the DOC variability α.But its variations during the pre-storm pe-  riod are similarly large, thus pointing to a negligible effect of forest disturbance on DOC leaching.Numerous studies have identified the forest floor as a DOC source (Borken et al., 2011;Michalzik et al., 2001).Windthrow generally causes a (short-term) pulse of above-and belowground litter (Harmon et al., 2011).Thereby, mineralization of the surplus litter input concurrent with improved soil climatic conditions likely increased the leaching of DOC from the forest floor (Fröberg et al., 2007;Kalbitz et al., 2007).Concurrent, increased soil water, surface and shallow subsurface flow may favour increased soil DOC leaching to downslope surface waters (Monteith et al., 2006;Neff and Asner, 2001;Sanderman et al., 2009).In mountainous catchments the latter flow paths are likely due to the steepness of the catchment slopes (Boyer et al., 1997;Sakamoto et al., 1999;Terajima and Moriizumi, 2013).The missing signal of forest disturbance on DOC concentrations at weir 1 even shortly after the disturbance may be due to the minor extension of the disturbed area, the minor increase of surface and shallow subsurface flow due to the relative low slope of the disturbed area, the buffering of increased topsoil DOC leaching due to absorption of DOC within the subsoil (Borken et al., 2011;Huber et al., 2004), missing DOC-rich riparian source areas (i.e.wetlands, floodplains) and the reduction in pre-disturbance organic matter input to soil (i.e.litter, root exudates; Högberg and Högberg, 2002).Theoretically, hydrological pro-cesses such as a decrease in transpiration or an increase of groundwater recharge may also occur.But these superficial changes are probably minor considering the typically high karstic infiltration capacities that remove surface water quite rapidly (Hartmann et al., 2014b(Hartmann et al., , 2015)).Therefore, hydrological impacts of windthrow on karst systems (for instance on transpiration) may not be as pronounced as in non-karstic domains because a large fraction of the infiltration during high flow periods will not be available for transpiration anyway.Consequently, a disturbance-caused impact on DOC availability could also be hidden because increased infiltration and DOC leaching during strong rainfall events may just not be detectable considering the weekly to monthly sampling of DOC.For a better understanding of disturbance-induced changes in DOC, more sampling in high temporal resolution of DOC concentrations at the weir (Fig. 1) should be undertaken to elucidate the effect of forest disturbance on DOC dynamics and to improve the simulation of DOC production and transport within the studied ecosystem

N leaching from the soil
Adapting the DIN solute transport parameters through use of an automatic calibration scheme resulted in an increased KGE DIN (Fig. 7).However, it did not resolve the bias of simulated and observed DIN concentrations during the wind disturbance period since the overall improvement of KGE DIN was reached by an improvement of r and α (Table 2).Adjusting the DIN parameters manually resulted in a more acceptable decrease in the bias β that also went along with an increase of the overall KGE DIN .An increase of the DIN production rate of ∼ 2 mg L −1 indicates a massive mobilization of DIN and a reduction in its seasonal amplitude by ∼ 1.1 mg L −1 .Even though there may be some correlation between mean annual production and amplitude (see previous section), the annual maximum of 2.80 mg L −1 (P DIN + A DIN ) indicates an increase of the DIN concentrations in the soil of at least ∼ 0.8 mg L −1 (from 1.95 to 2.01 mg L −1 at the pre-disturbance period).
We identified the beginning of the impact at 1 May 2007 and its end by the end of the hydrological year 2010/11.This is more than 2 years after the last storm in 2008, which indicates how long the ecosystem takes to recover from the disturbance.Other studies have shown comparable recovery times (Katzensteiner, 2003;Weis et al., 2006) or longer (Huber, 2005).Considering the deviations between DIN simulations by the pre-disturbance calibration and the DIN simulations obtained by the manual adjustment, they sum up to an additional release of 9.9 kg ha −1 of DIN over the whole period of ∼ 3.7 years, or 2.7 kg ha −1 a −1 in addition to 5.8 kg ha −1 a −1 that would have been released without the wind disturbance.These values only correspond to inorganic N. Other studies showed that dissolved organic N can also contribute to vertical percolation but only in small ratios from 2 to 5 % (Solinger et al., 2001;Wu et al., 2009).
The apparent shift of S PH,DIN towards an earlier maximum of DIN release (7 days) is most probably be due to the earlier onset of snowmelt in open areas as compared to forests because snowmelt is a major driver of DIN leaching from the soils in our study area (Jost et al., 2010).However, due to the rather slow melting rates, most of the melting water will slowly/diffusively enter the groundwater system rather than flowing rapidly through the karst conduits.Therefore, a slightly earlier beginning of snowmelt may not be visible at the system outlet due to the slow reaction of the groundwater storage.

N propagation through the hydrological system
The virtual tracer injections that we applied with the beginning of the disturbance period elaborate the hydrological system's filter and retardation capacity.Due to their higher dynamics, the soil and the epikarst system adapt more rapidly to the change within weeks and months.Similar behaviour was also found in previous studies (Hartmann et al., 2012a;Kralik et al., 2009).The majority of the simulated flow paths adapt to the virtual tracer signal within a few months, which is in accordance with water isotope studies at the weir (Humer and Kralik, 2008;Kralik et al., 2009).However, using age dating (CFC and SF6) and artificial tracer experiments at individual springs within the study area, Kralik et al. (2009) also found ages from several days to several decades.Hence, the majority of transit times found by the virtual tracer experiment reflect the average behaviour of the sub-catchment drained by the weir, which can be regarded as more dominant than observations at individual the springs that rather represent fast and slow flow paths of minor importance.The retardation is also visible from the dynamics of the DIN concentrations just after the end of the disturbance period (beginning of 2011/12, Fig. 7).Even though DIN production is set to predisturbance conditions, it takes almost 4 months for the DIN simulations (by manual calibration) to adapt to their undisturbed concentrations (pre-disturbance calibration).Due to their small contribution (< 5 %), the slower flow paths do not have a significant impact on the retardation capacity of the hydrological system.

Implications
Our results corroborate findings from many other studies that extreme events such as during the wind disturbance period in our study can result in a significant increase in DIN in the runoff, despite the area impacted being relatively small (5-10 % of the watershed).Particularly in karst catchments, such changes can happen quickly and prevail for a significant duration, in our case more than 2 years after the last storm.Due to subsurface heterogeneity, the impact did not travel uniformly through the system.Instead, it split into different pathways and mixed with old water that percolated prior to the impact.In our system, large parts of the water travelled rapidly through the system.However, a smaller number of pathways had large storages of old water and slow flow velocities, resulting in significant retardation.Taking into account that forest disturbances will most probably increase with climate change (Seidl et al., 2014), DIN mobilization as observed in our study may occur more often and become more intense.The hydrological system may dilute and delay rapid shifts of N concentration, and it will "memorize" the impacts for some time.However, our present analysis showed that the timescale of the wind disturbance on DIN production and leaching from the soil exceeds the timescale of transit of the disturbance through the system.This is most probably due to the small size and the subsurface karstic behaviour of our study site favouring faster flow paths and lower system storage than hydrological systems with larger extent or with other types of geology.

Conclusions
In our study we used a process-based semi-distributed karst model to simulate DOC, DIN and SO 2− 4 transport through a dolomite karst system in Austria.We calibrated and validated our model during a 4-year time period just before a series of heavy storms caused strong wind disturbance to the study site's ecosystem.To quantify its impact we ran the model for the entire disturbance period using the parameters we found at the pre-storm period.The deviations between the simulations and the observations gave us an indication that there was a significant shift in DIN mobilization, its seasonal amplitude and its timing.In estimating the beginning and end of the disturbance period we applied a continuous virtual tracer injection to obtain the mean transit times of the karst system.The transit time distributions showed us how the hydrological system filtered and retarded the impact of the disturbance at the system outlet.
Even though our study is only considering one site and one wind disturbance period, it provides some generally applicable conclusions: (1) hydroclimatic extremes such as storms not only create droughts or floods but can also affect water quality; (2) a hydrological system can filter and delay surface impacts, but it may also memorize past impacts, although only at a limited timescale; and (3) water quality models that have been calibrated without consideration of such external impacts will provide poor predictions.For these reasons we believe that future large-scale simulations of water resources need to include water quality simulations that take into account the impact of ecosystem disturbances.Even without anthropogenic contamination, climate change will strongly affect water quality in our aquifers and streams, and we need to understand processes that affect water quality and prepare ourselves in order to avoid threats to future water supply.
or pref-Published by Copernicus Publications on behalf of the European Geosciences Union. A. Hartmann et al.: Model-aided quantification of dissolved carbon and nitrogen release erential recharge by cracks in the soil or fractured rock outcrops

Figure 2 .
Figure 2. Intra-annual and interannual variations in (a) DOC concentrations, (c) DIN concentrations and (e) discharge, and relation between discharge and (b) DOC and (d) DIN before and during the wind disturbance period.

Figure 3 .
Figure 3. Diagram of model structure; it is assumed that discharge and hydrochemistry at the two weirs is composed by different mixtures of diffuse recharge (green), concentrated recharge (red), diffuse groundwater flow (blue) and concentrated groundwater flow (purple).

Figure 4 .
Figure 4. Observed vs. simulated discharges for the entire karst system and weir 1.

Figure 6 .
Figure 6.Individual components of the KGE: (a) ratio of simulated and observed variabilities, (b) ratio of simulated and observed average values, and (c) their correlation for the wind disturbance period; for comparison the KGE components and their interannual variability are also shown for the pre-storm period and after the correction of the DIN production model parameters during the wind period.

Figure 7 .
Figure 7. Observed and simulated DIN dynamics using the pre-storm parameters (red line), the scenario 1 parameters derived from the deviations assessed by the KGE components (orange line), and the scenario 2 parameters derived by systematic variation (dark-red line).

Figure 8 .
Figure 8. Mean transit times for (a) the soil and epikarst and (c) the groundwater storage derived by an infinite virtual tracer injection starting with the beginning of the wind disturbance period, and (b) the reaction of the soil and epikarst and (d) the groundwater storage as the impact ends.

Table 1 .
Model parameters, description, ranges and calibrated values with KGE performances for the calibration and validation samples.

Table 2 .
Calibrated pre-storm parameters for DIN dynamics and two scenarios for adapting the model at the stormy period.