On the use of spring baseflow recession for a more accurate parameterization of aquifer transit time distribution functions

Baseflow recession analysis and groundwater dating have up to now developed as two distinct branches of hydrogeology and have been used to solve entirely different problems. We show that by combining two classical models, namely the Boussinesq equation describing spring baseflow recession, and the exponential piston-flow model used in groundwater dating studies, the parameters describing the transit time distribution of an aquifer can be in some cases estimated to a far more accurate degree than with the latter alone. Under the assumption that the aquifer basis is sub-horizontal, the mean transit time of water in the saturated zone can be estimated from spring baseflow recession. This provides an independent estimate of groundwater transit time that can refine those obtained from tritium measurements. The approach is illustrated in a case study predicting atrazine concentration trend in a series of springs draining the fractured-rock aquifer known as the Luxembourg Sandstone. A transport model calibrated on tritium measurements alone predicted different times to trend reversal following the nationwide ban on atrazine in 2005 with different rates of decrease. For some of the springs, the actual time of trend reversal and the rate of change agreed extremely well with the model calibrated using both tritium measurements and the recession of spring discharge during the dry season. The agreement between predicted and observed values was however poorer for the springs displaying the most gentle recessions, possibly indicating a stronger influence of continuous groundwater recharge during the summer months.


Introduction
Spring baseflow recession analysis has a long history in hydrology, starting more than a century ago when Boussinesq and Maillet proposed using quadratic or exponential laws to describe spring recession (Boussinesq, 1904;Maillet, 1905).Following this seminal work, a number of more complex models have been developed combining two or more reservoirs and considering non-linear responses (Horton, 1933;Brutsaert and Nieber, 1977;Brutsaert, 1994;Coutagne, 1948;Mangin, 1970;Padilla et al., 1994).Dewandel et al. (2003) give an excellent review of the subject, and classify all recession studies into two approaches: the first considering drainage in both saturated and unsaturated zones, and the second concentrating on the recession of the saturated zone only.For the latter, only Boussinesq's quadratic solution is both analytically exact and interpretable in terms of hydraulic parameters of the aquifer (hydraulic conductivity and effective porosity).As was shown from the statistical analysis of the recessions of 100 karstic springs (Drogue, 1972) and using numerical techniques (Dewandel et al., 2003), a quadratic law describes much more truthfully spring recession than an exponential law, although the latter is more popular in groundwater hydrology.Furthermore, the quadratic law, albeit derived from a number of simplifying assumptions, proved robust for more realistic aquifers.Acknowledging this, Boussinesq quadratic law was preferred over an exponential law to describe spring flow recession in this paper.
One of the reasons to study baseflow recession is that since its slope is controlled by the hydrodynamic properties and geometry of the aquifer, it is possible to estimate from the discharge recession an averaged (so-called effective) hydraulic conductivity and storage coefficient (Boussinesq, 1904;Brutsaert and Nieber, 1977;Brutsaert, 1994;Szilagyi et al., 1998;Mendoza et al., 2003).Since the equation relates the volume of water in storage and spring discharge, the mean hydraulic transit time in the saturated zone can also be derived from fitting the Boussinesq quadratic solution to an observed spring recession (a fact that to our knowledge has not received much attention until now).
Mean groundwater residence times are usually estimated using lumped-parameter models (Maloszewski and Zuber, 1982) calibrated on the measurement of environmental tracers such as tritium or CFCs.Since tritium infiltrates conservatively with rainwater, the estimated transit time is the sum of transit times in the unsaturated and saturated zones.Due to the shape of the tritium input function, the solution of the inverse estimation can be non-unique, especially for transit time distributions with more than one parameter or when the tritium record in the outlet is short.If the model is used to predict solute transport time, this non-uniqueness is propagated to the results of the solute transport model.
In this paper, we show that the discharge recession can be used to reduce parameter uncertainties in a model predicting atrazine concentration in spring water over time.The parameters of the model are estimated from tritium measurements and baseflow recession, and predictions then compared with the observed atrazine time series in a second testing step.

