Karst spring recession and classification: efficient, automated 1 methods for both fast and slow flow components 2

. 10 Analysis of karst spring recession hydrographs is essential for determining hydraulic parameters, geometric characteristics and 11 transfer mechanisms that describe the dynamic nature of karst aquifer systems. The extraction and separation of different fast 12 and slow flow components constituting karst spring recession hydrograph typically involve manual and subjective procedures. 13 This subjectivity introduces bias, while manual procedures can introduce errors to the derived parameters representing the 14 system. To provide an alternative recession extraction procedure that is automated, fully objective and easy to apply, we 15 modified traditional streamflow extraction methods to identify components relevant for karst spring recession analysis. 16 Mangin’s karst-specific recession analysis model was fitted to individual extracted recession segments to determine matrix 17 and conduit recession parameters. We introduced different parameters optimisation approaches of the Mangin’s model to 18 increase degree of freedom thereby allowing for more parameters interaction. The modified recession extraction and 19 parameters optimisation approaches were tested on 3 karst springs in different climate conditions. The results show that the 20 modified extraction methods are capable of distinguishing different recession components and derived parameters that 21 reasonably represent the analysed karst systems. We record an average KGE >0.85 among all recession events simulated by 22 the recession parameters derived from all combinations of recession extraction methods and parameters optimisation 23 approaches. While there are variabilities among parameters estimated by different combinations of extraction methods, 24 optimisation approaches and seasons, we find even much higher variability among individual recession events. We provide 25 suggestions to reduce the uncertainty among individual recession events and raise questions on how to improve confidence in 26 system’s attributes derived from recession parameters.


Introduction 28
Groundwater from karst aquifers are essential water sources globally (Stevanović 2018;Goldscheider et al. 2020). Karst 29 aquifers are characterized by high degree of heterogeneity and complex flow dynamics resulting from the interplay of conduit 30 and matrix drainage systems (Kiraly 2003;Goldscheider and Drew 2007). Groundwater flow is rapid in the highly-conductive 31 5 storm runoff influence at the beginning of recession, the first 30% of each recession segment is deleted. Additionally, the 126 difference between two consecutive streamflow values must be ≤ 30% for the pairs to be accepted. All recession segments that 127 satisfy these conditions are then accepted as slow flow recessions segment. 128 129 In order to objectively determine streamflow recession that is derived purely from a dry or low flow period, Brutsaert (2008) 130 introduced a more strict recession extraction method. For streamflow Q, during time t, the Brutsaert method eliminates data 131 point with increased or zero values of dQ/dt, as well as points with abrupt dQ/dt values. To exclude data points that might be 132 influenced by storm runoff, three data points after a positive or zero dQ/dt are removed -in major events, an additional fourth 133 data point is removed. While Brutsaert (2008) did not provide a description for a majors event, Stoelzle et al. (2013) applied 134 the Brutsaert method in their study and defined major event as streamflow values exceeding 30% streamflow frequency. 135 Therefore, our study uses this definition of major event from Stoelzlz et al. (2013). Furthermore, the Brutsaert method also 136 excludes last two data points before a positive or zero dQ/dt and spurious data points with larger -dQ/dt values. 137 The three recession extraction approaches (Vogel and Kroll 1992;Brutsaert 2008; Aksoy and Wittenberg 2011) were 156 specifically developed to extract streamflow recessions that are mainly from slow flow contribution. Thus, rules and exclusion 157 criteria specified by each method aim at eliminating the quickflow influences from the extracted recession segments. In karst 158 systems, concentrated rapid flow through the conduit networks constitutes the quickflow, while the contribution from slow-159 velocity drains through the matrix pores constitutes the slowflow. The quick and slow flow represents flows from two different 160 drainage structures and both contribute to karst spring recession (Fiorillo, 2014;Ford & Williams, 2007;Padilla et al., 1994). 161 162 Adapting the streamflow methods for karst spring recession analysis requires both slow and quick flow components to model 163 matrix and conduit spring discharges. So, to adapt the traditional REMs, we (i) extract spring flow recession curve based on 164 the specific method approach, (ii) attribute part of the recession curve that satisfies the specified method's exclusion criteria 165 as slowflow (matrix) component, and (iii) assign the remaining part that is excluded as quickflow (conduit) component.  where Qro is the baseflow contribution at the beginning of recession when t = 0, α is the recession coefficient with a unit of T -183 1 and t is the lapsed time between discharge at any time t, Qt and initial discharge at t = 0, Q0; and Eq. 4 describes the non-184 linear relationship during quickflow recession from the unsaturated zone. 185 where q0 is the difference between Q0 and Qro, parameter η describes the infiltration rate through the unsaturated zone. The 189 parameter is defined as 1/ti for the duration of quickflow recession between t = 0 and ti = 1/η. ε in T -1 unit describes the 190 regulating capacity of the unsaturated zone during infiltration and characterises importance of concavity of quickflow recession 191 (Padilla et al. 1994). The algebraic sum of Eq. 3 and 4 gives Eq. 5, which defines the discharge at time t during the recession 192 period.

