# Assessment of the value of groundwater age time-series for characterizing complex steady-state flow systems using a Bayesian approach

Arash Massoudieh, Sarah Leray, Jean-Raynald de Dreuzy

► To cite this version:

Arash Massoudieh, Sarah Leray, Jean-Raynald de Dreuzy. Assessment of the value of groundwater age time-series for characterizing complex steady-state flow systems using a Bayesian approach. Applied Geochemistry, 2014, 50, pp.240-251. 10.1016/j.apgeochem.2013.10.006. insu-01083782

HAL Id: insu-01083782

https://insu.hal.science/insu-01083782

Submitted on 28 Apr 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Assessment of the value of groundwater age time-series for characterizing complex steady-state flow systems using a Bayesian approach

Arash Massoudieh<sup>a,*</sup>, Sarah Leray<sup>b,1</sup>, Jean-Raynald de Dreuzy<sup>b,c</sup>

In this work, the effectiveness of transient environmental tracer data in reducing the uncertainty associated with the inference of groundwater residence time distribution was evaluated. A Bayesian Markov Chain Monte Carlo method was used to infer the parameters of presumed residence time distribution forms—exponential and gamma—using concentrations of five tracers, including CFC-11, CFC-12, CFC-113, SF<sub>6</sub>, and <sup>85</sup>Kr. The transient tracer concentrations were synthetically generated using the residence time distributions obtained from a model of the Ploemeur aquifer in southern Brittany, France. Several measures of model adequacy, including Deviance Information Criteria, Bayes factors, and measures based on the deviation of inferred and true cumulative residence time distribution, were used to evaluate the value of groundwater age time-series. Neither of the presumed forms of residence time distributions, exponential and gamma, perfectly represent the simulated true distribution; therefore, the method was not able to show a definitive preference to one over the other in all cases. The results show that using multiple years of tracer data not only reduces the bias of inference (as defined by the difference between the expected value of a metric of inferred residence time distribution and the true value of the same met-ric), but also helps quantify the uncertainty more realistically. It was found that when one year of data is used, both models could almost perfectly reproduce the observed tracer data, even when the inferred residence time distributions differed substantially from the true one. When the number of years of tracer data is increased to four years, the uncertainty associated with the distribution parameters and the model structural uncertainty increased, as the presumed forms were not able to reproduce all the data accurately. This resulted in a more realistic assessment of model uncertainty due to structural error. It was also found that regardless of the prescribed age distribution form, the Bayesian method does a better job of capturing the cumulative ages at older ages; however, it is not able to reproduce the early ages well. The ability of the model to capture older ages improves as a greater number of years of tracer data is used, in cases of both presumed exponential and gamma distributions.

## 1. Introduction

Environmental tracers have been increasingly used for characterizing groundwater residence time distribution (RTD) (Leray et al., 2012; Massoudieh et al., 2012), and calibrating hydrogeological models (Portniaguine and Solomon, 1998; Reilly et al., 1994; Robertson and Cherry, 1989; Solomon et al., 1992). They characterize the underground flow patterns from recharge to discharge zones and provide an integrated insight into the system. Concentrations of CFCs, SF<sub>6</sub>, <sup>3</sup>H, <sup>39</sup>Ar, <sup>14</sup>C, He, and <sup>85</sup>Kr, for instance,

provide essential information over a broad range of temporal scales, from several months to several decades (IAEA, 2006; Newman et al., 2010). The interest in this type of data has promoted the quest for tracers with different atmospheric concentration histories and decay rates, in order to expand the range of ages that can be effectively characterized (Busenberg and Plummer, 2008; Newman et al., 2010; Sturchio et al., 2004). The ultimate goal is to characterize the full RTD beyond some of its moments, which requires more tracers than are often available. The set of environmental tracers that can be used for this purpose is limited. The tracers must be relatively stable with respect to the ages they are used to characterize, and they must not react or undergo adsorption into the solid matrix (Glynn and Plummer, 2005). They should be ideally inert tracers with low or no biodegradability or decay rates that are known with high confidence.

<footnote><sup>a</sup>Department of Civil Engineering, Catholic University of America, Washington, DC, USA</footnote>

<footnote><sup>b</sup>Géosciences Rennes (UMR 6118 CNRS), Université de Rennes 1, Campus de Beaulieu, 35042 Rennes Cedex, France</footnote>

<footnote><sup>c</sup>IDAEA (CSIC), c/Jordi Girona, 08034 Barcelona, Spain</footnote>

<footnote>* Corresponding author. Tel.: +1 (202) 3195671.</footnote>

<footnote>E-mail address: massoudieh@cua.edu (A. Massoudieh).</footnote>

<footnote n="1">¹ Present address: IFP Energies nouvelles, 1-4 Avenue de Bois-Préau, 92852 Rueil-Malmaison, France.</footnote>
Along with seeking additional potential environmental tracers, we assess here the effectiveness of the temporal variation of the concentrations of the commonly used ones. This variation might come from a transient evolution of the flow patterns, such as what occurs when starting a pumping (Leray et al., submitted for publication), and from the transient tracer atmospheric concentrations (Troldborg et al., 2008). The time series of environmental tracer concentrations have been used in a limited number of studies, and mainly for interpreting the field data (Troldborg et al., 2007, three-year period); (Long and Putnam, 2009, 16-year period); (Osenbrück et al., 2006, five-year period); (Sültenfuß et al., 2011, five-year period); (Bauer et al., 2001, three-year period). Manning et al. (2012) compiled a 13 year time series of multiple dissolved-gas age tracers from springs to constraint the form of RTD and concluded that the time series of age tracers is potentially a valuable tool for constraining age distribution in sampled water. In catchment hydrology, time-series of seasonally variable tracers in precipitation including ¹⁸O, ²H and ³H have been widely used to estimate transit time distributions using direct deconvolution and indirect methods based on lumped parameter models. A comprehensive review on different techniques used in estimation of transit time distributions in catchment hydrology can be found in McGuire and McDonnell (2006). Parametric and non-parametric methods have been used in the past to infer the transfer functions from tracer tests (Fienen et al., 2006; Luo et al., 2006). Because of inherent limitations of knowledge of the natural flow patterns, although the value of temporal chronicles has been emphasized, it has not been quantitatively evaluated. Temporal tracer data can both further corroborate and scrutinize previous inferences about the forms of age distribution. Temporal chronicles have been used to help identify the parameters of an *a priori* known form of the RTD (Ivey et al., 2008). The form of the distribution, however, is itself part of the identification problem. It has also been proposed that the transient variations of tracer age might help predict the general form of the RTD (Zhang, 2004). However, assessment of this idea is difficult, as it requires an independent knowledge of the flow patterns, which is out of reach in most natural sites.