Atrazine transfer function
The model predicting atrazine concentration in spring water over time is based on the transit time distribution of the aquifer, representing the sum of flow through all flow lines connected to the outlet and having different transit times.Besides being used to estimate groundwater transit times, convolutions have been applied in hydrogeology to estimate groundwater recharge (Besbes and de Marsily, 1986) and study water table fluctuations (Olsthoorn, 2008).The physical meaning of the transfer function used depends on the geological and hydrogeological setting of the study site, and may even be chosen purely empirically.The transit time distribution function used here is the exponential piston-flow model (EPM) proposed by Maloszewski and Zuber (1982), a combination of two simpler models: piston-flow and exponential.The piston-flow component simulates an unsaturated zone of near-constant thickness where all flow lines are approximately vertical and of equal length and transit time (Farlin et al., 2013), and the exponential component simulates the transit time distribution in the saturated zone.Transit time distributions in real aquifers, which are by nature heterogeneous, most probably differ from simple analytical transit time distributions (Etcheverry, 2001).However, since no data is available concerning the spatial heterogeneity of the fracture network and in the absence of empirical evidence of extreme heterogeneities on the study site (Farlin et al., 2013), the simplest model of a pseudo-homogeneous aquifer was adopted.As was shown by Haitjema (1995) and Etcheverry (2001), ground water transit times are exponentially distributed in a homogeneous aquifer with constant saturated thickness.The EPM has two fitting parameters (which can be expressed as the mean transit time in the unsaturated and saturated zones respectively).
Atrazine concentration in spring water C out is predicted from the input leaching concentration C in by (Farlin et al., 2013) with x cropland /x = ratio of cropland to total catchment area, g (t) = atrazine transit time distribution and λ = atrazine half-life in the sandstone in years (both saturated and unsaturated zones).The model assumes that sorption processes are negligible in the aquifer.The transit time distribution g (τ ) can be estimated as follows: first, the parameters of the transit time distribution of water g (τ ) are estimated from tritium measurements.g (τ ) is defined by Maloszewski and Zuber (1982) as follows: with η = ratio of total volume of water in the stored groundwater system (V ) to the volume of water stored in the reservoir with exponentially distributed transit times (V EM ) [-].
where t epm = total mean transit time of the tracer in the system (in a double porous medium), t pf = mean transit time of the tracer in the unsaturated zone and t em = mean transit time of the tracer in the saturated zone The parameters of g (τ ) are then calculated from those of g (τ ) using the land-use distribution (Farlin et al., 2013).The second step is necessary to take into account the fact that atrazine is only applied on agricultural surfaces, whereas tritium infiltrates homogeneously over the entire recharge area.g (τ ) only differs from g (τ ) by the value of the parameter η and t em , which are a function of the proportion of arable land in the catchment.The greater this proportion, the closer g (τ ) will be to g (τ ).Details can be found in Farlin et al. (2013) and are omitted here for space reasons.Since the relationship between g (τ ) and g (τ ) is bijective and the parameters η and t em can be converted uniquely without ambiguity, we will refer in the rest of the text to the parameters of g (τ ) only.Both parameters of the EPM are a priori unknown and are estimated from environmental tracer measurements (tritium in the present study; see Farlin et al., 2013, for details).The goodness of fit is calculated from the misfit between predicted and observed tritium concentrations: The best fit is obtained by minimizing ε.Theoretically, both t epm and η can be estimated from tritium data.However, in practical applications, η is not always a sensitive parameter and in some cases more than one minimum of Eq. ( 4) may be obtained for η using tritium measurements only, even when t epm can be estimated unambiguously.In other words, the total mean transit time t epm can be estimated, but not separated into its components t em and t pf (respectively transit time in the saturated and unsaturated zone).For that reason, we try to estimate η from the recession hydrograph.Boussinesq (1904) derived an exact analytical solution of the diffusion equation for a set of simplifying assumptions and an aquifer with a horizontal basis.Recession follows then a quadratic law:

Recession curve analysis
The volume of water in storage at any time t is Taking the limit lim t→t 0 =0 V (t), we obtain (Drogue, 1972) According to Maloszewski and Zuber (1982), the mean transit time in the saturated zone (the turnover time) is Comparing Eqs. ( 6b) and ( 7), we see that Knowing t em from the recession and t epm from tritium observations, Eq. ( 3) can be used to estimate η.
Mean transit times in the saturated zone calculated from a tracer and from the recession will be approximately equal if and only if the groundwater volumes concerned are approximately the same.This is not the case for double porous systems where diffusive exchanges of tritium mass take place between a mobile and an immobile water domain.In that situation, water flows only through the hydraulically active part of the system while tritium additionally diffuses in and out of the stagnant water zones.A second possible situation is when the basis of the aquifer is convex, causing the trajectory of some flow lines to plunge below the aquifer outlet and creating a dead volume (referred to by Zuber, 1986, as the minimum volume) that does not influence the discharge rate because it is situated deeper than the outlet itself.In that case, the calculated transit times depend upon the volume of water stored above the datum of the outlet (the dynamic volume), whereas tracer transit times will be calculated for the total volume of water stored in the aquifer (i.e. the sum of dynamic and minimum volumes).Hence, two assumptions have to be made when combining tracer information with hydrograph recession.The first is that, in the case of double porous aquifers, the volume of the immobile domain is negligibly small for both water and tracer transport, the second that the minimum volume is negligibly small compared to the dynamic volume.

Predicting atrazine concentrations in spring water
The history of atrazine soil application was unknown except for two important milestones.The first was the introduction of combination products in the early 1990s, which made overdosing atrazine impractical, as this simultaneously increased the dose of other products present in the formulation, damaging the treated crops.The second date was the nationwide ban on atrazine, which was enforced in 2005.Hence, it is possible to reconstruct a schematic history of atrazine leaching to the groundwater consisting of three stages: a first stage of higher application causing higher leaching until the mid-1990s, followed by decreased leaching up to approximately 2005 (taking into account a certain lag until farmers had exhausted their atrazine stock and the atrazine reservoir in the soil had been sufficiently depleted), and a third stage consisting of atrazine-free recharge water.As a first approximation, we assumed the leaching during the combination product period to be much smaller and negligible compared to the leaching of the pure atrazine phase.Thus, the predictive model consisting of Eq. ( 1) was used with an atrazine time series C in consisting of a single step function with a break placed in the mid-1990s.The parameters of g(τ ) were estimated both from tritium measurements alone, yielding two different models (models "tritium 1" and "tritium 2"), and tritium combined with spring baseflow recession (models "recession").Depending on the data available, one or two recession models were used (corresponding to a parameter k estimated from the recession in 2008 and 2009 respectively).
www.hydrol-earth-syst-sci.net/17/1825/2013/ Finally, the break in C in was shifted until the best agreement was reached between model predictions and observations.This additional fitting step was necessary for two reasons: the change from pure atrazine to combination product probably did not take place instantly, and the soil acts as an additional reservoir that reacts to changing application practices with a lag of a few years (Farlin et al., 2013).

Study area
A series of springs draining the Steinsel plateau, a sandstone cuesta situated 10 km north of Luxembourg City, were sampled regularly over 3 yr (2008 to 2011).The plateau is 10 km long from north to south and 2 km from east to west at its widest.Its central part is farmed intensively while the slopes are covered by mostly deciduous forests.Mean annual precipitation is 800 mm, of which 200 mm recharge the aquifer.The aquifer is part of a fractured sandstone formation named the Luxembourg Sandstone, which provides about half of the country's drinking water.The sandstone is densely fractured, unconfined and is drained by numerous free-flowing contact springs.The mean thickness of the formation on the Steinsel plateau is 70 m with a large unsaturated zone of up to 60 m.Few measurements of matrix and fracture porosities exist for the Luxembourg Sandstone, and none on the study site.Matrix porosity is known to vary with the degree of dissolution of the calcareous cement, from 5 % to 35 % according to Colbach (2005).Fracture porosity is estimated by the same author to be approximately 1 %, but Farlin et al. (2013) calculated (from the mean groundwater transit time, the mean annual recharge and the thickness of the saturated and unsaturated zones) fracture porosity values between 5 and 6 % for the Steinsel plateau nearly accounting for the total porosity.We hypothesize that the mobile water domain corresponds to the fracture network, whereas stagnant water zones could be present in the calcareous matrix.
Twenty springs were sampled either weekly or monthly from 2008 to March 2010, and thrice in 2011.pH, electric conductivity, water temperature as well as spring discharge were measured in the field.Samples were taken to the laboratory for standard chemistry and pesticide analysis.Details of the analytical method for the atrazine samples can be found in Farlin et al. (2013).Tritium was measured twice each year in July and September.The water quality of many springs draining the Luxembourg Sandstone, including those sampled, is impacted by intensive farming practices taking place in recharge areas.Equation (1) was developed to predict the evolution of atrazine concentration in spring water following its countrywide ban in 2005.Atrazine concentrations were stable until 2011, when the beginning of a decreasing trend was observed.This prompted returning to the original predictions made based solely on the data for the period 2008-2010 and assessing with hindsight the model's predictive power.The original model presented in Farlin et al. (2013) was sufficient to explain the inertia of the aquifer system and estimate the time to trend reversal.It suffered however from one major flaw: the tritium data were not sufficient to differentiate between two models predicting different mean transit times in the saturated and unsaturated zones (for the same total transit time within the formation).