Mangin classification framework 204
Following the estimation of recession parameters α, η and ε using Eqs 3 -5 above, Mangin proposed a classification scheme 205 for karst systems based on two additional parameters: (1) aquifer regulation capacity, K, and (2) infiltration delay, i. To 206 determine K, the dynamic volume, Vdyn, which is defined as the volume of water stored in the phreatic zone at the peak 207 discharge time t0 is calculated using Eq. 6. The average volume of water, Vann, discharged through the spring's outlet over one 208 hydrological year is also calculated. The regulation capacity K, is therefore given by the ratio of Vdyn and Vann as expressed 209 with Eq. 7. This parameter represents the extent of the phreatic zone and its ability to regulate groundwater release from 210 8 storage. While porous aquifers can have values of K > 0.5, a typical karst system is expected to have K < 0.5 (Marsaud 1997

Comparison and evaluation of REMs and POAs 274
The three REMs (Vogel, Brutsaert and Aksoy) are combined with the three POAs (M1, M2 and M3) of the recession model 275 to derive slow and quick flow recession parameters of selected karst springs for a total of nine possible methods. The recession 276 11 parameters are derived separately for both summer and winter recession events. The overall performance of the different REM 277 and POA combination is determined by calculating the goodness of fit between observed spring recession discharges and ones 278 simulated with the derived parameters using Kling Gupta Efficiency (KGE) measures (Gupta et al. 2009). We use KGE because 279 it considers the common model error types -the mean error, variability and the dynamics. The mean and interquartile ranges 280 of the derived parameters are compared among different method pairs and seasons. The estimated recession parameters were 281 used to identify the dynamic of the systems according to Mangin's karst system classification described in subsection 2.2.2. 282 The Mangin classification scheme describes the aquifer drainage characteristics, conduit development and speleological 283 network (Mangin 1975; El-Hakim and Bakalowicz 2007). Therefore, this is used to evaluate the representativeness of recession 284 parameters estimated for the selected karst springs aquifer systems. 285

Test springs and data 286
The REMs and POAs were tested using three karst springs; Lehnbachquellen, Saivu and Qachquoch located in Austria, 287 Switzerland and Lebanon respectively (Figure 2). The selection of these springs were based on geographical spread which 288 covers different climate and hydrological settings, availability of discharge hydrograph in high resolution as well as literature 289 reference on hydrological characterisation of aquifer systems drained by the spring. Daily and sub-daily spring discharge time 290 series of the selected springs were obtained from WoKaS database (Olarinoye et al. 2020). Important characteristics of the 291 spring hydrographs as well as the catchments in which they are sited are presented in (Table 3) Lehnbachquellen is sited in snow-dominated, Saivu in humid and Qachquoch is in the Mediterranean catchment. It should be 294 noted that in snow catchment, recession behaviour will be externally influenced by snow storage. However, we have included 295 snow-dominated catchment in this study to access this impact of this external influence. The spring discharge time series 296 measured at a uniform time-step for each spring span between 3 and 13 years. All discharge time series were aggregated to 297 daily temporal resolution, and missing data values which were only found (<0.01%) in Lehnbachquellen spring discharge data 298 were excluded.

Extracted recessions and performance of POAs 306
The adapted recession extraction methods adequately identified karst spring conduit and matrix flow components. The 307 parameters obtained with the different REM-POA pairs also produced a satisfactory simulations of recession events. Only 308 complete recession events >= 7 days period were considered for analysis. Here, complete recession referred to events that 309 featured both conduit and matrix components. For each spring hydrograph, a different number of recession events are identified 310 by the REMs. As shown on Table 4, Vogel method captured the highest number of recession events across all springs, followed 311 by Brutsaert (except for Lehnbacquellen spring) and Aksoy showed the least ability to capture recession periods from the 312 observed spring discharges. However, the average length of the recession events varied among the different REMs in no 313 13 particular order (see Fig. A2 in appendices). Based on the number of recognizable recession events, the REMs were defined 314 as permissive (Vogel), less permissive (Brutsaert) and restrictive (Aksoy). 315 316