In this study a Bayesian approach is used to evaluate the additional information content of tracer data collected at different times in comparison to the use of single samples obtained at a single time. The Bayesian approach allows us to quantify the uncertainty associated with estimation of the parameters of RTDs and therefore can provide an insight into how additional data help reduce these uncertainties. A synthetic aquifer case for which the flow patterns can be both quite complex and perfectly controlled is used for this study. It is common to test methodologies first on synthetic controlled cases to examine their relevance to natural cases. Steady-state flow conditions are considered and the capacity of the same environmental tracer, sampled at different times at the same position, to provide the needed information on the hydrosystem is evaluated. The steady-state assumption is acceptable over the long term if the overall withdrawal or recharge patterns at a site do not substantially change over the study period. Under a steady-state condition, the concentration of a tracer sampled at different dates is controlled solely by its past atmospheric chronicle (Leray et al., submitted for publication). If the historic atmospheric concentrations vary strongly enough and the sampling dates are sufficiently separated, samples of the single tracer taken at multiple times, in effect, can contain the equivalent information content of a number of different tracers. In reality, the same element taken at different times might not provide as much information as many unrelated tracers, due to the gradual changes in the atmospheric inputs, as well as the noise that results from observation error and other spatial and temporal variability. Nonetheless, it might still be valuable as a complementary tool for evaluating the relevance of classically used RTDs (Haitjema, 1995; Maloszewski and Zuber, 1982; Amin and Campana, 1995) and providing ways to cross-validate the inferred RTDs.

## 2. Methodology

Our analysis and methods are organized into four stages: construction of the hydrogeological model; simulation of solute transport and the derivation of the time-dependent environmental tracer concentrations; independent interpretation of these data to obtain estimated RTDs; and the ultimate comparison between the estimated and simulated RTDs. The tracers concentrations used for the analysis was obtained in one case from a hypothetical gamma RTD and in the second case from simulated RTD from a hydrogeological model of the Plømeur aquifer in North West France. The simulated tracer concentrations as a result of the hydrogeological model were treated as observed tracer data for inferring the RTD using several presumed lumped parameter models. Several measures of model sufficiency have been used to evaluate and compare the performance of the approach in determining the age distribution, as well as uncertainties associated with the estimated RTDs when different periods of tracer data are used for the RTD inference.

### 2.1. Simulated aquifer characteristics

The hydrogeological model is based on the Plømeur aquifer, a highly heterogeneous, hard-rock aquifer located on the south coast of Brittany (France). This aquifer has been exploited for water supply of the nearby city at a rate of $1.1 \times 10^6 \text{ m}^3/\text{year}$ since 1991. This aquifer offers a realistic hydrodynamic context for our methodological study, which is to assess the value of groundwater age time-series in inferring age distributions.

Based on the inversion of gravimetric data (Ruelleu et al., 2010) and observations (Touchard, 1999), the geological conceptual model of the Plømeur site is composed of two transmissive structures at the kilometric scale, which are embedded in less permeable micaschists. The first transmissive structure, named the contact zone, is a shallow-dipping fractured zone that constitutes the interface between an underlying leucogranite, Plømeur granite, and the overlying micaschists. The second structure is a North 20° normal fault. The almost impervious Plømeur granite forms the southern hydraulic barrier and part of the substratum. Another leucogranite, Guidel granite, forms the northern hydraulic barrier. The mean thickness of the system composed of the shallow-dipping structure and the micaschists has been estimated at 180–280 m (Ruelleu et al., 2010) (Fig. 1).

We will use the steady-state model developed by Leray et al. (2012) for which the transmissivity of the contact zone has been constrained both by hydraulic tests and by the mean piezometric level at the pumping well, while the hydraulic conductivity of the micaschists should not be too small to ensure that the area around the pumping does not seep. The hydraulic properties of the different rock units are summarized in Table 1.

The potential recharge rate, $R$, is assumed to be approximately 200 mm per year, based on previous estimations (Carn, 1990; Leray et al., 2012; Touchard, 1999). The supplying area to the pumping well located in the shallowest part of the aquifer then amounts to a few square kilometers and is limited to a North–South direction by the two granites. The long-term pumping well was found to induce a radial flow all along the fractured contact zone almost parallel to the dip plan of the contact zone because of its high conductivity contrast with the overlying micaschists. Flows through the micaschists are vertical near the surface and slightly bend perpendicular to the fractured zone at depth. As the goal of this paper is clearly methodological and without detailed information about the porosity field, we used an homogeneous porosity
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Conceptual model of the Plœmeur aquifer. The red lozenge is the pumping well where the steady-state residence time distribution is simulated. The North 20° normal fault does not outcrop and is represented by black stripes. Lengths are in meters. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<caption>**Table 1**</caption>

Parameters of the conceptual flow model based on the Plœmeur aquifer and chosen as the synthetic model.

<table><thead><tr><th></th><th>Synthetic model based on Plœmeur site</th><th>Source</th></tr></thead><tbody><tr><td>Structural parameter – mean thickness</td><td>180 m</td><td>Ruelleu et al. (2010)</td></tr><tr><td>Potential recharge rate R</td><td>200 mm/year = 6.3 × 10<sup>-9</sup> m/s</td><td>Touchard (1999)</td></tr><tr><td>Pumping rate Q<sub>w</sub></td><td>3.32 × 10<sup>-2</sup> m<sup>3</sup>/s</td><td>Ruelleu et al. (2010)</td></tr><tr><td>Hydraulic conductivities K</td><td></td><td></td></tr><tr><td>– Granite K<sub>G</sub></td><td>10<sup>-11</sup> m/s</td><td>Le Borgne et al. (2006)</td></tr><tr><td>– Micaschists K<sub>MS</sub></td><td>10<sup>-6</sup> m/s</td><td></td></tr><tr><td>Transmissivities T</td><td></td><td></td></tr><tr><td>– Contact zone T<sub>CZ</sub></td><td>2.27 × 10<sup>-3</sup> m<sup>2</sup>/s</td><td>Le Borgne et al. (2006)</td></tr><tr><td>– Normal fault T<sub>N20</sub></td><td>1.13 × 10<sup>-3</sup> m<sup>2</sup>/s</td><td></td></tr></tbody></table>

(1%). The calibration of a model would require a more complex representation of the porosity field in case of complex systems such as fractured ones (Cook et al., 2005). However, in our case, the non-trivial shape of the residence time distribution obtained from the synthetic model shows that the model complexity is sufficient to perform our analyses. Particularly, the form of RTD is determined by the aquifer structure. The smaller depth of the aquifer close to the pumping well, skews the RTD toward younger ages by shortening locally the flow paths Leray et al. (2012). The RTD is overall controlled by the superposition and interactions between different reservoirs and circulation scales; i.e., by local as much as global effects. Thus, we do not control the shape of the distribution, nor do we expect a classical reservoir model with an exponential or dispersion model with inverse-Gaussian distribution to reproduce these specific features. This particularly makes the synthetic example relevant to real cases.

## 2.2. Derivation of environmental tracer concentrations

Solute concentration is classically determined by solving the advection dispersion equation (Bear, 1991; de Marsily, 1986). In

this study, we disregarded local dispersion and diffusion, as we expected that mixing occurs mainly at the pumping well, where sampling is carried out, and is much larger than the mixing due to the diffusion-dispersion processes in this media (LaBolle and Fogg, 2001). The advection equation is solved backwards in time using a particle-tracking scheme (Neupauer and Wilson, 1999, 2001). 5 × 10⁶ particles are injected at the well according to a flow-weighted distribution. The RTD, $p(t)$, is obtained by simply looking at the distribution of the equally likely travel times that it takes for particles released at the well to arrive at the recharge locations. Mean concentration, $c_{i,j}$, which is the concentration of tracer $i$ in sample $j$ at the pumping well, is determined by the convolution of $p(t)$ with the atmospheric chronicle, $c_{in,i}(t_{w,j} - t)$ (Kreft and Zuber, 1978; Maloszewski and Zuber, 1982):