Spring recession
Most springs displayed a clear recession lasting from February till as late as June.Table 1 summarizes the t em estimated from the recession periods in 2008 and 2009 and Eq. ( 8).The parameter estimation was performed by least-square fitting, with both Q(t 0 ) and k allowed to vary.For some springs, both years yielded comparable mean transit times in the saturated zone, but for others (K17, K21 and K21a) estimates differ by a factor of four between 2008 and 2009, with the shorter transit time closer to those in the other springs.Note that the recession in 2008 nearly systematically yields longer mean transit times than the recession in 2009 (with the exception of K17 and to a much lesser extent K9).Discharge time series and quadratic fit are shown exemplarily for three springs in Fig. 1, each illustrating one type of recession (nearly identical in both years, modified slightly by additional recharge during the dry period in one of the two years, and modified significantly by continuous recharge during the dry period).

Atrazine concentration
The mean groundwater transit times were estimated from tritium measurements using the EPM, and were found to be approximately equal to 15 yr.The predicted atrazine concentration in spring water and the reconstructed history of atrazine input to the groundwater are shown exemplarily for the spring K17 in Fig. 2. The model tritium 2 displays a quicker reaction to change than model tritium 1, reflecting the larger piston-flow component of the former (η = 5.9 and 4.4 respectively).The end of atrazine leaching to the groundwater corresponding to the observed decrease in spring water is predicted to have taken place in 1996, which is consistent with the end of the first application phase discussed in Sect.2.3.A comparison of the different models predicting atrazine concentration over time is shown for the same three springs as above in Figs. 3, 4 and 5.For the models fitted using the recession, the shortest t em value of each spring (Table 1) was used for parameter estimation to reduce the disturbance in the signal due to additional recharge during the wet period.The uncertainty interval shown is calculated from the uncertainty of ±10 % in the proportion of arable land present in the catchment.The selected springs illustrate three different cases.In the first case (Fig. 3), η can be estimated uniquely from tritium data, and model predictions made using either tritium only or tritium in combination with the recession  curve are nearly identical.In the second case (Fig. 4), only one of the two tritium models predicts correctly the atrazine decrease.The η value of this model is also close to the estimate made from the recession coefficient.For the third example (spring K2), atrazine concentration was still within the range observed in the period 2008-2010.Consequently, because the information content is poor, and although η estimates from the recession and one tritium model are close to one another (2.4 versus 2.9), none of the models performs clearly better in predicting the atrazine evolution over time.

