MODELING OF GROUNDWATER AGE ANDRESIDENCE-TIME DISTRIBUTIONS

In this chapter, we try to mathematically describe how groundwater ages obtained by various methods—as described in the previous chapters—can be different from or similar to the real ages of each groundwater molecule/par-ticle we sample in any dating study. Such a chapter is unprecedented in the groundwater textbooks, especially those that target Earth science specialists. Therefore, in contrast to all other chapters, it is necessary to start with a broad overview of the subject and to proceed step by step to help the average reader have a general impression of the content of the chapter.
The migration of groundwater particles through an aquifer is governed by the physical transport processes of advection, dispersion, and diffusion, which are described by the well-known flow and transport equations (see, e.g., Freeze and Cherry, 1979; Wang and Anderson, 1982; Huyakorn and Pinder, 1983; de Marsily, 1986; Fetter, 1999). As an intrinsic property of these moving particles, their respective age can therefore be assimilated to an extensive quantity evolving in space depending on the structure of the flow field, itself dictated by aquifer heterogeneity and recharge and discharge conditions. Hence, in any aquifer, groundwater particles move at various speeds along flow paths of various shapes and lengths. The morphology of the flow paths characterizes the geometric nature of a given flow system and determines the time water and rock are in contact. In regional systems, ages or residence times are often found between very short and very large values. The major consequence of this is that a groundwater sample taken in the aquifer, or at any of its outlets,includes a range, or a distribution of particle ages. Hence, full groundwater age information at a given location must be represented by a frequency distribu-tion, in the statistical sense, and not by a single value. This distribution may possibly exhibit a complex shape and cover a wide variety of age values depending on the aquifer nature and its hydrodispersive components. For example, multimodal, asymmetric, and skewed distributions are expected in cases of multilayer aquifers, fracture porosity, and fault systems. Figure 7.1 illustrates the residence-time variations in relation to the many possible flow paths leading to the discharge zone of a regional multilayer system. Depend-ing on where infiltration takes place, these flow paths may reach different depths, and imply extremely different travel times to the outlet, affecting water chemistry and mineralization accordingly.
Groundwater age-dating by various tracers described in Chapters 4–6 gen-erally gives us an average value for a given sample. This value as described in Chapter 3 is certainly very instructive and useful but gives no information at all on the actual age distribution and its possible complex features, which remain unknown. In fact, many age distributions can have the same average value (statistical mean, or first moment), and this average is sometimes far from the modes of the unknown distribution, namely the most frequent age occurrences in the sample. Using the mean ages to evaluate aquifer charac-teristics such as recharge rates and flow velocities should therefore be made very carefully, and in principle only for relatively simple hydrogeological con-figurations for which standard age distributions can be enforced. For other, more complex flow systems, the mean age interpretation may lead to erro-neous conclusions with sometime severe consequences when water quality and water protection issues are addressed.
To date, the only ways to determine groundwater age and residence-time distributions are those of mathematical modeling based on the solution of age transport equations. To perform accurately, these approaches require full hydrodynamical description of the aquifer under investigation and the knowl-edge of its flow field. In general, the flow field is derived from the solution of the steady-state groundwater flow equations for calibrated hydrodynamical data sets. The calibrated flow field is then used as input, as well as other hydrodispersive parameters, to calculate the age transport processes and the resulting age distributions.
The following sections provide an overview on groundwater age modeling techniques from the 1960s up to the present state-of-the-art, and some mathematical basics on groundwater age transport, groundwater age, and residence-time distributions, as well as illustrative typical examples.
7.1 OVERVIEW AND STATE-OF-THE-ART
To quantify the distribution of ages in aquifers, several types of mathematical models have been developed during the past decades, and research in this field is progressing rapidly. The classical age modeling practice consists of using simple analytical models, commonly called lumped-parameter models, to interpret environmental tracer data. This type of approach has been exten-sively used and is still very popular (e.g., Maloszewski and Zuber, 1982, 1993; Campana and Simpson, 1984; Balderer, 1986; Richter et al., 1993; Amin and Campana, 1996). With lumped-parameter models, simple specific age distri-butions describing piston flow, exponential mixing, combined piston-exponential mixing, or dispersive mixing are assumed beforehand, sometimes very arbitrarily. The inverse problem is then solved by fitting the model response to tracer measurements by adjustment of some parameters (see Figure 7.2). This procedure calls for significant simplifications, which are often not justified, such as neglecting the reservoir structure, as well as the spatial variability of infiltration rates (Campana, 1987) and aquifer flow and transport parameters. Amin and Campana (1996) proposed to model the mixing process effects by means of a three-parameter gamma function ranging between no mixing (piston-flow model) and perfect mixing (exponential model).
Lumped-parameter models are often used to make inferences on ground-water flow velocities, aquifer turnover time, and recharge rates. However, sound verifications of the reliability of such models can hardly be found (e.g.,see Haitjema, 1995; Luther and Haitjema, 1998; Etcheverry and Perrochet, 2000). Haitjema (1995) provided some guidelines for the use of lumped-parameter models and derived an analytical solution of the exponential model for semiconfined stratified aquifers. Luther and Haitjema (1998) claimed that the conditions for the validity of the exponential model (horizontal flow and homogeneous aquifer with respect to porosity φ, recharge rate i, and saturated thickness e) can be extended to configurations where the ratio φe/i (namely the only free parameter in the exponential model) remains constant through-out the aquifer. This ratio is the system turnover time, or average residence time, which means that a water sample taken anywhere in the reservoir must give a mean age corresponding to the aquifer turnover time. Such a configu-ration and these uniform mean age conditions may hardly be found in nature, making the assumptions of Luther and Haitjema (1998) rather limiting. As demonstrated by Etcheverry (2001) a simple linear variation of the thickness e significantly influences the shape of the theoretical exponential age distri-bution. Classically, and for simplicity, age transport has also often been assumed driven by the only process of advection. In doing so, average ages can indeed be calculated easily, either analytically if the flow configuration is simple, or numerically by postprocessing the results of a groundwater flow model with particle-tracking techniques. Given a known velocity field, parti-cle-tracking algorithms basically track down a certain number of released particles along path lines and during selected times. This can be done either forward by following the flow lines downstream or backward by following the flow lines upstream. However, the resulting advective ages ignore the effect of dilution, dispersion, and mixing on age transport (Cordes and Kinzelbach,1992) and often reveal to be ill posed in complex heterogeneous systems (Varni and Carrera, 1998), for which the three-dimensional implementation is subject to severe technical problems. Moreover, particle-tracking techniques do not allow the calculation of residence-time distributions, since the moving groundwater volumes are not associated to the simulated ages. Nevertheless, the particle tracking method is very popular and may provide some help to validate dating methods, mainly isotopic methods, which consider the disinte-gration of radioactive isotopes as the age index (Smith et al., 1976).
An example of advective age simulation by particle tracking is shown in Figure 7.3(a), representing a plan view of a horizontal aquifer with a pumping well. Backward particle tracking is applied from the well, and isochrone markers at specific times (black dots) give an idea of the life expectancy of water particles in the well capture zone (i.e., the time left before absorption of these particles by the well). The advective capture zone of the well can therefore be delineated by the envelope of all isochrone markers. Following this approach, water particles are either in or out of the capture zone. However, it is clear that water particles also move laterally under dispersion and diffusion processes and, thus, depending on their distance to the well, may eventually exit or enter the capture zone. When comparisons are made between modeled and measured ages, the importance of accounting for age diffusion and dispersion has been pointed out by many authors (e.g., Sudicky and Frind, 1981; Maloszewski and Zuber, 1982, 1991).
More elaborated quantitative approaches consider groundwater age as an extensive property, where transport through an aquifer can be described by volume-averaged field equations of the advective-dispersive type (e.g., Spalding, 1958; Harvey and Gorelick, 1995; Goode, 1996; Ginn, 1999). Com-paring various modeling results with the experimental data of a complex system, Castro and Goblet (2005) found that the mean age field resulting from advective-dispersive age mass transport (Goode, 1996) was the most consis-tent. Varni and Carrera (1998) derived a set of recursive, temporal moment equations that were also compared to radiometric age measurements. Accord-ing to Harvey and Gorelick (1995), the first four moments (characterizing the mean, the variance, the skewness, and the kurtosis of a breakthrough curve, respectively) may provide sufficient information to evaluate the entire age dis-tribution at a given location. However, this finding only holds in cases where the principle of maximum entropy applies. Since many natural systems are likely to reveal multimodal age distributions of complex shapes, within the reservoir and at the discharge zones, this approach to construct entire age dis-tributions (a priori unknown) remain somewhat uncertain, or nonunique. For instance, it is fairly obvious that calculations of mean age, or of age variance, at a given location (with first and second temporal moment equations) do not represent the full age picture since many different age distributions may exhibit the same mean and variance. Examples of mean age and mean life expectancy simulations including dispersion and diffusion effects are given in Figure 7.3(b) and (c). The mean age and mean life expectancy fields are obtained by two solutions of the advective-dispersive mean age equation (Goode, 1996), with forward and backward flow field, respectively. The mean life expectancy field in Figure 7.3(c), where the isolines indicate the time required by groundwater particles to reach one of the two outlets, can be com-pared to the advective backward particle-tracking solution in Figure 7.3(a), and the effects of dispersion on the size of the well capture zone can be assessed. However, none of these advective or advective-dispersive calcula-tions provides any insight on the entire age or life expectancy distributions at a given location. For example at the extraction well, where many groundwater fractions with different ages mix, the knowledge of the entire age distribution shown in Figure 7.3(d) is of particular interest. However, approaches based on temporal moment analysis are not suitable for this purpose, and improved models must be enforced to achieve this task in a deterministic manner for any arbitrary hydrogeological configuration (Cornaton, 2004; Cornaton and Perrochet, 2005a, b).
Travel time probabilities have also been a subject of high interest in many studies characterizing solute transport in subsurface hydrology (e.g., Dagan, 1982, 1987, 1989). The travel time probability is commonly defined as the response function to an instantaneous unit mass flux impulse (Danckwerts, 1953). In their transfer function approach of contaminant transport through unsaturated soil units, Jury and Roth (1990) model tracer breakthrough curves with one-dimensional travel time probability functions. Shapiro and Cvetkovic (1988) and Dagan and Nguyen (1989) derived the forward travel time proba-bility for a mass of solute by using the Lagrangian concept of particle dis-placement in porous media. Using a stochastic, longitudinal exponential model for the travel time density function, Woodbury and Rubin (2000) combined this approach with full-Bayesian hydrogeological parameter inference by inverting the travel time moments of solute breakthroughs in heterogeneous aquifers.
The derivation of forward and backward models for location and travel time probabilities has become a classical mathematical approach for contam-inant transport characterization and prediction (e.g., van Herwaarden, 1994; van Kooten, 1995; Neupauer and Wilson 1999, 2001). The spreading of a contaminant mass is analyzed by following the random motion of solute par-ticles, and, to do so, the advection-dispersion equation is assimilated to the Fokker–Planck, or forward Kolmogorov equation. The expected resident con-centration of a conservative tracer is taken as the probability density function for the location of a particle, at any time after entry into the system. LaBolle et al. (1998) recall that the standard diffusion theory, which relates the dynam-ics of a diffusion process to Kolmogorov equations, may apply at local scale to advection-dispersion equations if the porosity and the dispersion tensor are functions varying smoothly in space. To solve these equations, random-walk procedures are often chosen as alternatives to more standard discretization methods (see, for example, Kinzelbach, 1992; LaBolle et al., 1998; Weissmann et al., 2002).
The derivation of adjoint state equations to define backward location and travel time probabilities was recently revisited by Neupauer and Wilson (1999, 2001). They derived specific forward and backward probabilities to predict contaminant sources from concentrations measured at monitoring and extrac-tion wells. Compared to the forward modeling approach used to predict the future evolution of contamination when the source is known, the backward modeling approach reveals to be useful when the contamination has been detected and the possible sources have to be identified.
Following a different line of research, Eriksson (1961, 1971), followed by Bolin and Rhode (1973), introduced the reservoir theory in the field of envi-ronmental problems. Originating from chemical engineering, this theory con-siders a reservoir globally and is particularly useful to address residence-time issues. In essence, when applied to conservative particles flowing through a bounded system, the theory links the distribution of particle ages in the inte-rior of the system to the residence-time distribution of the particle population. However, in spite of its fundamental features and capabilities, this approach has not received much attention in hydrogeology until recently (Etcheverry and Perrochet, 2000; Etcheverry, 2001). These authors proposed a method to directly calculate the groundwater residence-time distributions at the outlet of an aquifer of arbitrary heterogeneity by combining the age mass transport equation of Goode (1996) to the reservoir theory. One significant advantage of the method is that full temporal information is recovered at the outlet based on the solution of only two steady-state equations. Moreover, based on the internal organization of ages in the entire reservoir, the method minimizes the loss of age information due to mixing in the vicinity of local outlets where flow paths converge. The potential of this approach was investigated further by Cornaton (2004) and Cornaton and Perrochet (2005a, b) who generalized the reservoir theory to hydrodispersive systems with multiple inlets and outlets, and combined it to forward and backward Laplace-transformed equations for location and travel time probabilities. Achieving a particularly high computa-tional efficiency, these recent works open a range of new research and appli-cation perspectives, not only in the field of age-dating, but also in associated contaminant hydrogeology issues involving large time scales. Facing the math-ematical modelers’ community, sustainable aquifer management with respect to climatic changes, accurate delineation of well-head protection zones, risk and safety assessments of nuclear and chemical waste repositories, among others, are now all 21st-century crucial issues in which groundwater age, life expectancy, and residence-time distributions play a major role in the long term.
7.2 BASICS IN GROUNDWATER AGE TRANSPORT
This section introduces the basics for modeling groundwater age and residence-time distributions by combining the reservoir theory to determinis-tic mathematical flow and transport models. This is made using a formalism kept as simple as possible, so that the reader can understand the principles of the method with standard mathematical knowledge and calculus techniques. However, notions on the modeling of groundwater flow and transport processes, whose presentation is not the scope of this chapter, are recom-mended and can be found in influential textbooks such as those by Freeze and Cherry (1979), Wang and Anderson (1982), Huyakorn and Pinder (1983), de Marsily (1986), and Fetter (1999).
7.2.1 The Reservoir Theory
The reservoir theory was introduced in environmental sciences by Eriksson (1961, 1971) to address atmospheric and surface hydrology issues. In the groundwater context, this theory can be developed from the divergence theorem applied to particle age-fluxes through a bounded domain (Etcheverry, 2001; Cornaton, 2004). This global approach enables us to quan-tify the intrinsic relation between age and residence-time distributions in any hydrodispersive and mixing flow regimes, and for multidimensional configu-rations of arbitrary heterogeneity. However, the essential features of the reser-voir theory are more easily derived and interpreted for averaged steady-state flow systems, as generally assumed in groundwater modeling practices. In addi-tion, the notion of local average ages is adopted in the following sections. This last restriction is enforced for the sake of mathematical simplicity without altering the generality and the usefulness of the approach to handle full, locally dispersed age and residence-time distributions.
Relation Between Internal Ages and Residence Times The simple system in Figure 7.4 represents a homogeneous aquifer of total pore volume Vo and total discharge rate Qo. The various flow paths taken by water infiltrating over the inlet imply a range of residence times over the outlet and a range of ages in the aquifer. Given one particular, average age contour τ in the system, ground-water volumes can be arranged in a cumulative manner, so that the function V(τ) represents the resident volumes where ages are smaller than τ. Similarly, over the outlet, the discharge can be arranged in a cumulative manner, so that the function Q(τ) represents the volumetric rate of water leaving the system after a residence time smaller than τ. Considering now mass balance, the total recharge entering the system during the interval τ, namely the infiltrated volume Qoτ, must correspond to the sum of the groundwater volume V(τ) and the volume exfiltrated, during the interval, with a residence time smaller than τ. Hence, the relation between Q(τ) and V(τ) can be written as