$$c_{i,j} = \int_{0}^{\infty} c_{in,i}(t_{w,j} - t) p(t) e^{-\lambda_i t} dt \quad (1)$$

where $c_{in,i}$ represents the atmospheric concentration of tracer $i$ and $\lambda_i$ is its decay rate. Eq. (1) shows the time dependency of the sampled concentration, which, under steady-state assumption, comes exclusively from the non-linear variation of the atmospheric concentration of the given environmental tracer (Leray et al., submitted for publication; Zhang, 2004).

The steps to infer the residence time distribution, $p(t)$, in this work include first assuming a functional form representing $p(t)$ [i.e., $p(t) = f(t; \phi_1, \phi_2, \dots)$], where $\phi_1, \phi_2, \dots$ are the parameters needed to fully determine the distribution, and then estimating the parameters $\phi_1, \phi_2, \dots$ using a statistical inverse modeling approach. Various sources of uncertainty associated with the inferred parameters include observation error due to measurement error and the spatial and temporal variability of tracer concentrations; model structural error (epistemic uncertainty) due to the fact that the actual RTD cannot be represented accurately by any presumed functional forms and the uncertainties associated with environmental factors, such as tracer atmospheric concentration records and parameters affecting the transport of tracers in the aquifer. Assuming that an RTD function with given parameters $\phi_1, \phi_2, \dots$ results in tracer concentrations $C = [c_{ij}]_{n \times m}$ (a matrix) and that the observed tracer concentrations can be represented as $\tilde{C} = [\tilde{c}_{ij}]_{n \times m}$, while assuming (1) no prior information is present
regarding the parameters defining age distribution (non-informative priors), (2) there are no uncertainties associated with the decay rate of the tracers, and (3) the tracers travel with pore water without undergoing any reactions or adsorption into the aquifer material, the Bayes theorem can be expressed as:

$$ \pi(\phi_1, \phi_2, \dots, \Gamma | \tilde{\mathbf{C}}) = \frac{\pi(\tilde{\mathbf{C}} | \phi_1, \phi_2, \dots, \Gamma)}{\pi(\tilde{\mathbf{C}})} \quad (2) $$

where $\Gamma$ is the variance-covariance matrix of the observation error (or a transformation of it), $\pi(\tilde{\mathbf{C}}|\phi_1, \phi_2, \dots, \Gamma)$ is the likelihood function, and $\pi(\tilde{\mathbf{C}})$ is a normalizing factor that is equal to the integral of the numerator over the entire parameter space. Using a log-normally distributed observation error structure, and assuming that the observation error variance for log-transformed observed concentrations, given the true concentrations of all the tracers, are all the same, the variance-covariance matrix can be written as a scalar value $\sigma$ multiplied by an identity matrix, and the likelihood function can be written as:

$$ \pi(\tilde{\mathbf{C}}|\phi_1, \phi_2, \sigma) = \frac{1}{\sqrt{2\pi^m \prod_{i=1}^m \prod_{j=1}^n \tilde{c}_{ij} \sigma^{mn}}} e^{-\sum_{j=1}^n \sum_{i=1}^m \frac{[\ln(\tilde{c}_{ij}) - \ln(c_{ij})]^2}{2\ln(\sigma^2)}} \quad (3) $$

Combining Eqs. (2) and (3) results in the following equation for the posterior distribution of the parameters defining the RTD:

$$ \pi(\phi_1, \phi_2, \dots, \Gamma | \tilde{\mathbf{C}}) \propto \frac{1}{\sigma^{mn}} e^{-\sum_{j=1}^{n} \sum_{i=1}^{m} \frac{[\ln(\tilde{c}_{ij}) - \ln(c_{ij})]^2}{2\sigma^2}} \quad (4) $$

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Historical equivalent atmospheric concentrations of the tracers used in this study.</figcaption></figure>

<table><caption>**Table 2**
Mean concentration $C_w$ for CFC-11, CFC-12, CFC-113, SF$_6$ and $^{85}$Kr at the pumping well.</caption><thead><tr><th>Sampling time</th><th></th><th>CFC-11</th><th>CFC-12</th><th>CFC-113</th><th>SF6</th><th><sup>85</sup>Kr</th></tr></thead><tbody><tr><td rowspan="2">2003</td><td>Mean concentration, C<sub>w</sub></td><td>249.74</td><td>501.07</td><td>74.33</td><td>3.66</td><td>51.68</td></tr><tr><td>Apparent age</td><td>14.89</td><td>12.18</td><td>12.76</td><td>7.13</td><td>6.07</td></tr><tr><td rowspan="2">2004</td><td>Mean concentration, C<sub>w</sub></td><td>250.27</td><td>506.07</td><td>75.02</td><td>3.88</td><td>52.44</td></tr><tr><td>Apparent age</td><td>15.84</td><td>12.75</td><td>13.88</td><td>7.26</td><td>6.46</td></tr><tr><td rowspan="2">2005</td><td>Mean concentration, C<sub>w</sub></td><td>250.49</td><td>510.44</td><td>75.51</td><td>4.11</td><td>53.13</td></tr><tr><td>Apparent age</td><td>16.82</td><td>13.40</td><td>14.76</td><td>7.40</td><td>6.82</td></tr><tr><td rowspan="2">2006</td><td>Mean concentration, C<sub>w</sub></td><td>250.27</td><td>513.78</td><td>75.85</td><td>4.34</td><td>53.66</td></tr><tr><td>Apparent age</td><td>17.84</td><td>14.21</td><td>15.68</td><td>7.42</td><td>7.11</td></tr><tr><td rowspan="2">2007</td><td>Mean concentration, C<sub>w</sub></td><td>249.92</td><td>516.86</td><td>76.04</td><td>4.56</td><td>54.15</td></tr><tr><td>Apparent age</td><td>18.88</td><td>15.04</td><td>16.64</td><td>7.42</td><td>7.02</td></tr><tr><td rowspan="2">2008</td><td>Mean concentration, C<sub>w</sub></td><td>249.06</td><td>518.97</td><td>76.09</td><td>4.78</td><td>54.48</td></tr><tr><td>Apparent age</td><td>19.96</td><td>15.77</td><td>17.62</td><td>7.46</td><td>6.93</td></tr><tr><td rowspan="2">2009</td><td>Mean concentration, C<sub>w</sub></td><td>248.15</td><td>520.65</td><td>76.02</td><td>5.00</td><td>54.84</td></tr><tr><td>Apparent age</td><td>21.03</td><td>16.46</td><td>18.64</td><td>7.55</td><td>6.83</td></tr><tr><td rowspan="2">2010</td><td>Mean concentration, C<sub>w</sub></td><td>246.85</td><td>521.58</td><td>75.84</td><td>5.22</td><td>55.03</td></tr><tr><td>Apparent age</td><td>22.11</td><td>17.27</td><td>19.68</td><td>7.67</td><td>6.78</td></tr></tbody></table>