Discussion
The purpose of the study presented here was to combine information concerning mean tracer transit times and mean groundwater transit times calculated from spring recession into a consistent model.As could be shown, this approach is useful to reduce the number of possible models obtained from tritium measurements by rejecting those that do not agree with independent observations of spring baseflow recession.
In the case study presented here, continuous recharge during the recession period constitutes the main limitation of the method, as this additional inflow modifies the pure recession that would be observed if the groundwater reservoir were only emptying itself.Such an influenced recession may explain some of the variations observed in the discharge of K17 in 2009 for instance (Fig. 1,bottom).Another interesting result is that transit times calculated from the recession are with two exceptions longer in 2008 than in 2009 (Table 1).Since 2008 was a wetter year than 2009, it most probably indicates that significant recharge took place during the recession period in 2008 (and maybe to a lesser extent in 2009 as well, as illustrated by the large fluctuations in the discharge of K17 mentioned above).Springs larger than those draining the Luxembourg Sandstone or situated in more arid climates may be less sensitive to a small recharge rate during the summer months, as for instance in the ophiolite hard-rock aquifer presented by Dewandel et al. (2003).
A pragmatic solution is to inspect first the shape of the recession and look for irregularities not explainable by measurement error, and secondly to adopt the steeper recession observed for a spring as the closest approximation to a perfect uninfluenced recession (which may never be observed in shallow aquifers in temperate climates).Deviations due to fast preferential infiltration and limited in time are less problematic than continuous recharge, since fast-flowing water might not reach the groundwater table but travel laterally downslope as interflow without modifying the hydraulic gradient in the saturated zone, and thus lead to a temporary increase in the discharge without shifting the entire recession limb upwards.
Applying model-based corrections may constitute a workable alternative, with the effective recharge calculated from a bucket model for instance serving as input for a simple reservoir model parameterized using Eqs.( 5) and (6b).The fit could thus be performed on the entire time series, and not just on the recessions, a particular advantage for longer records.By that means it would be possible to include recharge dynamics into the model, which would make it more comprehensive, but also more complex.In the present case however, the agreement between the recession model and one of the tritium models show that, for most springs, the recession constant calculated from the data probably comes close to its pure value, and warrants keeping the analysis as simple as can be supported by available data.
The agreement between η estimates from recession and tritium also demonstrates that the initial assumptions of negligible matrix porosity and minimum volume still hold.This means in particular that, on the study site, the Luxembourg Sandstone is a good approximation of a simple porous medium where both water flow and water storage take place in the fracture network that nearly completely accounts for the total porosity of the formation.This agrees with previous estimates made from tritium measurements alone by Farlin et al. (2013).In groundwater systems where one of these assumptions is not respected however (either because diffusion into the rock matrix cannot be neglected or because of the presence of groundwater stored below the outlet elevation), the comparison between mean transit times in the saturated zone calculated from tritium and from the recession can allow the estimation of both the ratio of mobile to immobile volume and the ratio of water stored above the outlet to the total volume stored in the aquifer.

Conclusions
We have shown here that a quantitative analysis of aquifer baseflow recession can constrain the parameter range of transit time distribution function of an aquifer, contributing to reduce the prediction uncertainties of a model describing atrazine evolution in spring water over time.Simple to measure and thus often readily available even for scantily monitored sites, spring discharge can provide extremely useful information on the integrated hydraulic response of the aquifer to recharge dynamics.Groundwater systems dominated by slow flow are particularly suited, since modification of the pure baseflow recession caused by fast flow is minimal.

Fig. 1 .
Fig. 1.Discharge measurements and best fit of the quadratic recession.

Fig. 2 .
Fig. 2. (a) Predicted atrazine leaching from the soil and concentration time series in spring water (K17).The duration of the leaching period was adjusted to agree with the atrazine decrease observed in the springs (Figs. 3, 4 and 5).Observed concentrations between 2008 and 2010 are shown in the insert.(b) Effect of atrazine degradation on the breakthrough curve for half-lives of 5, 10 and 30 yr (model 1).All concentrations are normalised by the soil leaching concentration.

Fig. 3 .
Fig. 3. Model prediction of atrazine concentration (top graph) and model error as function of η (bottom graph) for K1.The models with t em estimated from tritium only or from the recession curve are nearly identical.η predicted from the discharge recession is shown by a red circle, very close to the minimum (certainly the global minimum, since η values higher than 5.5 are physically unlikely).

Fig. 4 .
Fig. 4. Model prediction of atrazine concentration (top graph) and model error as function of η (bottom graph) for K17.η predicted by tritium model 1 is shown by a dark blue circle, by tritium model 2 by a light blue circle, and from the discharge recession by a red circle.Based on the local minima of the error function, two models are nearly equally likely (η = 3.7 and η = 5.2).The first model however is closer to η predicted by the recession in 2008, and agrees better with observations.The recession observed in 2009 (recession 2) was obviously influenced by additional recharge, and lead to an estimate for η (1.5) much lower than estimates gained from tritium data.

Fig. 5 .
Fig. 5. Model prediction of atrazine concentration (top graph) and model error (bottom graph) for K2.η predicted by tritium model 1 is shown by a dark blue circle, by tritium model 2 by a light blue circle, and from the discharge recession by a red circle.As for K17 (Fig.3) two tritium models are equally likely (η = 3.7 and η = 5.2) with the first model closer to η predicted by the recession.In spite of this, none of the models is clearly better at predicting the atrazine evolution over time.

Table 1 .
t em calculated from Eq. (8) for the recession constants k estimated from the observed recessions in 2008 and 2009.