Using environmental tracers to determine the relative importance of travel times in the unsaturated and saturated zones for the delay of nitrate reduction measures

Groundwater quality in many regions with intense agriculture has deteriorated due to the leaching of nitrate and other agricultural pollutants. Modified agricultural practices can reduce the input of nitrate to groundwater bodies, but it is crucial to determine the time span over which these measures become effective at reducing nitrate levels in pumping wells. Such estimates can be obtained from hydrogeological modeling or lumped-parameter models (LPM) in combination with environmental tracer data. Two challenges in such tracer-based estimates are (i) accounting for the different modes of transport in the unsaturated zone (USZ), and (ii) assessing uncertainties. Here we extend a recently published Bayesian inference scheme for simple LPMs to include an explicit USZ model and apply it to the Dünnerngäu aquifer, Switzerland. Compared to a previous estimate of travel times in the aquifer based on a 2D hydrogeological model, our approach provides a more accurate assessment of the dynamics of nitrate concentrations in the aquifer. We find that including tracer measurements (H/He, Kr, Ar, He) reduces uncertainty in nitrate predictions if nitrate time series at wells are not available or short, but does not necessarily lead to better predictions if long nitrate time series are available. Additionally, the combination of tracer data with nitrate time series allows for a separation of the travel times in the unsaturated and saturated zone. 2018 The Authors. Published by Elsevier B.V. This is anopenaccess article under the CCBY license (http:// creativecommons.org/licenses/by/4.0/).


Introduction
Rising nitrate concentrations are negatively impacting the quality of drinking water around the world. In some cases, production wells are ultimately abandoned (e.g. Alikhani et al., 2016); in others, agricultural management practices are changed to reduce the input of nitrate (e.g. Bailey et al., 2015;Hansen et al., 2011;Wassenaar et al., 2006). Due to the costs involved with such nitrate reduction measures, it is important to assess their effectiveness. Such assessments entail not only an estimation of the nitrate reduction potential, but also an estimation of the time lag between implementation of the measures and the decrease of nitrate concentrations in the pumping wells. This time lag is determined mainly by the travel time of water from the point(s) of recharge to the pumping well (Böhlke, 2002;Böhlke and Denver, 1995). Under ideal conditions, this time lag can be determined by comparing observed nitrate concentrations in wells and historical inputs of nitrate into groundwater. However, measured nitrate time series in wells are often too short, the nitrate input history is not known well enough, or other processes affect the nitrate concentrations, such as denitrification (see e.g. Alikhani et al., 2016) or dilution with water that is low in nitrate (e.g. Baillieux et al., 2014). Thus, in most cases, travel times of groundwater, as well as these additional processes affecting nitrate are determined by means of hydrogeological modeling (2D or 3D) and particle tracking, environmental tracers, or a combination of both (Alikhani et al., 2016;Green et al., 2008;Jeffrey Starn et al., 2014;Kaown et al., 2009;MacDonald et al., 2003;McMahon et al., 2006;Osenbrück et al., 2006;Visser et al., 2009;Wang et al., 2012;Zoellmann et al., 2001). Environmental tracers are a complementary tool to hydrogeological models and time series of nitrate or other tracers for determining travel time distributions (or age distributions) of groundwater and dilution (Turnadge and Smerdon, 2014;Visser et al., 2009). Measurements of d 18 O, SF 6 , 3 H, 85 Kr, CFCs, 39 Ar,14 C,or 4 He are commonly interpreted with lumped-parameter models (LPM) to obtain travel time distributions: The concentrations of the tracers in a well are calculated as the convolution integral of the input history of the tracers and a relatively simple shape of the age distribution (Zuber and Maloszewski, 2001). These shapes can be described by one to just a few parameters (fitted to best reproduce measured concentrations of the environmental tracers) and their pre-selection is often based on the hydrogeological characteristics of the studied area (Cook and Böhlke, 2000). Knowing the age distribution and the nitrate input history can help reconstruct nitrate concentrations in a given well in the past and predict how they will evolve in the future. Alternatively, the nitrate input history can be reconstructed from measured nitrate time series in combination with the age distribution (Visser et al., 2007). Two current difficulties for such assessments are that (i) there is a lack of environmental tracers that can reliably determine travel times of water in the unsaturated zone and (ii) very few studies to date have included a rigorous assessment of the uncertainties associated with the determined age distributions and predictions of nitrate concentrations (Alikhani et al., 2016).
In the unsaturated zone (USZ), nitrate is transported in the water phase, while most environmental tracers are predominately transported in the gas phase, with the exception of 2 H, 3 H, and 18 O, which are part of the water molecule itself. Because the two modes of transport have different velocities, travel times for nitrate and environmental tracers from the surface to the well are different, especially if the USZ is thick. Because diffusive transport of gaseous tracers in the USZ is typically faster than the advective transport in the water phase, the USZ is often neglected for gaseous tracers (Zoellmann et al., 2001), i.e. the travel time estimated with these tracers is assumed to represent the travel time in the saturated zone (SZ) only. However, for aquifers with a thick USZ, the travel time of water -and thus nitrate -in the USZ can be a significant portion of the total travel time to the well (Sousa et al., 2013;Wang et al., 2012) and several studies have pointed out that it is problematic to ignore thick USZs (>10 m, Johnston et al., 1998;Schwientek et al., 2009;Zoellmann et al., 2001). To our knowledge, a study by Osenbrück et al. (2006) is the only study that has considered the time lag of gaseous tracers in the USZ (for CFCs) in conjunction with predicting the evolution of nitrate concentrations, but they provided only a limited assessment of the uncertainty of their results.
Several sources of uncertainty need to be assessed when determining age distributions and predicting nitrate concentrations. These include but are not limited to (i) uncertainty in defining a suitable conceptual model of the system (Bredehoeft, 2005), (ii) uncertainty of the tracer measurements, and (iii) uncertainty in additional aquifer parameters that are needed to calibrate models or estimate the time lag in the USZ for example. While uncertainties of tracer measurements are often reported, the resulting uncertainties in estimated parameters, age distributions, mean ages, and predicted nitrate concentrations are rarely calculated rigorously. Approaches to deal with uncertainty typically include sensitivity calculations and calculations of how quickly the goodness of fit decreases for a change in a parameter value (e.g. Corcho Alvarado et al., 2007;Onnis, 2008;Schwientek et al., 2009) and comparing results for different LPMs or different tracers with each other (Corcho Alvarado et al., 2007;e.g. Johnston et al., 1998;Osenbrück et al., 2006). Green et al. (2014) assessed the bias due to model selection and found that, while the obtained age distributions can be very different for different models, they tend to produce similar mean ages and even more similar predictions of solute concentrations, implying that uncertainty associated with model selection is of minor importance for the prediction of the evolution of nitrate concentrations. A new Bayesian modeling approach to estimate age distributions and their uncertainties resulting from measurement errors and additional model parameters in a more rigorous way was introduced by Massoudieh et al. (2012) and applied and extended in two follow-up studies of relatively simple systems (Alikhani et al., 2016;Massoudieh et al., 2014).
In this study, we measured multiple environmental tracers and nitrate to estimate age distributions and predict nitrate concentrations with the Bayesian statistical framework of Alikhani et al. (2016) in an aquifer with a thick unsaturated zone. Multiple pumping wells and piezometers in the Dünnerngäu, Solothurn, Switzerland, were sampled to determine how quickly already implemented nitrate reduction measures will result in a trend reversal of nitrate concentrations. Multiple tracers were used to characterize the groundwater age distribution on a time scale of a few years to decades: 3 H, tritiogenic 3 He (Poreda et al., 1988;Schlosser et al., 1988), radiogenic 4 He (e.g. Torgersen and Stute, 2013), 85 Kr (Althaus et al., 2009;Loosli et al., 2000;Smethie et al., 1992) and 39 Ar (Loosli, 1983;Loosli et al., 2000). The Bayesian statistical framework of Alikhani et al. (2016) allows us to estimate confidence intervals for the age distributions and predicted nitrate concentrations. Due to the thick USZ of the study area, a special focus is given to the role of the USZ, which is explicitly included in the LPM. Measurements in the USZ were performed to corroborate our model of the USZ. Results are discussed with respect to the effect that explicitly including the USZ in the LPM has on the estimated age distributions and modeled nitrate concentrations. Furthermore, we assess the relative importance of tracer measurements and nitrate time series to obtain reliable and accurate estimates. Finally, the combination of several environmental tracers, including 3 H, with the nitrate time series enables us to estimate the travel times in the unsaturated and saturated zones separately.

Study area
The Dünnerngäu aquifer, located in the Canton of Solothurn, Switzerland, is an important source of drinking water in an intensely farmed region of the Central Lowlands of Switzerland affected by high nitrate concentrations. The main crops are grains, corn, potatoes, vegetables, and hay. Extensive application of mineral fertilizer and manure led to increasing nitrate concentrations in leachate and ultimately in groundwater in the 80s and 90s. Nitrate concentrations today are generally above the environmental quality target of 25 mg/L set by the Swiss legislation and even close to the Swiss tolerance value for drinking water of 40 mg/L at some of the wells (Fig. 1). In 1999, the Canton of Solothurn initiated a program of nitrate reduction measures in order to meet the water quality targets for drinking water; primarily by adapted cultivation techniques and secondarily by transforming cropland to extensive meadows. The area of cropland farmed with adapted cultivation techniques increased steadily until 2010 to approximately 890 ha and by then 120 ha of cropland had been converted to extensive meadows (see Fig. A1 in the Electronic Appendix; Hunkeler et al., 2015). Together, this constitutes 89% of the agricultural area in the project perimeter and 26% of the total area of the aquifer (41 km 2 ). Based on a 2-D particle-tracking model implemented at the start of the nitrate reduction measures in 2000, nitrate concentrations were expected to decrease within a few years of the changes in cultivation (Geotechnisches Institut/TK Consult, 1999). While there was a small initial decrease in nitrate concentrations in the pumping wells, concentrations remained constant or increased again after 2005 ( Fig. 1). This unexpected behavior could either be due to ineffective nitrate reduction measures or due to an underestimation of the travel times, possibly because the particle-tracking model disregarded the thick unsaturated zone.
The nitrate input into the groundwater ( Fig. 1) was estimated based on mineral fertilizer application rates in Switzerland (Cornaz et al., 2005), which were scaled to a maximum nitrate concentration of 60 mg/L (corresponding to $64 kg of nitrate-N/ha/yr, depending on the recharge rate assumed for conversion), similar to maximum concentrations that were previously observed in a similar study area in Switzerland (Baillieux et al., 2014). For the time before 1975 no data were available and an exponential increase from 10 mg/L (11 kg N/ha/yr) in 1910 to 37 mg/L (39 kg N/ha/yr) in 1975 was assumed. The average nitrate concentrations of recharge in agricultural parcels where land use changes had already been implemented for several years were quantified based on nitrogen balances and measurements of vadose zone nitrate concentrations. An average nitrate concentration of 35 mg/L (37 kg N/ha/yr) was obtained. After the start of the nitrate reduction program in 2000, nitrate input concentrations were thus assumed to decrease proportionally to the increase of the area with adapted cultivation techniques; reaching a concentration of 35 mg/L in 2010, when most of the farmers were participating in the program (see Fig. A1 in the Electronic Appendix). Nitrate concentrations are expected to decrease to 30 mg/L (32 kg N/ha/yr) by 2020 due to the conversion of cropland to grassland. The hills and mountains surrounding the Dünnerngäu are mostly covered by forests and grasslands and nitrate loads are thus considerably lower there.
Geologically, the Dünnerngäu aquifer is bounded by the Jura Mountains (limestone) in the north and the Gäu Hills (molasse) in the south (Fig. 2). It consists mostly of unconsolidated gravel deposited after the glacial retreat at the end of the Würm glacial. In the eastern part, the aquifer is in direct contact with limestone from the late Jurassic (Pasquier, 1986). In the western part, the aquifer and the limestone are separated by marly molasses from the Oligocene and Miocene (Chattien-Aquitanien) and clayey silt (Pasquier, 1986). In large parts of the aquifer, the gravel is overlain by 1-3 m of clayey silt but the gravel may contain significant proportions of finer material in other areas as well (Hunkeler et al., 2015).
The general groundwater flow direction is from Oensingen in the southwest towards Olten in the northeast along the main axis of the system (Fig. 2). The saturated zone (SZ) is typically between 30 and 55 m thick along the main axis and thins out towards the Jura Mountains and the Gäu Hills (Hunkeler et al., 2015). The unsaturated zone (USZ), in turn, decreases in thickness from the upstream part of the aquifer towards the downstream part (see Fig. 2). Water balance modeling showed that direct diffuse recharge by precipitation, infiltration from streams, and lateral inflows from the Jura Mountains and Gäu Hills each contribute roughly a third to the total groundwater recharge of around 1000 mm/yr (Geotechnisches Institut/TK Consult, 1999;Hunkeler et al., 2015;Pasquier, 2000). The lower nitrate concentrations in infiltrating stream water and lateral inflows dilute the high nitrate concentrations in the direct recharge. The lateral inflows of mountain block recharge (MBR) from the Jura Mountains might exhibit a wide range of travel times from fast, episodic shallow subsurface flows to slow flows along much longer and deeper flow paths with long travel times (Pasquier, 1986).

Sampling and analytical methods
Samples were collected over the course of two years, from September 2011 to October 2013. Compared to the long term average, groundwater levels were low in 2011 and rose to average levels in 2013. In the pumping wells, water was pumped from near the top of the water column with the installed pumps; in the piezometers, water was extracted with submersible pumps (MP1). The wells were flushed several times before sampling and field parameters (pH, temperature, electrical conductivity, and oxygen saturation) were monitored with multi-parameter portable meters (HACH HQ30d). In total, 6 pumping wells and 50 shallow piezometers were sampled, most of them multiple times. However, this paper focuses on the six pumping wells characterized in Table 1.
Water samples for major ion measurements (including NO 3 À ) were filtered in the field and the cation samples acidified with concentrated nitric acid (10% HNO 3 ) to stabilize them. Subsequently, samples were measured by ion chromatography (DIONEX DX-120) at the University of Neuchâtel. For the analysis of HCO 3 À , the samples were titrated with 0.1 N HCl. The analytical uncertainty of the chemical analyses is ±5% and the error in the charge balance was <5% for all samples. The d 15 N and d 18 O stable isotopes of NO 3 À were analyzed by Isotope Ratio Mass Spectrometry (IRMS) at the University of Waterloo (Canada).
For 85 Kr and 39 Ar measurements, 500-2500 L of water were degassed in the field. From the extracted gas, Kr and Ar were isolated by preparative gas chromatography (Lu et al., 2014;Purtschert et al., 2013;Riedmann and Purtschert, 2016;Yokochi, 2016). The isotopic abundances of 85 Kr and 39 Ar were then determined by low-level gas proportional counting at the University of Bern (Loosli, 1983;Loosli et al., 1986;Loosli and Purtschert, 2005). Results for 85 Kr are reported in the unit of dpm/cc Kr , which stands for the number of 85 Kr decays per minute per mL (STP) of Kr gas. Activities of 39 Ar are given in %modern Ar, relative to the atmospheric activity concentration of 1.78 mBq/L Ar (Loosli et al., 2000). Typical relative measurement uncertainties are 5-10% for 39 Ar and <5% for 85 Kr.
Samples for dissolved noble gas analyses ( 3 He, 4 He, 20 Ne) were collected in copper tubes (Aeschbach-Hertig and Beyerle et al., 2000) and analyzed by mass spectrometry at the University of Bremen in Germany (Sültenfub et al., 2009). The analytical error is typically better than 1% for the noble gas concentrations and better than 0.5% for the He isotopic ratio. Tritium was electrolytically enriched and measured by liquid scintillation counting (LSC) at the Hydroisotop GmbH (Germany) on 1L of water with an uncertainty of typically 1 TU (Tritium Unit).
To characterize the water content and nitrate concentrations in the USZ, six vertical soil cores of up to 10 m depth were taken and analyzed. The cores were collected by rotary drilling in double core barrels without using water to avoid modification of the pore water composition. The cores were extruded from the inner core barrel with a piston with minimal disturbance. Samples of $100 g were immediately placed in air-tight centrifuge containers. Then, a known amount of water was added to the filled containers for nitrate extraction. The containers were placed on a shaker for 24 h before they were centrifuged to separate the solid and the liquid phases. The liquid phase was filtered and analyzed for NO 3 À by ion chromatography. These measurements were then corrected for the dilution by the water added for the extraction. The wet weight of the soil was compared to the weight of the dried solid phase to determine the gravimetric water content of the soil. Two additional multi-level boreholes were drilled in the USZ near soil core locations to extract soil gas for 85 Kr analyses (see Fig. 2). These were installed using the direct-push technique of the sealed multiport sampling system (SMPS, Ducommun et al., 2013). Thereby, the individual levels are isolated from each other by an injected mixture of cement and bentonite. Gas tightness to the atmosphere and between the screen levels was verified by measuring CO 2 and CH 4 concentrations of the extracted air. The latter was injected at one level and monitored in the adjacent levels.

Estimation of travel times
In lumped parameter models (LPM), the observed concentration c i (t) is related to the input concentration c i,0 (t) of tracer i by means of a simplified analytical function, q i (t), that represents the travel time distribution of the groundwater (e.g. Zuber and Maloszewski, 2001): where the last term represents a first order reaction such as radioactive decay (for a stable tracer, the decay constant k i is set to zero). Travel time distributions are characterized by one or several parameters such as the mean travel time or a measure of dispersive mixing. When solving the inverse problem of estimating these parameters (and others), given measured tracer concentrations in groundwater and the tracer input history, it is desirable to propagate uncertainties from, for example, the concentration measurements to obtain uncertainties for the estimated parameters and travel time distributions.

Tracer input functions
Input concentrations c i,0 (t) for 85 Kr are based on measurements from the BfS (Bundesamt für Strahlenschutz) station in Freiburg i. Br. (Winger et al., 2005). For 39 Ar, a constant atmospheric concentration was assumed (Loosli, 1983). The 3 H input is based on data from the closest monitoring station of the global network for isotopes in precipitation (GNIP) -Constance -(IAEA/WMO, 2014). Local emissions of 3 H from the watch-making industry in the studied area (Althaus et al., 2009) were considered with a constant offset from the Constance values. The best match between measurements and model was obtained for an offset of +7 TU. For pre-bomb waters (before 1953), a constant 3 H concentration of 5 TU was assumed (Roether, 1967).
Nitrate can be included in the model as an additional tracer sensitive to the travel time of water in the USZ, if denitrification is negligible or can otherwise be quantified, providing an additional constraint (Alikhani et al., 2016). The nitrate-rich water from intensely farmed areas is diluted with water low in nitrate (typically $10 mg/L in the study area) originating from streams or recharge from the surrounding hillslopes. To account for this, the nitrate input function (see Fig. 1) is scaled by a dilution factor d. A side effect of that dilution factor is also, that an error in the assumed maximum of the nitrate input history will mainly affect the dilution factor, but not other parameters, as long as the shape of the input history is approximately correct. For the old component (see below), a constant nitrate input of 10 mg/L is assumed.

Travel time distributions
Conceptually, the study area is divided into three compartments as illustrated in Fig. 3: The saturated zone (SZ), the unsaturated zone (USZ), and a compartment contributing old groundwater from below. Simple and commonly used LPMs such as the Exponential Model (EM) represent the saturated zone compartment only. While the model used here does explicitly include the other two compartments, LPMs require a high degree of simplification and idealization for each compartment. All three compartments, as well as recharge to the unsaturated and saturated zone are assumed to be homogeneous. The thickness of the saturated zone is assumed to be constant. In contrast, the USZ is represented by a wedge shape with linearly decreasing thickness from 30 m to 0 m along the main axis of the aquifer as a simple approximation of the decreasing thickness of the USZ along the groundwater flow direction in the study area, (see Fig. 2). The model only considers the vertical dimension (depth z) and one horizontal dimension (distance x along main axis of the aquifer, see Fig. 2). The upstream boundary is a no-flow boundary (there is a groundwater divide south-west of well 1). This simplified conceptual model does not distinguish explicitly between recharge from streams, through agricultural fields, or adjacent hill slopes, as this would result in an overly complex LPM relative to the available data. The following paragraphs explain in detail how the three compartments were included in the mathematical model of the residence time distributions.
The saturated zone in isolation meets the conditions of the exponential model first proposed by Vogel (1967), resulting in an exponential age distribution for wells screened over the whole aquifer thickness. For partially screened wells, as is the case in this study, the youngest and/or the oldest water components are not extracted (partial exponential model, PEM) . The cut-off age(s) are a function of the (known) depth of the screened interval (see Table 1 for well screens). For the EM and the PEM, the travel time of groundwater entering the well at depth z below the water table is given by Cook and Böhlke (2000) as where T 0 is the mean residence time, and Z is the thickness of the saturated zone. Furthermore, we relate the depth z below the water table to the horizontal position x where the flow path originated (see Fig. 3) by (Cook and Böhlke, 2000): where Z j is the total saturated thickness and X j is the distance of well j from the upstream boundary of the aquifer. This allows us to express the travel time along individual flow paths in the idealized saturated zone as a function of x, which will be useful for integrating the unsaturated zone. Moreover, Eq.
(3) also defines the contributing area (along distance x) of each well given its screen depth. There are partial overlaps between the contributing areas of several wells. However, the overlaps are smaller in reality than in the model (Geotechnisches Institut/TK Consult, 1999) because the wells are offset perpendicular to the main axis of the aquifer to various degrees. There is no comprehensive yet simple way of considering the USZ in the travel time distributions of LPMs. The USZ can either be incorporated through (apparent) travel time distributions which differ for different tracers, depending on their mode of transport in the USZ or, alternatively, the travel time can be defined as the time the water spends in the SZ only, correcting the tracer input functions to the water table. Cook and Solomon (1995) presented the governing differential equations for transport in both phases in the USZ and derived an analytical solution for exponentially increasing atmospheric concentrations of gaseous tracers under specific conditions. However, for the general case, no analytical solution of these equations is available. Consequently, simplified approximations have to be used (e.g. Osenbrück et al., 2006) or the equations have to be solved numerically (e.g. Corcho Alvarado et al., 2009Alvarado et al., , 2007Rueedi et al., 2005). All these studies with LPMs assumed a uniform thickness of the USZ for the entire study area (Corcho Alvarado et al., 2007;Rueedi et al., 2005) or for each well individually (Osenbrück et al., 2006). So far, only Schwientek et al. (2009) have considered a spatially variable thickness of the USZ of single wells, apart from 2D or 3D hydrogeological models, where this is more common (Onnis, 2008;Zoellmann et al., 2001).
Here, we present a new approach that allows for an USZ with variable thickness to be included in an extended LPM model. For parameter estimation, typically, Eq. (1) is integrated numerically. In contrast to this, and similar to Schwientek et al.'s (2009) moving particle approach, we numerically integrate a large number of flow paths over the screen depth. The total travel time for tracer i along a single flow line intersecting the well at depth z, t tot,i (z), is given by the sum of the tracer-specific travel time (or apparent time lag for gaseous tracers, further discussed below) in the unsaturated zone, t USZ,i (z), and the travel time t GW (z) in the saturated zone (see Eq. (2)). The tracer-specific travel time distributions can be obtained by aggregating all flow paths intersecting the contributing depth interval (well screen plus a margin of 4 m above and below for production wells). In addition to the flow paths through the USZ and SZ, slow mountain block recharge is included in the model as an admixed old component with proportion m and a piston-flow age T old (Fig. 3). Overall, the concentration c i (t) of tracer i at time t at the well can thus be calculated by: where N is the total number of flow paths f intersecting the contributing depth interval, z = z(f) is the (equidistantly-spaced) depth at which flow path f intersects the well, and k i is the radioactive decay constant of tracer i (k i = 0 for a stable tracer). One important factor affecting the travel time in the USZ, t USZ,i (z), is the spatially variable thickness of the USZ, h(x). Using Eq. (3), h(x) can easily be written as h(x(z)). The wedge shape of the USZ results in a wider age distribution of the water in the wells than would be obtained for a USZ of constant thickness, especially for wells with a long screen. For tracers transported in the aqueous phase of the USZ, a dispersion model (DM) with mean travel time t USZ,w (h(x(z))) = h(x (z))/v USZ,w and dispersion parameter d is used with v USZ,w the mean seepage velocity (same for all flow paths).
Transport of gases in the USZ is modeled with a simplified analytical solution of the governing equations described in Cook and Solomon (1995). They derived an analytical solution for stable gaseous tracers with exponentially increasing atmospheric concentrations in the USZ. This solution can be adapted for tracers with a first-order decay and constant atmospheric concentrations to give the apparent time lag t USZ,g,i of tracer i at the water table relative to the surface: where h(x) is the thickness of the unsaturated zone, D 0 i the diffusion coefficient in air of tracer i, D eff i the effective diffusion coefficient, h g the gas filled porosity, and 0 < s g < 1 the tortuosity coefficient of the gas filled pore space. Because the diffusion coefficient is different for every tracer, the apparent time lag is also tracer-specific. Fortunately, the diffusion coefficients are well known and do not introduce additional parameters that need to be estimated. The assumption of a constant atmospheric concentration is correct for 39 Ar and approximately valid for 85 Kr during the last 10-20 years (Winger et al., 2005). For tritiogenic 3 He and the radiogenic 4 He, Eq. (5) is not applicable because they accumulate in the subsurface. Because of the large diffusion coefficient of He, 3 He and 4 He are assumed to be in equilibrium with the atmosphere throughout the USZ (Solomon and Cook, 2000) so that tritiogenic 3 He and radiogenic 4 He are zero at the water table.

Parameter and uncertainty estimation
In total, there are nine fit parameters for each pumping well (Table 2), which are constrained by two to six tracer measurements per well and the nitrate time series. All fit parameters are assumed to be log-normally distributed and temporally constant, that is, transient effects on residence times and the functioning of the flow system are not considered. The 4 He accumulation rate a, the dispersion parameter d, and the gas filled porosity in the unsaturated zone h g are chosen as global parameters (i.e. the same for all wells). Furthermore, the age of the old component T old was only considered for wells where 39 Ar or 4 He were measured (wells 1, 2, and 5). All other parameters are well-specific. Thus, in total there are 3 + 3 + 5*6 = 36 parameters which are constrained by 24 tracer measurements, and several hundred nitrate measurements. Without the nitrate data, this model is underdetermined without further constraints on the parameter space.
Additional constraints can be obtained from relationships between different parameters of the model or from a priori knowledge of the hydrogeological situation. Firstly, the seepage rate in the USZ (R USZ ) and the groundwater recharge rate in the SZ (R GW ) can be calculated from the model parameters: where h is the total porosity, h w is the water filled porosity in the USZ, and Z the thickness of the SZ. In a perfect system, R USZ and R GW would be exactly the same, but here they are only assumed to be approximately the same because of considerable uncertainties in some of the parameters in these equations, which are not fitted. Secondly, an empirical relationship exists between the tortuosity coefficient s g and the gas-filled porosity h g (Millington and Quirk, 1961) which is also assumed to be only approximately true. Thirdly, R GW is constrained to R 0 % 1000 ± 700 mm (1r) based on estimates of the recharge rate from previous studies (Hunkeler et al., 2015). The uncertainty of R 0 is assumed to be so large not only because these estimates were quite uncertain, but also because the recharge rate might vary for the contributing area of the different pumping wells and because the estimate of 1000 mm would be too high if the old component, which does not contribute to R GW in the model is a large fraction of the water originating from the Jura Mountains (30%). Overall, these relations add a total of 3*6 = 18 constraining equations. The overlap of some contributing areas was not used as an additional constraint on the well-specific parameters of the respective wells, because it is unclear how a partial overlap should be included (simply making all parameters global for wells with only a small overlap in contributing areas clearly does not seem appropriate).
To fit the parameters and estimate their uncertainty, the Metropolis-Hastings algorithm (MCMC, Metropolis et al., 1953) was used with the following acceptance probability of a new state x' given state x: where where r ji is the standard deviation of the observations. The first term represents how well the model tracer concentrations agree with the measured tracer concentrations. The second term represents the constraints of Eqs. (6) and (7) and the estimate of the recharge rate. How strict these constraints are followed is determined by the value of r R1 , r R , and r h (r ? 1 disregards the constraints completely, whereas r ? 0 enforces them perfectly).
Here, the r were chosen to reflect the uncertainty of a priori parameters which are part of the respective equations. To obtain a more balanced weighting of the different tracers, especially when combining tracers (few measurements) and nitrate time series (many measurements with relatively low measurement uncertainty), the uncertainty of 39 Ar measurements (few measurements with high uncertainty) was halved to increase their weight and the weight of nitrate measurements was decreased by a factor of 4 to 10, depending on the length of the time series. New proposals for the MCMC algorithm were generated in a similar way to Turner and Sederberg (2012). An initial burn-in period before reaching convergence was discarded. Six different scenarios were investigated to assess the value of tracers and nitrate time series and the importance of explicitly including the USZ for modeling and predicting nitrate concentrations (see Table 3). For the models without USZ, the nitrate and tracer input was directly applied at the water table (i.e. t USZ,i = 0) and the parameters only relevant for the USZ were neglected (s g , h g , v USZ,w , and d). Depending on the scenario, burn in periods of up to a hundred thousand samples were required, after which a hundred thousand samples were generated using the MCMC algorithm. For the scenario NU, we found that the number of fit parameters was too high and the correlation between them too Table 2 Parameters of the model. The prior (log-normally distributed) is given in terms of the geometric mean (GM) and the geometric standard deviation (GSD, Limpert et al., 2001 strong to give meaningful results. To reduce this problem, the mean seepage velocity in the USZ was also made a global parameter in this scenario.

Chemistry and stable isotopes
Signatures of major ion composition and stable isotopes of the water are similar for all pumping wells with the exception of well 1. For example, measured bicarbonate concentrations vary only slightly from 359 to 377 mg/L (Table 4). Only nitrate concentrations exhibit a larger relative variability ranging from 27.1 to 39.2 mg/L (Fig. 4). Little temporal variation in concentrations of major ions was observed over the period of the three sampling campaigns. Measured stable isotopes of the water are in the range found in local precipitation (Schotterer et al., 2000). Well 1 differs from the others by having lower nitrate and bicarbonate concentrations and a slightly more negative stable isotope signature. Measurements in the shallow piezometers near well 1 confirm this observation (see Fig. 4 and Fig. A2). The more negative isotope signature of the water from well 1 is likely related to a higher recharge altitude, because the well is located close to a transverse valley intersecting the Jura mountain chain and draining water from the Jura Mountains. This water is low in nitrate and bicarbonate, thus also explaining the observed differences in the chemical signature.
Nitrate scales approximately linearly with bicarbonate in pumping wells and shallow piezometers, indicating a mixing of at least two water types (Fig. A3 in the Electronic Appendix): (i) water high in nitrate and bicarbonate (>40 mg/L and > 400 mg/L, respectively), and (ii) water low in nitrate and bicarbonate concentrations (approximately 10 mg/L and 270 mg/L, respectively). Type (i) derives from diffuse recharge of local precipitation; water of type (ii) is contributed by unpolluted water from the Jura Mountains, recharge from the streams in the study area, or water from the transverse valley.
The measured stable isotope signature of nitrate (d 15 N and d 18 O, Table 4) agrees best with the isotopic signature of soil N (Kendall and Caldwell, 1998;see Fig. A4 in the Electronic Appendix). However, the nitrate could also derive from mineral fertilizer and manure, with a larger contribution from manure than from mineral fertilizer. The isotopes of the nitrate show no sign of denitrification with the exception of one shallow piezometer. This is in accordance with the oxic conditions, which prevail throughout the aquifer (see Table 4). Consequentially, nitrate is considered to show a conservative behavior in this study.
Six profiles in the USZ show large variations of water content and nitrate concentrations between profiles and as a function of depth. The profiles generally show a decreasing water content with depth, which is related to the progressively coarser material with depth. Depth-averaged water contents range from 10 to 19% wt. There is no clear trend in nitrate concentrations with depth, which could confirm the effectiveness of the nitrate reduction measures (for more details, see Hunkeler et al., 2015).  Table 3 The six scenarios that were modeled. Parameters not relevant for the primary fit target were also optimized in a separate second step.

Tracers of travel time
Measured activities of 85 Kr in the USZ decrease significantly with depth (Fig. 5) and are up to 28% lower than the modern atmospheric value (2011) of approximately 75 dpm/cc Kr . Concentration profiles of 85 Kr in the USZ were calculated for several values of the gas phase tortuosity coefficient, s g , assuming a constant porosity with depth, a fixed depth of the water table (see Fig. 5), and a no-flux boundary at the water table (Cook and Solomon, 1995).
Values of s g from 0.06 to 0.10 agree best with the measured concentrations (Fig. 5).
Measured concentrations at the wells are shown in Fig. 6 (see also   fraction of the water is younger than 30 years. Well 1 also contains a considerable fraction of old water (>60 years and thus 85 Kr and 3 H-free) as indicated by an 39 Ar activity of 77.3% modern and an elevated 4 He rad concentration. Simple piston-flow or exponential model ages calculated from 3 H-3 He are systematically smaller (on average by about 4 years) than the same model ages calculated from 85 Kr. This time difference is attributed to the effect of the USZ and the corresponding delay of 85 Kr at the water table, whereas 3 He trit and 4 He rad in the USZ are strongly diluted by the additional He in the gas phase so that they their accumulation in the USZ is negligible.

Concentrations of travel time tracers
The six scenarios differ in their ability to reproduce the measured tracer concentrations. Generally, scenarios only relying on tracer measurements and scenarios including the USZ perform better in this respect (Fig. 6). As expected, measured tracer concentrations are best reproduced by scenario TU (using tracer concentrations only and including the USZ): Modeled tracer concentrations are close to or within the 1r (67%) confidence interval of the measured concentrations (and vice versa) for all 85 Kr, 4 He, and 3 He measurements with the exception of 85 Kr in well 1. The model is not as successful in reproducing the 3 H and 39 Ar data. For 3 H, this might be due to spatially and temporally varying local emissions, which are not considered by the constant offset of +7 TU. For scenario TU, the estimated model uncertainty of tracer concentrations is generally 2-3 times smaller than the reported measurement uncertainty for all tracers.
Including nitrate as a primary fit objective (e.g. scenario NTU) reduces the relative weight of the tracer measurements, resulting in an increased discrepancy between modeled and measured tracer concentrations. When comparing scenario TU with scenario NTU, this is especially the case for 3 H (Fig. 6). Tritium is affected most because it is transported in exactly the same way as nitrate. In contrast, transport of the other tracers depends on additional parameters (such as the 4 He accumulation rate a), which can compensate for changes in the parameters that determine the transport of nitrate. However, when nitrate alone is the primary fit objective (scenario NU), this compensation mechanism is not sufficient anymore and discrepancies between measured and modeled concentrations increase for 3 H, 3 He, and 85 Kr (Fig. 6). Scenario NU also has 2-5 times larger uncertainties of the modeled tracers, indicating that only using nitrate as a fit objective has limited ability to constrain the transport of these tracers.
Not including the unsaturated zone in the model generally results in a poorer reproduction of the measured tracer  (Table A1). concentrations, especially for 3 H in scenarios NT and T compared to NTU and TU (Fig. 6). However, if only nitrate is used as a primary fit objective, measured tracer concentrations are sometimes better reproduced if the USZ is not included in the model (scenario N). This is attributed to the fact that the mean seepage velocity in scenario NU was used as a global parameter, whereas in scenario N it was a local one. The uncertainty of modeled tracer concentrations is relatively independent of whether the model includes an USZ or not.

Range and sensitivity of fit parameters
Uncertainties of estimated fit parameters are smaller than the range of the prior, indicating the usefulness of the data for parameter estimation. The range of obtained values for each fit parameter and for the recharge rate are given in Table 5. The complete results for each scenario and pumping well are given in the supplementary material (Table A2). Uncertainties are similar for scenarios NTU, NT, TU, and T, and larger for the nitrate only scenarios (NU and N), indicating that nitrate measurements alone are not able to constrain the parameters as much as the tracer data.
Estimated values of a fit parameter for an individual pumping well are often similar for all six scenarios. For example, wells 2 and 5 consistently have the longest mean groundwater residence time T 0 . Scenarios without the USZ generally result in mean groundwater residence times T 0 that are 2-5 years longer, in compensation for not considering the delay of the USZ.

Age distributions
All pumping wells contain a significant fraction of water that is only a few years old, and much of the water of the young component is less than 10 to 15 years old (Fig. 7). Exceptions are wells 2 and 5, which are not containing any water younger than 2 and 5 years, respectively. This is in agreement with the fact that in these two wells the distance between the water table and the top of the screen is largest ($20 m), whereas in wells 9 and 10, for example, it is only $ 5 m. The cumulative age distributions for the investigated pumping wells vary between the scenarios, as Fig. 7 shows. Only scenarios NT and T have almost identical age distributions and related uncertainties for all wells except well 10. Uncertainties of the cumulative age distributions are smallest for scenario NT, followed by scenario NTU and largest for scenarios NU and N. Part of this increased uncertainty when only using nitrate as a fit objective comes from a higher uncertainty in the fraction of the old water component m. The estimated fraction m also is considerably larger, when only fitting nitrate: 47-75% for well 1 and typically 10-23% for the other pumping wells, compared to 26-38% and typically <10%, respectively, for scenarios including tracers. The problem is that in scenarios NU and N, m and d are highly correlated because, contrary to the other four scenarios, m is not constrained by additional tracer measurements.
Total mean travel times (T tot , USZ and SZ combined) across all pumping wells and scenarios range from 0.6 to 25.4 years ( Table 6). For individual wells, values obtained for different scenarios show less spread but still vary by 5-16 years (see Table A3 in the Electronic Appendix). Inter-scenario variability is thus considerably larger than the uncertainty estimated for T tot (e.g. 2-4 years for scenario NTU, Table 6). Travel times for wells 2 and 5 tend to be larger when including the tracer data and smaller when including the nitrate data (Fig. 7). This pattern applies not only to total travel times, but also to travel times in the unsaturated zone for models where it is included. Disregarding the old component, 30-58% of the total travel time is spent in the unsaturated zone (Table 6). Including the old component, which does not have an unsaturated zone in the model, reduces the mean travel time in the unsaturated zone to 1-25% of the total mean travel time. With a few exceptions in wells 2 and 5, the relative proportions of travel time in the USZ for a specific pumping well (T USZ /T tot ) are similar for all scenarios ( Table 6). The apparent time lag of 85 Kr in the unsaturated zone varies from 25% to 100% of T USZ . The travel time of water through the USZ is thus larger than the time lag for gases.

Nitrate
Modeled nitrate concentrations for all scenarios agree similarly well with observed time series for all pumping wells except well 2 and 5 (Fig. 8). In well 2, the observed time series has two distinct concentration peaks and modeled concentrations for the nitrate only scenarios (NU and N) are aligned with the first peak, whereas modeled concentrations for the tracer only scenario TU are aligned with the second peak. Modeled concentrations for the scenarios including both nitrate and tracers have the peak concentration in the middle between these observed peaks. Well 5, which has a single peak, shows a similar shift to a later peak concentration for tracer-based scenarios. The degree to which the model can explain the measured data (R 2 values) depends strongly on the individual wells and ranges from a few % to over 80%. The strikingly low correlation coefficient in wells 2 and 10 for scenario NU compared to other scenarios is caused by making v USZ,w a global variable in that scenario, whereas it is a well-specific variable for the other scenarios.
While median modeled nitrate concentrations are similar for different scenarios, differences between the six scenarios are larger concerning the width of the credible intervals of the modeled concentrations and the modeled peak nitrate concentration, especially if the observed time series begins when nitrate concentrations are already declining (wells 6 and 10 in Fig. 8). Credible intervals are generally largest for the tracer-only scenarios and smallest for scenarios NTU and NT, which use both tracer and nitrate data, with the nitrate-only scenarios having intermediate credible intervals. However, for the pumping wells where the peak nitrate concentration is not part of the observed time series, nitrate-only scenarios Table 5 Mean values and range of estimated fit parameters and recharge rates summarized over all scenarios and pumping wells. Full details of the results are given in Table A2  have the largest credible intervals. Especially in well 6, there is also a large difference between the two nitrate-only scenarios with scenario NU and N estimating peak nitrate concentration of approximately 32 mg/L and 42 mg/L, relatively. Most scenarios predict that none of the pumping wells has crossed or will cross in the future the tolerance value for drinking water set by the Swiss authorities thanks to decreasing nitrate inputs.
All scenarios indicate that, once the nitrate reductions measures become fully effective, nitrate concentrations should fall below the legislative target value of 25 mg/L for all wells. For wells 6, 9, and 10, all scenarios predict similar timing, with concentrations falling below the target concentration in approximately 2018, 2020, and 2015, respectively. For wells 2 and 5, the tracer-only scenarios indicate a much slower response to the nitrate reduction measures than the other scenarios. In these wells, achieving the target value will not happen before 2030 and 2020, respectively, but could be delayed until 2050 for well 2. No matter how different the age distribution and thus the timing of the decline in nitrate concentrations, the concentration that can ultimately be reached once the nitrate measures are fully effective is similar for most scenarios, as it only depends on the estimated dilution factor and to a lesser extent on the fraction of old water m. An exception is scenario N, which predicts a higher final nitrate concentration for well 6 than the other scenarios (Fig. 8).

Discussion
In this study, we employed environmental tracers and nitrate time series to estimate travel time distributions of groundwater and predict the future evolution of nitrate concentrations. Modeled nitrate concentrations agree with observed nitrate time series for most pumping wells and scenarios. All scenarios that use tracer Fig. 7. Cumulative age distributions (median and 67% (1r) credible intervals) of all pumping wells for 5 of the scenarios. Scenario T is not shown because the age distributions are almost identical to those of scenario NT. Dashed lines show the median age distribution at the water table for scenarios with an explicit USZ.

Table 6
Overview of estimated mean travel times in the unsaturated and saturated zone in years. The upper half of the table shows statistics for the 6 scenarios (3 in the case of T USZ and T USZ /T tot ) where Min, Mean, and Max refer to the shortest, mean, and longest mean travel time of all 6 (3) scenarios. The lower half of the data can also reasonably well reproduce the tracer measurements (Fig. 6). By combining tracer data and nitrate time series with a LPM that explicitly considers the USZ, we are able to separate travel times in the SZ and USZ. On one hand, the low computational costs of LPMs allow for a thorough assessment of uncertainties associated with calculated travel time distributions and nitrate predictions. On the other hand, the simplifications inherent in the LPM approach can lead to questionable results in a few specific cases. The main implication for the Dünnerngäu aquifer is that total travel times of water and nitrate are longer than previously thought, mainly due to the delay in the USZ, which has been neglected in the particle tracking model used for the previous estimate (Geotechnisches Institut/TK Consult, 1999). For nitrate, including the travel time in the USZ doubles the total travel time for most wells (Table 6, disregarding the old component), resulting in a slower reaction to nitrate reduction measures. Consequentially, nitrate concentrations at some of the wells will exceed the target value until 2040. Considerable retardation of gaseous tracers in the USZ was validated by 85 Kr measurements in the USZ (Fig. 5). None of the scenarios can fully reproduce some of the multi-annual fluctuations of the nitrate time series (Fig. 8), indicating that other processes are also modifying nitrate concentrations in the pumping wells. These could range from climatic effects on the nitrate input function (dry years leading to an enrichment of nitrate in the soil column) to shifting groundwater flow paths and mixing under changing hydrological conditions. An assessment of such processes would require a more detailed and transient hydrogeological model (3D or at least 2D+USZ).

Reliability and suitability of the model
Lumped parameter models are very simplified representations of flow dynamics and processes in real aquifers. Therefore, some aspects that might be important cannot be included, which may lead to biased results. For example, the limitation to two dimensions means that lateral inflows of groundwater cannot be considered appropriately, which might partially explain the discrepancy between R USZ and R GW (see Table 5). These lateral inflows could also have infiltrated through an USZ with different characteristics than local direct infiltration. Furthermore, the model assumes that parameters such as porosity and recharge rates are homogeneous for the contributing area of each well. This results e.g. in the dispersion parameter not representing actual hydrodynamic dispersion, but rather locally varying recharge rates, porosity, or deviations of USZ thickness from the idealized model in Fig. 3.
Despite this note of caution, in the following, model estimates of field observables such as porosity are compared to alternative sources of information. In the saturated zone, R GW (not including the old component) varies spatially, ranging from 680 to 4570 mm/yr for well 1 (average of all six scenarios: 2810 mm/yr, scenario NTU: 2760 mm/yr) to 310-680 mm/yr for well 2 (average: 460 mm/yr, scenario NTU: 480 mm/yr, Table 5). Well 1 is located close to the upstream boundary of the aquifer where the assumption of a constant thickness of the aquifer is likely not met. Additionally, the chemistry and stable isotopic signature suggest that a sizeable proportion of the water in well 1 derives from the mountains. For well 1, the model might thus not be completely suitable. Averaged over all other wells and all six scenarios, R GW is 860 mm/ yr. The highest and lowest recharge rates for a well often result from one of the nitrate-only scenarios (NU and N), whereas scenario NTU consistently results in recharge rates that are in the midrange of all scenarios. This indicates that the results of scenario NTU are more robust than those of scenarios relying on nitrate alone. Differences in recharge rates between wells might indicate different proportions of local recharge from precipitation versus lateral groundwater inflows from the hills and recharge from streams. However, values of some wells might also be biased due to simplifications of the lumped parameter model, such as assuming that the thickness of the saturated zone at the well is representative of the thickness of the saturated zone over the whole contributing area.
Seepage rates in the unsaturated zone (as darcy velocities) are generally smaller than R GW , with an overall average (excluding well 1) of 720 mm/yr (17% less than R GW ). Taking into account the average dilution factor of 0.53, these estimates suggest that direct local recharge from precipitation is on the order of 380-460 mm/yr. A direct comparison with the recharge rates estimated in previous studies (total recharge $1000 mm/yr, local recharge from direct precipitation $330 mm/yr) is difficult, because of the different spatial scales and because these studies distinguish local direct recharge from precipitation, infiltration from streams and lateral inflows from the mountains and hills surrounding the aquifer (Hunkeler et al., 2015;Paratte, 2013;Pasquier, 2000;TK Consult, 2001), whereas the current model only distinguishes between an old and a young component. Nonetheless, our results are of a similar magnitude as in these reports.
In the USZ, tortuosity coefficients and the partition of the porosity into water and gas phase also have a relative direct physical meaning. Fitted tortuosity coefficients for wells 2-10 and scenarios NTU and TU range from 0.08 to 0.18 with a typical uncertainty of 0.03 and an unweighted mean of 0.13 (Table 5). These values are close to the coefficients derived independently from the 85 Kr profiles (Fig. 5) of approximately 0.06-0.10. Scenario NU (nitrate only) results in considerably higher values and a large uncertainty. Well 1 has high tortuosity coefficients (0.22-0.46) for all three scenarios, which together with the highest v USZ of all pumping wells indicates that the thickness of the unsaturated zone in the contributing area of well 1 is smaller than the wedge shape suggests. This is in agreement with the location of the well at the upstream boundary of the system, where groundwater from the Klus of Balsthal, which has a thinner USZ, still constitutes a significant part of the whole water mass (see Fig. 2).
The measured soil water content of the two profiles in the USZ (Section 4.1) provides an estimate of h w (20%-30%) and together with Eq. (7) and the tortuosity values from the 85 Kr profiles in the USZ also of h g (15%-20%) resulting in a total porosity of 35%-50%. The modeled h g of 13%-19% are in a similar range. The modeled h w , in turn, are smaller ranging from 11% to 17% only, most likely because the total porosity h was fixed in the model to 30% ±6% in the USZ, which is lower than the total porosity indicated by the profile data (35%-50%). We did not adjust the porosity in the model for several reasons: Firstly, the profiles were taken from the top 8-10 m of the USZ and porosity generally decreases with increasing depth, resulting in a lower average porosity for the USZ as a whole. Secondly, the two profiles are only two points in space which do not necessarily represent the larger spatial average represented by the tracer based model results because soils with a high porosity often have a lower hydraulic conductivity (e.g. clays). Finally, any stagnant water, for example in local clay lenses, would be included in h w as derived from the USZ profiles, whereas the fitted h w represents the smaller effective porosity.
An important aspect not considered in our model is the effect of different travel times for nitrate-rich and nitrate-poor young water components. Comparing results of scenario TU with scenario NTU reveals that the performance of modeled tracer concentrations decreases most strongly for 3 H when including nitrate as a primary fit objective. This suggests that 3 H and nitrate represent a different groundwater mix. The dilution factor already implies the presence of at least two young water components, one nitrate-rich (local direct infiltration) and one nitrate-poor (e.g. lateral inflows or infiltration from streams). Tracer-based estimates refer to the travel time distribution of the whole water sample, whereas nitrate time series-based estimates are biased towards the nitrate-rich water component(s). An example of this effect can be seen at wells 2 and 5, where a comparison of scenarios NU, NTU, and TU imply a shorter mean travel time for the nitrate-rich component than for the nitrate-poor component (Fig. 7) and consequently a different timing of nitrate trend reversal (Fig. 8). Further evidence for the presence of a water mass that is young, but has a low nitrate concentration is that the dilution parameter is smaller than the fraction of young water (1-m) and that m is much larger for the nitrate only scenarios (17-26% for wells 2 and 5, compared to 3-5% for scenarios including tracers). An explicit modeling of these two water components would further increase the complexity of the model beyond a reasonable level for an LPM and the available calibration data.
Implementing a 3D hydrogeological model would help to overcome some of the limitations of the LPM approach in these hydrogeologically complex settings McCallum et al., 2015;Turnadge and Smerdon, 2014;Visser et al., 2009), but would introduce its own limitations again. Hydrogeological modeling would enable a proper incorporation of lateral inflows, spatially varying recharge rates, and aquifer properties and would also allow determining mixing between several young water components. However, hydrogeological models require not only more detailed data about aquifer geometry and properties, but also enough suitable data for calibration of unknown parameters. Especially for the USZ and the mixing of different water masses in the SZ, a calibration based on groundwater heads might not be sufficient, and additional tracer data or measurements from the USZ might be necessary for robust results (Schilling et al., 2017;Visser et al., 2009). Furthermore, Eberts et al. (2012) showed that if analogous conceptual models are used, LPMs and hydrogeological models often result in very similar contaminant predictions. Finally, due to the high computational costs involved in hydrogeological models, they are not very suitable for an estimate of the uncertainty as thorough as done here for the LPM approach.

The value of an explicit unsaturated zone, tracer measurements, and nitrate time series
An explicit representation of the USZ in the model allows for the separation of travel times in the USZ and the SZ at the cost of higher model complexity and thus higher demand on available measurements. Separating travel times in the USZ and SZ in a LPM is useful for examining processes that only take place in one of the zones, but also for comparing results with results from other models which modeled USZ and SZ separately (or only modeled one of them). In contrast to hydrogeological models, where omitting the USZ leads to a serious underestimation of the total residence time, not explicitly including the USZ in an LPM model calibrated on tracer measurements leads to a much smaller bias. Similar mean travel times are obtained for both models because the missing USZ is compensated for by a larger mean residence time T 0 in the SZ (Tables A2 and A3 in the Electronic Appendix). Although omitting the USZ in a tracer-calibrated LPM results in a similar total travel time, it results in bias in other parameters. For example, the recharge rate of the SZ obtained from T 0 will be underestimated if the USZ is not included in the model. If only nitrate measurements are used for model calibration, the parameters determining the travel times in the USZ and SZ (mainly v USZ,w and T 0 ) are highly correlated. Consequently, one of the fit parameters (v USZ,w ) had to be changed from well-specific to global in scenario NU to obtain meaningful results. As this caused unrealistic results for some of the wells, the explicit modeling of the USZ is not advisable if only nitrate time series are available for calibration.
Combining nitrate time series with measurements of gaseous tracers enables a separation of travel times in the USZ and the SZ. Gaseous environmental tracers are excellent tools to determine travel times in the SZ, but do not provide any information on travel times of water in the USZ. Because transport of gaseous tracers in the USZ is different from transport of nitrate (or 3 H), the strong correlation between v USZ,w and T 0 is alleviated. While time series of 3 H might be similarly useful as time series of nitrate, long nitrate time series are much more common. Combining gaseous tracers and nitrate time series thus creates a powerful tool to separate the travel times in the USZ and SZ.
The combination of environmental tracers and nitrate time series also improves parameter estimates and reduces the uncertainty of nitrate predictions for most wells. Scenarios that include tracer information generally result in physically more realistic parameter sets and reduced uncertainty, for example for the recharge rates, compared to scenarios which only use the nitrate time series to constrain the model. Regarding the prediction of the future evolution of nitrate concentrations, the benefit of including tracer data is largest for pumping wells where nitrate time series are short and do not include the peak of the nitrate concentration history (well 6 and 10 in Fig. 8). In this case, the tracer data reduce uncertainty of the peak nitrate concentration, the uncertainty of how nitrate concentrations will evolve, and the uncertainty of what nitrate concentration is finally achievable. For some wells with longer nitrate time series, including tracer data does not result in a significant improvement (well 2 and 5).
Our study also has implications for future assessments of nitrate reduction measures in other aquifers. Extending the statistical framework of Alikhani et al. (2016) to include an explicit USZ, we have shown that also in aquifers with a thick USZ, environmental tracers can be useful for obtaining better and more robust estimates of travel times and predictions of future nitrate concentrations. Moreover, combining these tracers with a time series of any water-borne compound that is not degraded in a specific aquifer (such as nitrate in this study) allows for a separation of the travel times of USZ and SZ. This could be useful information, especially if processes occurring in the USZ or SZ only are important. Furthermore, nitrate measurements in profiles from the USZ or samples from shallow piezometers would show much earlier whether nitrate reduction measures are effective. The advantage of the latter is that the water sample originates from a larger area, whereas the former will give a more plot-specific answer. However, a few occasional measurements in shallow piezometers or the USZ do not enable a reliable determination of the travel time from the surface to the piezometer for the water phase. Therefore, a combination of a few shallow measurements with the approach presented here to provide an estimate of travel times in the USZ seems to be a promising avenue. Finally, this study has shown once more that hydrogeological modeling without ground-truthing the results with suitable measurements can result in misleading conclusions, especially if the data used for model calibration mostly reflects transport in the SZ (groundwater head data), whereas the USZ also is a considerable and important part of the system.

Conclusion
Environmental tracer measurements and nitrate time series were used to estimate travel times of groundwater and predict the future evolution of nitrate concentrations in an aquifer polluted by nitrate. Special emphasis was given to the thick unsaturated zone that decreases in thickness along the flow direction. Model parameters were estimated using a Bayesian inference scheme that allows for an estimation of the joint error structure of these parameters. Our results show that the apparent underestimation of travel times in an earlier study based on a 2D hydrogeological model can be attributed to a disregard of the USZ in that study. By combining tracer data and nitrate time series, a separation of the travel times in the unsaturated and saturated zone becomes possible. Different scenarios of model complexity (with or without USZ) and data availability (only nitrate time series, only tracer measurements, or both) were evaluated. We find that including tracer measurements reduces uncertainty in nitrate predictions if nitrate time series are not available or short, but does not necessarily lead to better predictions if long nitrate time series are available. However, including tracer data generally results in physically more correct parameter estimates of the LPM.