A Markov Chain Monte Carlo (MCMC) method (Gamerman and Hedibert, 2006; Kaipio and Somersale, 2004) was used to draw a large number of samples from Eq. (4). A C++ code is written to perform the MCMC simulation. The number of Markov chains can be determined by the user. In the example application presented in the next section, 8 chains were used and 1,000,000 samples resulted in convergence of the MCMC method. The first 100,000 samples were left out as “burn-in” period. The criteria suggested by Geweke (1992) and Geweke and Tanizaki (2001) was used to evaluate the convergence of the MCMC algorithm. For more details on the Bayesian inference, see Massoudieh et al. (2012).

Fig. 2 shows the time evolution of $C_{in}$ for the environmental tracers used in this study: CFC-11, CFC-12, CFC-113, SF$_6$ and, $^{85}$Kr in the Northern Hemisphere atmosphere (IAEA, 2006). CFCs have been exhibiting an increasing atmospheric concentration since the late 1940s, which stabilized in the mid-1990s. SF$_6$ has been monotonically increasing since the late 1950s. $^{85}$Kr began to increase in 1980 and flattened out in 2000.

Table 2 synthesizes the mean concentration, $C_w$ of the five different tracers used in this study and apparent ages for each tracer at each year. These concentrations are obtained using RTD $p(t)$, produced by the chosen calibrated flow model (Table 1 and previous section) with a porosity of 1%. They constitute the data time series that covers the eight-year period from 2003 to 2010. This eight-year period is consistent with medium-term planning relative to site studies and projects duration. One interesting observation is that for the three tracers with highly non-linear atmospheric concentration trend in the recent past (i.e. CFC-11, CFC-12 and CFC-113) the apparent ages not only deviate from the known mean age but are also non-consistent within years while the apparent ages are more consistent for SF$_6$ and $^{85}$Kr with approximately linearly changing atmospheric concentrations in the recent past.

In order to make the synthetic tracer concentration resemble real-world observations more realistically, a noise that was produced based on a log-normally distributed structure and a standard deviation of 5% was added to each calculated concentration to represent the error due to measurement, heterogeneity, and model structural error.

## 2.3. Comparison of model performance and uncertainty

Two measures, Deviance Information Criteria (DIC) and Bayes factors, were used for evaluating the goodness of different RTD forms. DIC (Spiegelhalter et al., 2002), a measure that is used to compare the goodness of each model structure, takes into account how good a model structure can reproduce the observed data,
while explicitly imposing a penalty for the level of complexity of the model:

$$DIC = 2\overline{D}(\bar{\Theta}) - D(\bar{\Theta}) \qquad (5)$$

where $\overline{D}(\bar{\Theta})$ is the mean deviance and $\bar{\Theta}$ is the deviance of the mean, respectively defined as:

$$\overline{D}(\bar{\Theta}) = -2E[\ln p(\tilde{C}|\Theta)] - 2\log f(\tilde{C}) \qquad (6)$$

$$D(\bar{\Theta}) = -2\ln p(\tilde{C}|\bar{\Theta}) - 2\log f(\tilde{C}) \qquad (7)$$

where f is a standardizing term that is only a function of the observed data, $\Theta$ represents the parameters that are random-vector distributed based on the posterior distribution, and $\bar{\Theta}$ is the expected values of the posterior distribution of parameters. Because $2\log f(\tilde{C})$ is only a function of the observed data, it is not affected by the model structure, and therefore, it becomes irrelevant when comparing different forms of RTD functions using the same data. The DIC was easily calculated using the results of the MCMC calculations. A smaller DIC means a better model.

To compare different age distribution forms in terms of their ability to reproduce the measured data, the Bayes factor (Jeffreys, 1935; Kass and Raftery, 1995) was used. The Bayes factor for comparing two models, $M_1$ and $M_2$, assuming an equal prior probability for the models, is defined as:

$$B_{12} = \frac{\int_{\Theta_1} p(\tilde{C}|\Theta_1, M_1)p(\Theta_1)}{\int_{\Theta_2} p(\tilde{C}|\Theta_2, M_2)p(\Theta_2)} \qquad (8)$$

$B_{12}$ represents the ratio of the odds that model 1 is the true model to the odds that model 2 is the true model, $p(\tilde{C}|\Theta, M)$ is the likelihood of observing concentrations $\tilde{C}$ given model $M$, and $p(\Theta)$ is the prior density of the parameters. The odds of each model being the correct model when the only choices of model structures are models 1 through m can be calculated as:

$$I_k = \frac{\int_{\mathbf{X}_s} p(\tilde{C}|\Theta_s, M_k)p(\Theta_s)}{\sum_{k=1}^{m} \int_{\mathbf{X}_s} p(\tilde{C}|\Theta_s, M_k)p(\Theta_s)} \qquad (9)$$

where m is the total number of RTD distribution forms being evaluated. The integrals of posterior distribution are estimated using the Monte Carlo method Carlin and Chib, 1995:

$$\int_{\Theta_s} p(\tilde{C}|\mathbf{X}_s, M_s)p(\Theta_s) \approx \frac{1}{n} \sum_k p(\tilde{C}|\Theta_s^{(k)}, M_s)p(\Theta_s^{(k)}) \qquad (10)$$

# 3. Results

## 3.1. Hypothetical gamma residence time

To evaluate the effect of transient tracer information on reduction of the uncertainty associated with inferred RTDs when there is no model structural error, the method was applied first to a case where the concentrations of four tracers—CFC-11, CFC-12, CFC-113, and SF$_6$—were obtained based on a known gamma residence time distribution, using Eq. (1), for eight consecutive years, 2003–2010. The gamma distribution was considered because first it is relatively simple and is determined using only two parameters while it also has the flexibility to reproduce skewness and tailing observed in commonly seen breakthrough curves. A 5% randomly generated observation error based on a log-normal and multiplicative error structure was added to each calculated tracer concentration to represent measurement and sampling error and the error due to spatial variability. The Bayesian approach was then used to estimate the parameters defining a gamma, inverse Gaussian (dispersion model) and a log-normal distribution and for cases

where one, two, four, and eight years of tracer data were used. The gamma and inverse Gaussian distributions were used to evaluate the capability of multiple years of data to identify the right form of distribution by discriminating between the correct (gamma) and incorrect (inverse Gaussian and lognormal) distributions. Fig. 3 shows the presumed (true) RTD (solid line) and 60 samples of inferred plausible RTDs when one, two, four, and eight years of data were used, while assuming gamma and lognormal distributions. As can be seen, as the number of years of tracer data used for inferring the RTDs increases, the uncertainty associated with the distribution parameters decreases. In order to quantify the uncertainty associated with the inferred RTDs, a measure representing the spread of the posterior probability distribution of a scalar measure extracted from the RTD around the true value can be used:

$$\begin{aligned} \delta &= E[(\xi - \xi_t)^2]^{1/2} \\ &= \left( \int_{-\infty}^{\infty} [\xi(\phi_1, \phi_2, \dots) - \xi_t]^2 \pi(\phi_1, \phi_2, \dots) d\phi_1 d\phi_2, \dots \right)^{1/2} \end{aligned} \qquad (11)$$