Aquifer characterization 397
To evaluate the overall representativeness of estimated recession parameters based on the modified REMs and different POAs 398 for the selected karst spring systems, we determine the drainage properties of the spring's aquifer using the parameters derived 399 from the individual recession event. As described in subsection 2.2.2, retardation between infiltration and output defined by 400 infiltration delay parameter, i, and aquifer regulation power, K, were calculated for individual recession event. With these high K values, the Lehnbachquellen system has the highest capacity to withhold groundwater among the three karst 426 springs used in this study. The wide dispersion of both K and i makes it impossible to confine the system's into a specific class. 427 The Lehnbachquellen system therefore falls within three classes; class II (well-developed system with large downstream flood 428 plains), class III and class V (poorly developed system).

Quality of extracted recessions 438
With the modification of the traditional REMs, we were able to establish a completely objective approach to distinguish 439 between slow and quick flow recession components. Furthermore, optimisation approaches (POAs) with more flexibility 440 showed better improvement over the conventional parametrisation procedure. The REMs tested use different statistical indices 441 to scrutinise genuine baseflow records, hence they have different level of tolerance. The ability of the extraction methods to 442 identify recession periods from hydrograph time series depend on the level of their restrictiveness. Vogel extraction method 443 defined by a 3-day moving average to smoothen the hydrograph and also allows for 30% increase in subsequent flowrates is 444 more permissive than Brutsaert and Aksoy methods that strictly enforce dQ/dt < 0. Hence, more recession events were 445 extracted by Vogel method. Study by Stoelzle et al. (2013) also showed the Vogel procedure to be more permissive, as it was 446 able to extract almost 50% more events than Brutsaert. Although main recession selection condition for Brutsaert and Aksoy 447 method is determined by decreasing dQ/dt, constraining real baseflow recessions to discharge data points with less than 10% 448 (CV ≤ 0.1) deviations makes the Aksoy more restrictive than Brutsaert method. 449 450 As previous mentioned, higher hydraulic head wold promote faster drainage of the aquifer resulting in higher values of α 512 parameter. 513 514 For quickflow recession parameters, seasonal variability is independent of the REM. The three springs showed different 515 seasonal patterns which could be directly linked to their hydroclimatic settings. Seasonal influence on quickflow recession 516 parameters was not clearly seen in the snow-dominated catchment. This could be attributed to the snow melting process 517 discussed above. Since snow melt compensate hydrologic flow during warmer period, there would be constant influx from the 518 surface throughout the year, also soil wetness conditions do not change significantly. This explains the lack of any evident 519 seasonal differences between parameters η and ε estimated for in Lehnbachquellen spring in the snow-dominated catchment. 520 But the Saivu spring in the humid and Qachquoch spring in the mediterranean catchment showed clear seasonal influences. 521 Estimated values of infiltration rates η for Saivu were higher in summer (lower in winter) and lower in summer (higher in 522 winter) for the Qachquoch spring. This pattern is believed to be controlled by the peculiarity of the different geographic and 523 climatic settings. In humid catchment, higher temperature in summer would result in dryer soil conditions which would 524 consequently facilitate faster infiltration. However, for the mediterranean settings, soil conditions are dry due to relatively 525 warmer temperature all year round. This makes precipitation is a limiting factor, and with more precipitation in winter, faster 526 infiltration through the unsaturated zone would occur. 527 528 5.3 How realistic are adapted REM-POA for karst system analysis? 529 Karst system classification proposed by Mangin (1975) is based on two parameters K and i (see subsection 2.2.2). These two 530 parameters were derived from the estimated recession parameters (α, η and ε), thus the variability found in the recession 531 parameters is expected to be propagated to K and i. Although, if the derived mean values of K were considered, some level of 532 coherency was found among all REM-POA combinations and between the seasons. But looking at the estimated standard 533 deviations, a large intra-event and seasonal variation can be found. In a study by Grasso & Jeannin (1994), the authors found 534 regulation power, K, to be more stable for various years and events. These findings did not agree with our analysis, the 535 outcomes of which showed a large variability among K for different events, most significantly in the snow-dominated 536 catchment. Regulation power is analogous to memory effect, and the periodic water release from an external snow storage that 537 is not captured within the saturated zone in real time makes K to fluctuate more in snow-dominated catchment. Considering 538 the standard deviations from the mean, infact the values of K exceeded the maximum value of 1 originally proposed in Mangin 539 karst classification scheme. Mangin (1975) set a maximum value of one for K, with assumptions that real karst systems would 540 not have a storage memory beyond one year. However, karst system in a snow catchment could have K values greater than 541 one due to snow accumulation and melting as found in Lehnbachquellen spring. Also, complex aquifer systems, as in the case 542 of Qachquoch spring could also have higher K values. 543 22 544 Infiltration delay, i, is strongly dependent on recharge type contribution as well as catchment size (Jeannin and Sauter 1998). 545 Recharge is control by climatic input (rainfall) which varies between seasons. However, the derived values of i were hardly 546 separated by season, but more varied among individual recession events. The complex interplay of REM and POA resulted in 547 a compensation phenomenon; whereby infiltration rate, η, was compensated by recession concavity parameter, ε. Since the 548 infiltration delay is defined by these parameters, it is difficult to explicitly infer the specific effects of REM and POA on 549 infiltration delay. 550