where $\xi$ is a random variable representing the scalar measure extracted from the RTD, and $\xi_t$ is the true value of $\xi$ that is known if the true RTD is known. Any single scalar quantity extracted from the RTD (e.g., mean, median, cumulative distribution at a certain age) can be considered to represent $\xi$. Here, we define $\xi$ as the fraction of water younger than a certain age $L$; therefore, $\xi_L = \int_0^L p(t)dt$. In addition, the bias of estimation can be calculated as $b_L = E(\xi_L) - \xi_t$. Using the samples drawn from the posterior distribution obtained from the MCMC algorithm, the values of the uncertainty measure, $\delta$, and the bias measure, $b_L$, can be calculated easily.

Fig. 4 shows the 95% credible intervals of modeled concentrations of the four tracers vs. the “observed” concentrations for the two cases were two years and eight years of data used for parameter estimation of RTDs using presumed gamma and log-normal distributions. The results for the inverse-Gaussian distribution have not been shown to keep the figures more readable. As it can be seen, as the number of years used for parameter estimation increases, the ability of the inverse-Gaussian model to reproduce the data is declined compared to the gamma distribution. Table 3 shows Bayes factors, uncertainty measure $\delta_5 = E[(\xi_5 - \xi_{t,5})^2]^{1/2}$, and bias $b_5$. The subscript 5 indicates that these measures were calculated using the cumulative RTD at 5 years. A smaller $\delta_5$ means a smaller uncertainty is quantified for cumulative age distribution at 5 years for the three presumed RTD forms. $I_k$ is the odd of one model being the correct one if the only possible models where the ones considered (i.e. gamma, inverse-Gaussian and lognormal). When one and two years of tracer data are used, the method is not able to definitively identify the right model as the odds of each model to be the true is not significantly different. However, when four years of tracer data is used a 90% odd is assigned to the correct (gamma) model and when eight years of tracer data is used the correct model is definitively selected. As is expected when the correct form of RTD is assumed (gamma), using samples from a larger number of years of tracer data results in a smaller uncertainty and a smaller bias of estimation. However, when the correct form of RTD is not used, using a larger number of years does not result in an improvement in uncertainty and bias. This is an important conclusion, as in practical cases, the form of RTD is not known and there is a high probability that the simple forms selected do not match the real distribution. In these cases, the uncertainty and bias of inference do not improve when using samples from multiple years. The implication of this outcome is that additional data can affirm or dispute the presumed form of distribution. The Bayes factors show the odds of each model structure to be the
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Hypothetical gamma residence time distribution (solid black line) and inferred gamma distributions using the tracer data calculated from it when (a) one year, (b) two years, (c) four years and (d) eight years of tracer data were used for inferring the residence time distribution.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** "Observed" tracers concentrations and 95% credible intervals obtained from the posterior distribution of gamma and exponential RTDs based on 8 years (top row), and two years (bottom row) of tracer data. The tracer concentrations are calculated based on a hypothetical gamma distribution and corrupted with 5% noise.</figcaption></figure>

correct form when the only choices are the ones considered (in this case, gamma, inverse Gaussian and lognormal). What can be learned from this experiment is that when one form of RTD describes the tracer data substantially better than other forms do, using transient data can strengthen the confidence toward the right form and more definitively reject the inappropriate forms.

## 3.2. Simulated residence time distribution from the conceptual flow model based on the Plœmeur aquifer

Gamma and exponential models were applied to the tracer concentrations calculated using the simulated RTD for the Plœmeur aquifer. Here, the simulated RTD is deemed the “true” RTD, and
<table><caption>**Table 3**
Bayes factors, uncertainty and bias of RTD inference for the hypothetical case using gamma and exponential forms.</caption><thead><tr><th></th><th>1 yr</th><th>2 yrs</th><th>4 yrs</th><th>8 yrs</th></tr></thead><tbody><tr><td><em>Gamma</em></td><td></td><td></td><td></td><td></td></tr><tr><td>BF</td><td>38.7%</td><td>27.0%</td><td>90.3%</td><td>100.0%</td></tr><tr><td>d5</td><td>2.02E-04</td><td>4.74E-05</td><td>1.68E-05</td><td>2.62E-06</td></tr><tr><td>bias</td><td>-1.27E-02</td><td>-3.32E-03</td><td>-2.74E-03</td><td>-2.48E-03</td></tr><tr><td><em>IG</em></td><td></td><td></td><td></td><td></td></tr><tr><td>BF</td><td>24.0%</td><td>38.5%</td><td>6.6%</td><td>0.0%</td></tr><tr><td>d5</td><td>6.42E-03</td><td>1.94E-03</td><td>2.21E-03</td><td>3.40E-04</td></tr><tr><td>bias</td><td>-1.27E-02</td><td>4.34E-02</td><td>4.68E-02</td><td>-3.36E-02</td></tr><tr><td><em>Lognormal</em></td><td></td><td></td><td></td><td></td></tr><tr><td>BF</td><td>37.3%</td><td>34.5%</td><td>3.1%</td><td>0.0%</td></tr><tr><td>d5</td><td>2.06E-02</td><td>1.68E-02</td><td>1.81E-02</td><td>6.39E-03</td></tr><tr><td>bias</td><td>1.42E-01</td><td>1.29E-01</td><td>1.34E-01</td><td>1.46E-01</td></tr></tbody></table>

the performance of the method is evaluated based on how well the inferred RTDs based on the tracer information resemble the “true” one. The concentrations of five tracers – CFC-11, CFC-12, CFC-113, SF₆, and ⁸⁵Kr – were calculated for eight years, ending in 2010, as shown in Table 2. These calculated tracer concentrations are considered observed tracer concentrations and are referred to as such henceforth. All the tracers were considered to be conservative (non-decaying) except ⁸⁵Kr, for which a decay rate of 0.0641 yr⁻¹ (equivalent to a half-life of 10.8 years) was considered. In order to compare the effect of one additional tracer against using multiple years, the analysis was performed once using all five tracers and once using only four of the tracers, excluding ⁸⁵Kr. A random 5% log-normal and multiplicative error was added to each of the calculated tracer concentrations to represent measurement error and the effect of temporal and spatial heterogeneities on measured concentrations.

### 3.2.1. Adequacy of analytical age distributions

Fig. 5 shows the “true” RTD and 60 samples produced based on random realizations generated from the posterior distributions of the RTD parameters, shown using gray lines. The 60 realizations can be thought of as plausible RTDs inferred by the algorithm, resulting in the true observed tracer concentrations. It can be seen that neither the gamma nor exponential distribution is able to reproduce the RTD on the full time range. The correct estimate at older ages (those over three years) corresponds with a poor match of the early parts of the “true” RTD (ages less than three years). More precisely, the very early part of the distribution for ages under half a year is generally overestimated, while the intermediary part, around one year, is underestimated. The shapes of the gamma and exponential distributions are too simple to represent the early part of the “true” distribution. In the gamma distribution, the large dispersion of distributions at early ages also shows that the early

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Simulated age distribution using the inverse particle tracking approach for the Plœmeur aquifer and 60 realizations from the inferred plausible age distributions based on five tracers including CFC-11, CFC-12, CFC-113, SF₆ and ⁸⁵Kr. The top row represents the results when a gamma RTD is assumed and the bottom row shows the results when the RTD is assumed to be exponential. The number of years used for inference of the RTD are one year in panels a and d (left), four in panels b and e (middle) and eight in panels c and f (right).</figcaption></figure>
part of the distribution is not very sensitive to the tested tracers, and that the information content of the environmental tracers is limited at very early ages. Differences between the "true" distribution and the analytical approximation come from the topographi-

cal and geological details close to the well. First, the well pumps water at a depth with a minimum time necessary for the water to reach the well from the surface. No such assumptions are taken into account in the gamma and exponential distributions. Second,

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Simulated age distribution using the inverse particle tracking approach for the Plœmeur aquifer and 60 realizations from the inferred plausible age distributions based on four tracers including CFC-11, CFC-12, CFC-113, and SF6. The top row represent the results when a gamma RTD is assumed and the bottom row shows the results when the RTD is assumed to be exponential. The number of years used for inference of the RTD are one year in panels a and d (left), four in panels b and e (middle) and eight in panels c and f (right).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** "Observed" tracers concentrations and 95% credible intervals obtained from the posterior distribution of gamma RTDs based on 8 years (top row), four years (middle row) and one year (bottom) of tracer data.</figcaption></figure>
the shape of the "true" distribution between zero and three years is sensitive to the details of the topography. More details about the effects of near well topography on the distribution of the very young (<3 years) of groundwater is provided in Leray et al., 2012. These details produce the irregular shape of the distribution, and they cannot be homogenized into a smooth curve for the later time that concerns a much larger part of the domain.

In fact, the part of the RTD greater than three years is captured much better using both the gamma and exponential forms. If the estimated distributions cannot be used for the very early ages, it does not preclude their applicability at later ages. In both cases of four and five tracers, the uncertainty associated with the earlier parts of the distributions increases as the number of years of data used is increased from one to four, and then it slightly decreases (Fig. 6). This may appear counter-intuitive, as more data would seem to reduce the uncertainty. However, the reason for the increase in uncertainty when the number of years of data increases from one to four can be explained better by looking at Fig. 7, which shows the "observed" tracer concentrations and the 95% credible intervals of posterior concentrations.

### 3.2.2. Information content of environmental tracer concentrations

The credible intervals of tracer concentrations when using the four-year time series are about twice as large as those obtained when using a single year (Fig. 7, middle line compared to bottom

line). When four years of data are used, the model is more constrained and the observation and model structural errors are reflected more appropriately in the uncertainty associated with the posterior parameters. The multiplication of the data introduces some possible tradeoffs, with some models possibly being better than others at some ages and worse at other ages. It does not mean that the models are more uncertain, but that the uncertainty of the model is better constrained.

Another observation from Fig. 7 is that the temporal trends in the concentrations of some of the tracers provide more information about the RTD than other tracers. For example, the temporal pattern of the concentration of $^{85}$Kr is completely disrupted by the addition of the 5% random noise, while the temporal trend is preserved for most of the other tracers. This is due to the fact that the historical atmospheric $^{85}$Kr concentrations have been relatively flat in recent years (Fig. 1). Conversely, in the case of SF$_{6}$, and to a lesser degree, CFCs, evaluating the 95% credible intervals of tracer concentrations for different years shows that the temporal variations in input concentrations have a statistically significant impact on the observed concentrations at the well. The smaller uncertainty associated with the posterior SF$_{6}$ compared to the CFCs can be attributed to the fact that observed SF$_{6}$ concentrations varied monotonically by year for the past eight years, while the variation in CFCs was not strictly monotonic.

<table><caption>**Table 4**
Measures of goodness of models including Normalized Bayes Factor and DIC and measures based on the ability of the method to predict cumulative residence time at 10 years for exponential and gamma RTDs as a result of using four and five tracers and one, four and eight years of tracer data.</caption><thead><tr><th rowspan="2"></th><th colspan="3">BF</th><th colspan="3">DIC</th><th colspan="3">100 × δ<sub>10</sub></th><th colspan="3">1000 × b<sub>10</sub></th></tr><tr><th>1 yr (%)</th><th>4 yrs (%)</th><th>8 yrs (%)</th><th>1 yr</th><th>4 yrs</th><th>8 yrs</th><th>1 yr</th><th>4 yrs</th><th>8 yrs</th><th>1 yr</th><th>4 yrs</th><th>8 yrs</th></tr></thead><tbody><tr><td>CFCs, SF<sub>6</sub></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>exp</td><td>29</td><td>48</td><td>7</td><td>38.60</td><td>66.10</td><td>131.51</td><td>5.72</td><td>3.16</td><td>2.09</td><td>56.03</td><td>6.86</td><td>12.99</td></tr><tr><td>gamma</td><td>71</td><td>52</td><td>93</td><td>39.59</td><td>65.50</td><td>136.22</td><td>5.16</td><td>3.97</td><td>1.50</td><td>49.78</td><td>19.99</td><td>1.88</td></tr><tr><td>CFCs, SF<sub>6</sub> <sup>85</sup>Kr</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>exp</td><td>48</td><td>61</td><td>48</td><td>36.20</td><td>99.28</td><td>177.35</td><td>7.25</td><td>1.70</td><td>1.15</td><td>68.57</td><td>7.62</td><td>0.25</td></tr><tr><td>gamma</td><td>52</td><td>39</td><td>52</td><td>34.56</td><td>97.84</td><td>177.11</td><td>8.66</td><td>3.57</td><td>1.10</td><td>55.85</td><td>4.52</td><td>0.66</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** True fraction of water with age greater than 10 years ($\bar{\epsilon}_{10}$) (vertical dashed line) and the inferred posterior distribution of $\bar{\epsilon}_{10}$ obtained using five tracers: CFC-11, CFC-12, CFC-113, SF$_{6}$, and $^{85}$Kr when one (left), four (middle), and eight (right) years of tracer data are used assuming gamma (top) and exponential (bottom) RTDs.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** True fraction of water with age greater than 10 years ($\bar{\xi}_{10}$) and the inferred posterior distribution of $\bar{\xi}_{10}$ obtained using four tracers: CFC-11, CFC-12, CFC-113, and SF$_6$ when one (left), four (middle), and eight (right) years of tracer data are used assuming gamma (top) and exponential (bottom) RTDs.</figcaption></figure>

Various measures of adequacy of the models are summarized in Table 4. The method does not seem to provide a definitive conclusion about the forms of RTDs. When four tracers are used, the method tends to assign a larger likelihood (although not substantially) to the gamma RTD (except in the case of four years of data), and when five tracers are used, the exponential and gamma models are given roughly equal odds. This can be due to the fact that neither the gamma nor the exponential model perfectly represents the “true” RTD, and neither model is able to fit all parts of it perfectly, especially the early part, as seen before. The DIC values also show that as more data are used, a smaller discrimination between the models is observed. It is interesting to point out that when four tracers are used, the exponential model is given preference in terms of DIC in all cases, while when five tracers are used, the gamma model is always better. Because DIC explicitly accounts for model complexity in determining the goodness of fit, this can be interpreted as meaning that four tracers are not adequate for characterizing the more complex gamma distribution, but when five tracers are used, the gamma distribution can be better determined.

It is also important to note that the $^{85}$Kr tracer that is added in the case of five tracers contains more information about the younger ages, due to its relatively short half-life.

### 3.2.3. Prediction capacity of analytical age distributions