551
The northern Alps karst system where the Lehnbachquellen spring is located has been defined as well karstified highly 552 permeable unit interlayered with less permeable Flysch formation (Goldscheider 2005 study place this spring in class 1, therefore coherently agreeing with the authors description. Taking in account the standard 556 deviations, the classification of Qachquoch spring ranged between medium to poorly karstified system. This is similar to a 557 recent study by Dubois et al. (2020) that categorises the system as poorly karstified with a very large regulation capacity. 558 559 Given that existing common karst spring recession extraction method involves visually supervised procedure and subjectively 560 determined duration of conduit infiltration, alternative faster, automated and objective approach is very useful. From our 561 analysis, resulting parameters of extracted recession segments are within reasonable ranges and derived systems classification 562 correspond to those found in literatures. The good performance recorded between simulated and observed flowrates during 563 recession events attests to the potential transferability of traditional extraction methods to karst systems. However, this good 564 performance does not necessarily translate to reliable parameter estimates. It is therefore important to choose REM methods 565 that gives reasonable parameters especially when paired with a less flexible optimisation approach. Furthermore, with prior 566 knowledge of the spring system, parameters ranges can be reasonably constrained during optimisation to achieve more 567 representative optimum parameters. 568 adapting them to also identify quick flow flow recessions. We fit individual extracted recession segments with Mangin's 576 recession model to determine the conduit and matrix drainages recession characteristics. We introduce new parameters 577 optimisation approaches (POAs) different from the conventional procedure to increase degree of freedom of parameter 578 interaction. 579 580 While we found that there were uncertainties in the estimated recession parameters resulting from the methodological choices 581 (REM and POA combinations) and seasonal influences, the uncertainties among individual recession events were much larger. 582 The large variability among individual event actually reflected the dynamic and heterogenous nature of the karst system. The 583 combination of this with REMs, POAs and seasons resulted in a more complex interplay and only amplified the uncertainties. 584 These uncertainties is actually useful to understand the dynamic nature karst system, but it is difficult to cope with and also 585 need to be systematically quantified. To avoid these large uncertainties, master recession analysis approach has being a popular 586 alternative for karst spring hydrograph analysis. But a single recession parameters' values derive from the master recession 587 approach does not reflect the highly dynamic nature of karst system. The uncertainty of karst recession parameters derived 588 from either single or master recession approach is presently not a discussion in karst hydrology. Maybe such discussion needs 589 to start to address the limitations and difficulties encountered in this study. Herein, we pose two major issues that need to be 590 addressed as seen from this study: (1) how can we do recession analysis more objectively with a single REM and separation 591 technique that accounts for all ranges and possible instances of slow and quick flow? and (2) how can we incorporate a more 592 robust parameters estimation and uncertainty quantification approach into individual recession analysis? Answering these 593 questions will help to expand confidence in the system's drainage characteristics that are derived from recession parameters. 594 595 Finally, this study has shown that there are a lot of potentials for extracting and separating karst spring recession components 596 by adapting the traditional REMs and introducing flexible parameter optimization approaches. The adaptation of the REMs in 597 combination with the different parameters estimation flexibility (POAs) provides a suit of automated tools that can be used for 598 karst recession study. This automated and multi approach for parameters optimization is essential to cope with the known 599 biases of single and visually supervised recession analysis methods. Different REMs has their specific advantages and there 600 are still room for improvement. For example, other extraction can could be tested and non-linear reservoir model can also be 601 considered for fitting the matrix model. The R codes for the different REMs and POAs used for the recession analysis can be accessed through our GitHub repository 609 here https://github.com/KarstHub/Karst-recession