The uncertainty measure, $\delta_{10}$, and the bias, $b_{10}$, decrease sharply as the number of years of data is increased (Table 4). Uncertainty $\delta_{10}$ is reduced by a factor of 3–7 when using eight years of data rather than one year, while $b_{10}$ is reduced by a factor of 3 to almost 100, depending on the number of tracers and distribution type (exponential or gamma). The reduction is systematic, except for the case of the exponential model from four to eight years of data. The distribution of the fraction of water age greater than ten years, $\bar{\xi}_{10}$, gives a more accurate illustration of the interest of the age data time series (Figs. 8 and 9). Using only one year of data, the distribution of $\bar{\xi}_{10}$ is both biased and uncertain, to the point that the “true” value is given a very low probability. With four years of data, the bias is significantly reduced. With eight years, the uncertainty

is further reduced. Several years of data reduce both the bias and the uncertainty. The good predictions on $\bar{\xi}_{10}$ are also confirmed by the small values of bias $b_{10}$ and uncertainty $\delta_{10}$ for the eight-year case, except for the exponential case without $^{85}$Kr (Table 4, first line). This highlights the effectiveness of the age data time series for the predictions of the relatively “older” water age fraction.

The very small influence of the distribution nature on the predictions is also noteworthy. Whether the exponential or the gamma distribution is used barely affects the distributions of $\bar{\xi}_{10}$, as predictions are much more sensitive to the fraction of ages greater than 10 years than to the distribution shape. This might be explained by the strong consistency of both distribution shapes at periods longer than ten years (Fig. 5). Depending on the targeted application, the shape of the distribution might not be an essential factor for the predictions.

## 7. Conclusion

In this paper, the value of age data time series for the inference of groundwater age distribution was studied. For this purpose, a hypothetical gamma and a simulated age distribution obtained using an aquifer model of the Plœmeur site were used to generate synthetic tracer concentrations for five tracers—3 CFCs, SF$_6$, and $^{85}$Kr—at the end of eight consecutive years ending in 2010. Bayesian inference by means of MCMC was used to infer the parameters of presumed mathematical age distribution forms, including exponential and gamma, when variable numbers of years of data were used. In addition, the effect of excluding $^{85}$Kr in the inference uncertainty was evaluated by repeating the analysis without $^{85}$Kr. Using DIC, the main conclusions that can be drawn from this study are:

- In the case where the tracer data were obtained using the hypothetical age distribution, using a larger number of data will result in a much more pronounced preference given to the right RTD form in terms of DIC.
- The two mathematical RTD forms considered here (exponential and gamma) are too simple to be able to capture the early part of the true RTDs. The fact that the earlier ages cannot be captured well is due to their dependency on local, non-homogeneous aquifer characteristics in the proximity of the pumping well.

- The actual bias and uncertainty of the inferred age distributions based on comparing the cumulative ten-year ages ($\xi_{10}$) decrease as the number of years of tracer data increases. Both bias and uncertainty evaluated based on 10 years cumulative age decrease when a larger number of years of data is used in the case when five tracers were used. When using four tracers one year of data resulted in a smaller spread for $\xi_{10}$ with a substantially larger bias compared to the case when four years of data was used, which might imply high confidence about the outcome in the absence of knowledge about the true RTDs. This means that more data does not necessarily lead to a reduction in the uncertainty as expressed by the spread of the parameters, but it leads to better assessment of the uncertainty. This is especially true when the true RTD cannot be fully captured by the assumed mathematical forms.

- In the case of the simulated RTD, the method cannot definitively pick one of the RTD forms over the other one. However, based on the DIC measure, when the additional tracer, $^{85}$Kr, is included in the analysis, the more complex form (gamma) is given higher odds, while the exponential form is given preference when only four tracers are used.

- Regardless of the distribution shape, the reliability of the inferred RTD improves with regard to predicting the ten-year cumulative age, and the uncertainties and bias associated with gamma and exponential models are close to each other.

- Shape-free RTDs can liberate us from presuming mathematical forms that restrict the age distribution to a limited range of shapes. However, due to larger degrees of freedom, they require a greater number of tracers to be constrained. Transient tracer information is a promising way to provide the amount of data needed for such an approach.

## Acknowledgments

Research and monitoring on the site of Płœmeur is funded by the Environmental research Observatory H+ (Network of hydrogeological sites) and by the European Interreg IV project Climawat. J.-R. de Dreuzy acknowledges the European Union for its additional funding through the IEF Marie-Curie fellowship (PIEF-GA-2009-251710).

## References

Amin, I.E., Campana, M.E., 1995. A general lumped parameter model for the interpretation of tracer data and transit time calculation in hydrologic systems. Journal of Hydrology 179, 1–21.

Bauer, S., Fulda, C., Schäfer, W., 2001. A multi-tracer study in a shallow aquifer using age dating tracers $^{3}$H, $^{85}$Kr, CFC-113 and SF$_{6}$ – indication of retarded transport of CFC-113. Journal of Hydrology 248, 14–34.

Bear, J., 1991. Modelling and Applications of Transport Phenomena in Porous Media. Kluwer Academic Publishers, Dordrecht.

Busenberg, E., Plummer, L.N., 2008. Dating groundwater with trifluoromethyl sulfur pentafluoride (SF$_{5}$CF$_{3}$), sulfur hexafluoride (SF$_{6}$), CF$_{3}$Cl (CFC-13) and CF$_{2}$Cl$_{2}$ (CFC-12). Water Resources Research 44, W02431.

Carlin, B.P., Chib, S., 1995. Bayesian model choice via Markov-chain Monte-Carlo methods. Journal of the Royal Statistical Society Series B – Methodological 57 (3), 473–484.

Carn, A., 1990. Mise en valeur des ressources en eau souterraine du socle breton. Département du Morbihan (35) – R31724 BRE 4S/90Rep., BRGM, Rennes.

Cook, P.G., Love, A.J., Robinson, N.I., Simmons, C.T., 2005. Groundwater ages in fractured rock aquifers. Journal of Hydrology 308 (1–4), 284–301.

de Marsily, G., 1986. Quantitative Hydrogeology: Groundwater Hydrology for Engineers, first ed. Academic Press, New-York, 440 pp.

Fienen, M.N., Luo, J., Kitanidis, P.K., 2006. A Bayesian geostatistical transfer function approach to tracer test analysis. Water Resources Research 42, W07426. http://dx.doi.org/10.1029/2005WR004576.

Gamerman, D., Hedibert, F.L., 2006. Markov Chain Monte Carlo – Stochastic Simulation for Bayesian Inference, second ed. Chapman & Hall/CRC.

Geweke, J., 1992. Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. In: Bernardo, J.M., Berger, J.O., Dawid, A.P., Smith, A.F.M. (Eds.), Bayesian Statistics. Oxford University Press.

Geweke, J., Tanizaki, H., 2001. Bayesian estimation of state-space models using the Metropolis-Hastings algorithm within Gibbs sampling. Computational Statistics & Data Analysis 37 (2), 151–170.

Glynn, P.D., Plummer, L.N., 2005. Geochemsity and the understanding of groundwater systems. Hydrogeology Journal 13, 263–287.

Haitjema, H.M., 1995. On the residence time distribution in idealized groundwathesheds. Journal of Hydrology 172 (1–4), 127–146.

IAEA, 2006. Use of Chlorofluorocarbons in Hydrology: A Guidebook. International Atomic Energy Agency, Vienna, 277 pp.

Ivey, S.S., Gentry, R.W., Larsen, D., Anderson, J., 2008. Inverse application of age-distribution modeling using environmental tracers $^{3}$H/$^{3}$He. Journal of Hydrologic Engineering 13 (11), 1002–1010.

Jeffreys, H., 1935. Some tests of significance, treated by the theory of probability. Proceedings of the Cambridge Philosophical Society 31, 203–222.

Kaipio, J., Somersale, E., 2004. Statistical and Computational Inverse Problems (Applied Mathematical Sciences), vol. 60. Springer.

Kass, R.E., Raftery, A.E., 1995. Bayes factors. Journal of the American Statistical Association 90 (430), 773–795.

Kreft, A., Zuber, A., 1978. On the physical meaning of the dispersion equation and its solution for different initial and boundary conditions. Chemical Engineering Science 33 (11), 1471–1480.

LaBolle, E.M., Fogg, G.E., 2001. Role of molecular diffusion in contaminant migration and recovery in an alluvial aquifer system. Transport in Porous Media 42 (1–2), 155–179.

Le Borgne, T., Bour, O., Paillet, J.-L., Caudal, J.-P., 2006. Assessment of preferential flow path connectivity, and hydraulic properties at single-borehole and cross-borehole scales in a fractured aquifer. Journal of Hydrology 328 (1–2), 347–359.

Leray, S., de Dreuzy, J.R., Bour, O., Labasque, T., Aquilina, L., 2012. Contribution of age data to the characterization of complex aquifers. Journal of Hydrology 464, 54–68.

Leray, S., de Dreuzy, J.R., Aquilina, L., Vergnaud-Ayraud, V., Labasque, T., Bour, O., Le Borgne, T., 2013. Temporal evolution of age data under transient pumping conditions. Journal of Hydrology, submitted for publication.

Long, A.J., Putnam, L.D., 2009. Age-distribution estimation for Karst groundwater: issues of parameterization and complexity in inverse modeling by convolution. Journal of Hydrology 376, 579–588.

Luo, J., Cirpka, O.A., Fienen, M.N., Wu, W.-m., Mehlhorn, T.L., Carley, J., Jardine, P.M., Criddle, C.S., Kitanidis, P.K., 2006. A parametric transfer function methodology for analyzing reactive transport in nonuniform flow. Journal of Contaminant Hydrology 83 (1–2), 27–41.

Maloszewski, P., Zuber, A., 1982. Determining the turnover time of groundwater systems with the aid of environmental tracers, 1. Models and their applicability. Journal of Hydrology 57 (3–4), 207–231.

Manning, A.H., Clark, J.F., Diaz, S.H., Rademacher, L.K., Earman, S., Plummer, L.N., 2012. Evolution of groundwater age in a mountain watershed over a period of thirteen years. Journal of Hydrology 460, 13–28.

Massoudieh, A., Sharifi, S., Solomon, D.K., 2012. Bayesian evaluation of groundwater age distribution using a radioactive tracers and anthropogenic chemicals. Water Resources Research 48, W09529. http://dx.doi.org/10.1029/2012WR011815.

McGuire, K.J., McDonnell, J.J., 2006. A review and evaluation of catchment transit time modeling. Journal of Hydrology 330 (3–4), 543–563.

Neupauer, R.M., Wilson, J.L., 1999. Adjoint method for obtaining backward-in-time location and travel time probabilities of a conservative groundwater contaminant. Water Resources Research 35 (11), 3389–3398.

Neupauer, R.M., Wilson, J.L., 2001. Adjoint-derived location and travel time probabilities for a multidimensional groundwater system. Water Resources Research 37 (6), 1657–1668.

Newman, B.D., Osenbrück, K., Aeschbach-Hertig, W., Solomon, D.K., Cook, P.G., Rozanski, K., Kipfer, R., 2010. Dating of ‘young’ groundwaters using environmental tracers: advantages, applications, and research needs. Isotopes in Environmental and Health Studies 46 (3), 259–278.

Osenbrück, K., Fiedler, S., Knöller, K., Weise, S.M., Sültenfuß, J., Oster, H., Strauch, G., 2006. Timescales and development of groundwater pollution by nitrate in drinking water wells of Jahna-Aue, Saxonia, Germany. Water Resources Research 42, W12416.

Portnaguine, O., Solomon, D.K., 1998. Parameter estimation using groundwater age and head data, Cape Cod, Massachusetts. Water Resources Research 34 (4), 637–645.

Reilly, T.E., Plummer, L.N., Phillips, P.J., Busenberg, E., 1994. The use of simulation and multiple environmental tracers to quantify groundwater flow in a shallow aquifer. Water Resources Research 30 (2), 421–433.

Robertson, W.D., Cherry, J.A., 1989. Tritium as an indicator of recharge and dispersion in a groundwater system in Central Ontario. Water Resources Research 25 (6), 1097–1109.

Ruelle, S., Moreau, F., Bour, O., Gapais, D., Martelet, G., 2010. Impact of gently dipping discontinuities on basement aquifer recharge: an example from Płœmeur (Brittany, France). Journal of Applied Geophysics 70 (2), 161–168.
Solomon, D.K., Poreda, R.J., Schiff, S.L., Cherry, J.A., 1992. Tritium and helium 3 as groundwater age tracers in the Borden aquifer. Water Resources Research 28 (3), 741–755.

Spiegelhalter, D.J., Best, N.G., Carlin, B.R., van der Linde, A., 2002. Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society Series B – Statistical Methodology 64, 583–616.

Sturchio, N.C. et al., 2004. One million year old groundwater in the Sahara revealed by krypton-81 and chlorine-36. Geophysical Research Letters 31, L05503. http://dx.doi.org/10.1029/2003GL019234.

Sültenfuß, J., Purtschert, R., Führböter, J.F., 2011. Age structure and recharge conditions of a coastal aquifer (northern Germany) investigated with <sup>39</sup>Ar, <sup>14</sup>C, <sup>3</sup>H, He isotopes and Ne. Hydrogeology Journal 19, 221–236.

Touchard, F., 1999. Caractérisation hydrogéologique d'un aquifère en socle fracturé: Site de Plœmeur (Morbihan). PhD Thesis, University of Rennes 1, France, 343 pp.

Troldborg, L., Refsgaard, J.C., Jensen, K.H., Engesgaard, P., 2007. The importance of alternative conceptual models for simulation of concentrations in a multi-aquifer system. Hydrogeology Journal 15 (5), 843–860.

Troldborg, L., Jensen, K.H., Engesgaard, P., Refsgaard, J.C., Hinsby, K., 2008. Using environmental tracers in modeling flow in a complex shallow aquifer system. Journal of Hydrologic Engineering 13 (11), 1037–1048.

Zhang, Y., 2004. Numerical simulations of dating young groundwater with multiple atmospheric tracers: CFC-11, CFC-12, SF<sub>6</sub>, <sup>3</sup>H/<sup>3</sup>He and <sup>85</sup>Kr. In: Paper presented at 15th International Conference on Computational Methods in Water Resources. Elsevier Science BV, Chapel Hill, NC (June 13–17, 2004).