Consistent Boundary Conditions for Age Calculations

: Age can be evaluated at any time and position to understand transport processes taking place in the aquatic environment, including for reactive tracers. In the framework of the Constituent-oriented Age and Residence time Theory (CART), the age of a constituent or an aggregate of constituents, including the water itself, is usually deﬁned as the time elapsed since leaving the boundary where the age is set or reset to zero. The age is evaluated as the ratio of the age concentration to the concentration, which are the solution of partial differential equations. The boundary conditions for the concentration and age concentration cannot be prescribed independently of each other. Instead, they must be derived from boundary conditions designed beforehand for the age distribution function (the histogram of the ages, the age theory core variable), even when this variable is not calculated explicitly. Consistent boundary conditions are established for insulating, departure and arrival boundaries. Gas exchanges through the water–air interface are also considered. Age ﬁelds ensuing from consistent boundary conditions and, occasionally, non-consistent ones are discussed, suggesting that the methodology advocated herein can be utilized by most age calculations, be they used for diagnosing the results of idealised models or realistic ones.


Introduction
Today's numerical models of geophysical and environmental fluid flows and the related (reactive) transport processes produce huge output files. Making sense of all these real numbers (that is, identifying key processes and establishing causal relationships between them) is no trivial task. Analysing primitive variables (e.g., velocity, pressure, temperature, concentrations) is not always conducive to the most fruitful interpretations. Evaluating auxiliary variables introduced for diagnostic purposes is an option worth considering. In this respect, diagnostic timescales (e.g., age, transit time, residence or exposure time) have been of use in the modelling of the atmosphere [1][2][3][4], various types of water bodies , and sediment transport [39,[50][51][52]. These diagnostic timescales paint a simplified picture of the impact that the phenomena under study have on the largest time and space scales of (reactive) transport.
The residence or exposure time [16,53] and the age are of fundamentally different natures, as may be seen, for instance, in Figure 13 of [54]. The former type of timescale looks into the future, whilst the assumed to be positive definite, that is, τ ∈ [0, ∞[. Function c(t, x, τ) is defined as follows: in the abovementioned control volume, at time t, the mass of the particles of the constituent under study whose age lies in the interval [τ, τ + δτ] tends to ρc(t, x, τ)δV(x)δτ as δV, δτ → 0. Clearly, the physical dimension of the age distribution function is time −1 , and this function may be viewed as the histogram of the particle ages at time t and location x. Advection and diffusion proceed independently of the age of the particles being transported. Then, by having recourse to mass conservation considerations alone, Delhez et al. [59] established the equation governing the age distribution function, which reads: where ∇, v(t, x), and K(t, x) denote the del operator, the fluid velocity, and the diffusivity tensor, respectively. Under the Boussinesq approximation, the velocity is divergence-free, that is, ∇ · v = 0. As was argued in Appendix A of Deleersnijder et al. [69], the diffusivity tensor must be symmetric and positive-definite. Equation (1) holds valid for a passive or inert constituent (also termed tracer). By adding suitable production-destruction terms, it can be extended to take reactions into account [59]. Doing so is, however, not necessary to serve the purpose of the present study. The last term in Equation (1) is related to ageing. It may be seen as an advection term related to a unit velocity, representing the fact that the age tends to increase at the same pace as time progresses. The mass of the constituent under study that is present in elemental control volume δΩ(x) is obtained by taking the sum over all age categories, that is, where ρδV(x) is the mass of the fluid present in δΩ(x) under the Boussinesq approximation. Then, the concentration of the constituent under consideration, defined as a mass fraction (that is, a dimensionless variable), reads: which means that the concentration is the 0th order moment of the age distribution function. By integrating (1) over the age, Delhez et al. [59] obtained the well-known advection-diffusion equation governing the evolution of the concentration: The age content of a particle [69] is defined as the product of its age and mass. Like mass, this quantity is of an additive nature. As a consequence, the age content of the particles of the constituent under study that are present in δΩ(x) is: This points to the importance of the first-order moment of the age distribution function, which, in the CART-related literature, is termed age concentration. Delhez et al. [59] established the equation obeyed by α(t, x) by multiplying (1) by the age and integrating over the age, eventually yielding: Owing to its relation to the age content, the age concentration is an extensive variable. This is why it is no surprise that it satisfies a reactive transport equation. The last term in the right-hand side of (7) is related to ageing. It is through this term that the equations for the concentration and age concentration are coupled. Concentration and age concentration Equations (4) and (7) are of a parabolic nature [72]. Therefore, to solve each of them, the initial value of their solution must be prescribed and one boundary condition has to be enforced at every point of the surface delimiting the domain of interest.
According to the abovementioned age-averaging hypothesis [69], the mean age of the particles under study that are present in δΩ(x) is: This type of averaging is not the only one that can be conceived. The choice that has been made is the only arbitrary ingredient in the developments leading to CART's age. As opposed to the concentration and age concentration, the age is an intensive variable, rather than an extensive one. In contrast with the equations for C and α, the equation obeyed by the age, which is obtained by combining (1), (7) and (8), cannot be cast into a conservative form. This is why in most, if not all numerical studies, the mean age has been computed as the ratio of the age concentration to the concentration rather than by solving the age Equation (9) [39,47,65,69,[73][74][75][76][77][78][79][80][81][82][83][84][85][86][87][88]. This equation may, however, prove to be useful in theoretical studies [89][90][91]. In the right-hand side of relation (9), the first term is associated with ageing, whilst the second one bears some similarity with an advection term. However, the expression that may be regarded as the velocity, v − 2C −1 ∇C · K, is not divergence-free. Its behaviour is at the root of the intriguing symmetry that the age field occasionally exhibits [89,[91][92][93].
A number of studies focused on seawater, which may be split into several water types, components, or masses (though the word "type" is likely to be more appropriate in this instance), which can be treated as inert tracers [15,60,62,69,94,95]. Obviously, the diagnostic strategy must be designed in such a way that the sum of all the water type concentrations must be equal to a constant at any time and location. With no loss of generality, this constant may be assumed to be equal to unity. For a tracer with unit concentration, age Equation (9) simplifies to the equation solved by England [58], and the corresponding age is sometimes called ventilation, or ideal age.

Consistent Insulating, Departure, and Arrival Boundary Conditions
Numerically solving the equation for the distribution function presents several challenges. First, there is an additional independent variable, namely, the age (τ). For example, if the three space dimensions are taken into account, Equation (1) must be discretised in a five-dimensional space, the corresponding independent variables of which are t, x, y, z, and τ. The necessary numerical developments are not trivial and the added computational cost cannot go unnoticed [71]. Furthermore, when the distribution function exhibits a long tail, Cornaton [96] argued that standard numerical techniques are no longer appropriate, which is why he developed a method involving the Laplace transform and classical space-time discretisations. This technique is computationally efficient, but its implementation is not straightforward.
Most authors did not find it necessary to compute the age distribution function and, instead, focused on the (mean) age, obtained from the solution of the concentration and age concentration equations. These equations are to be solved under boundary conditions that must be consistent with each other. They should be derived from the boundary conditions that the age distribution function would obey if the equation governing this function was to be solved explicitly.
Accordingly, we will derive a number of consistent boundary conditions, illustrate their impact on relevant problems and, when appropriate, show that opting for inconsistent boundary conditions would have a detrimental impact on the age field. Schematically, we will address insulating, departure, and arrival boundaries. We will also consider gas exchanges through the water-air interface.
Though the illustrations below essentially deal with passive constituents, the developments leading to consistent boundary conditions apply to any type of constituent, which includes constituents undergoing reactions.

Insulating Boundary
Consider surface Γ, whose outward unit normal is denoted n. Assume that this surface is impermeable, thereby insulating the neighbouring part of the domain of interest from its environment. As a result, the velocity satisfies boundary condition Since this boundary is impermeable, no particles of the constituent under consideration, irrespective of their age, cross it, leading to Combining (10) and (11) yields To derive the boundary condition for the concentration, we integrate (12) over the age and use relation (3): Multiplying (12) by τ, integrating over τ and using definition (6) of the age concentration, we obtain Unsurprisingly, C(t, x) and α(t, x) obey zero normal diffusive flux boundary conditions. Boundary conditions (12)-(14) may be cast into generic form [(−K · ∇ς) · n] x∈Γ = 0, where ς(t, x) is the variable whose diffusive flux is zero through the boundary. If n is parallel to one of the principal axes of the diffusivity tensor, then such a boundary condition simplifies to[∇ς · n] x∈Γ = 0. For instance, assume that the diffusivity tensor takes the widely-used form of where κ h and κ v are the horizontal diffusivity and the vertical one, respectively. If the boundary under consideration is horizontal, then the aforementioned zero normal flux boundary condition degenerates into [∂ς/∂z] x∈Γ = 0, where z denotes the vertical coordinate. Now, for illustration purposes, assume that the domain of interest, Ω, is completely insulated from its environment, implying that the abovementioned boundary conditions are to be enforced on the whole boundary ( Figure 1). Further assume that the age of all the particles of a passive tracer is zero at the initial time, that is, c(0, x, τ) = C 0 (x)δ(τ), where C 0 (x) and δ(τ) are the initial concentration field and the Dirac delta function, respectively. The corresponding initial age concentration is readily seen to be α(0, x) = 0. Then, at any point of Ω and time t ≥ 0, the age distribution function is c(t, x, τ) = C(t, x)δ(τ − t), implying that the age concentration and the age are α(t, x) = tC(t, x) and a(t, x) = α(t, x)/C(t, x) = t. This result is readily understood. No particles enter or leave the domain. As all particles age at the same pace as time progresses, their age is equal to the elapsed time. The aforementioned results, α(t, x) = tC(t, x) and a(t, x) = t, can also be obtained without computing the age distribution function and, instead, by solving only the concentration and age concentration equations under boundary conditions (13) and (14). This is readily seen.  In this case study, that the age is equal to the elapsed time is of little diagnostic value. However, it may be regarded as a piece of information supporting the well-foundedness of boundary conditions (13) and (14). In other words, these boundary conditions allow us to obtain the expected result. Opting for another set of boundary conditions would lead to an unacceptable age field.
We should note in passing that, had we studied the evolution of a tracer undergoing a first-order decay process associated with constant mean life λ −1 (and, hence, half-life log2λ −1 ≈ 0.7λ −1 ), we would have obtained the following result: c d (t, x, τ) = e −λt c(t, x, τ), C d (t, x, τ) = e −λt C(t, x, τ) and α d (t, x, τ) = e −λt α(t, x, τ), where subscript "d" identifies the fields related to the decaying tracer. As a consequence, the age would have been unchanged: a d (t, x) = a(t, x). This is chiefly because first-order decay proceeds at the same rate at every time and position [19,97].

Departure Boundary
If the age is defined as the time elapsed since leaving a given boundary, then the age of all the particles under consideration must be set or reset to zero at the moment they touch this boundary, leading to where, in this subsection, Γ refers to the departure boundary (that is, the boundary where the age is prescribed to be zero), which is usually a fraction of the domain boundary; C Γ (t, x) is the tracer concentration on departure boundary Γ, which is assumed to be known (Dirichlet boundary condition). The boundary condition for the age concentration is derived from (15) as follows: Thus, the concentration and age concentration obey Dirichlet boundary conditions. To help understand the role of these boundary conditions, consider a one-dimensional flow, in semi-infinite domain x ∈ [0, ∞] ( Figure 2). The age of a passive tracer particle is defined as the time elapsed since leaving the departure boundary (x = 0). At the initial instant (t = 0), there is no tracer in the domain and the tracer concentration is prescribed to be equal to constant C Γ on the departure boundary. If positive constants U and K represent the velocity and the diffusivity, the equation obeyed by age distribution function c(t, x, τ) is ∂c ∂t This equation is to be solved under the following initial and boundary conditions Without any loss of generality, C Γ may be assumed to be equal to unity. Accordingly, the solution reads ( Figure 3) where Υ is the Heaviside step function, that is, function Υ(t − τ) is equal to unity (resp. zero) according to whether t > τ (resp. t < τ).
The concentration of the tracer is whilst its age concentration reads The related age, a(t, x) = α(t, x)/C(t, x), may be seen to be smaller than the elapsed time, as it should be. It is noteworthy that the (mean) age could have been obtained without explicitly calculating the age distribution function. The tracer concentration and age concentration are governed by the following partial differential problems: where the initial and boundary conditions are derived from (18). Then, lengthy manipulations would allow us to show that the solutions to (22) and (23) are (20) and (21), as expected. At first glance, it is not obvious that concentration (20) and age concentration (21) satisfy boundary conditions C(t, 0) = 1 and α(t, 0) = 0. To remove doubts, we prove in Appendix A that these boundary conditions are actually obeyed.
In the limit t → ∞, the age distribution function, the concentration, the age concentration, and the mean age are Thus, in the steady-state limit, the age is the time needed to travel distance x at speed U. Unfortunately, this simple result obscures the fact that, because of diffusion, the time taken for a given particle to travel from the entrance of the domain to a point located at distance x to the inlet is generally not equal to x/U, as is illustrated in Figure 3. If diffusivity K is zero, then the age distribution function is implying that the concentration and age concentration are Therefore, the age is x/U for t > x/U, and is undefined otherwise. These expressions can also be obtained by setting K = 0 in (22) and (23) and dealing with the resulting equations in the sense of distributions.
x x=0 Dirichlet or Robin boundary conditions departure boundary On this boundary, we can either prescribe that all particles of the tracer under study have zero age (Dirichlet boundary condition) or that the age of the incoming particles is zero, which leads to a Robin boundary condition, as seen in Section 3.3.
. The dimensionless age distribution functions are plotted at x * = 3 as functions of the age at different instants.

Departure Boundary: An Alternative Approach
Having recourse to Dirichlet boundary conditions is not the only option to deal with a departure boundary. An alternative approach consists in imposing the incoming tracer flux and prescribing that the age of the particles entering the domain is zero. Accordingly, on boundary Γ with outward unit normal n, the age distribution function must satisfy where Φ(t, x) is the incoming tracer flux (ms −1 ), which we assumed to be known. In this case, the age is the time elapsed since entering the domain rather than the time elapsed since leaving boundary Γ. Integrating (28) over the age, we obtain which simplifies to the boundary condition for the concentration, that is, As for the age concentration, we integrate the product of the age and relation (28), yielding Then, on Γ, the age concentration satisfies: The above relations (30) and (32) are Robin boundary conditions. With the latter, the (mean) age on boundary Γ is unlikely to be zero. This is because, on Γ, there are particles that are entering the domain (their age is zero) and also particles that have been in the domain for some time (their age is positive). We now revisit the illustration of the preceding subsection. The only modifications to be made are related to the boundary conditions at x = 0, which must be transformed to We set Φ = U so that the steady-state concentration is equal to unity as in the previous illustration. Doing so entails little loss of generality. Then, the age distribution function is The present age distribution function has a longer tail than that associated with a Dirichlet boundary condition at the incoming boundary of the domain (Figure 3). The steady-state concentration and age concentration are and As a consequence, the steady-state age reads The difference between the age ensuing from the Robin boundary conditions and that obtained under Dirichlet boundary conditions is equal to a constant, namely, K/U 2 , which, unsurprisingly, is an increasing (resp. decreasing) function of the diffusivity (resp. velocity). Needless to say, the same concentration and age concentration fields could have been obtained by directly solving the equations for the concentration and age concentration under the abovementioned initial and boundary conditions.

Arrival Boundary
The particles of the constituent under study cease to be taken into account (that is, they are discarded) at the moment they touch an arrival boundary. If Γ is a boundary of this type, then the age distribution function must satisfy Therefore, on Γ, the concentration and age concentration must be zero (Dirichlet boundary conditions): Some may be left unconvinced by the above reasoning, mainly because it leads to the requirement that the concentration be zero on the boundary, causing uncertainties as to how the mean age is to be evaluated on the boundary. In Section 3 of Delhez and Deleersnijder (2006) [98], detailed Lagrangian and Eulerian developments led to a similar result. Though these calculations were made in a different context, they may be found to be helpful to grasp the issue at hand.
On the arrival boundary, the mean age appears as the ratio of two functions, the age concentration and the concentration, that are both zero. However, the age is unlikely to be arbitrarily large, for it is precisely on this boundary that the particles under study cease to be taken into consideration. In addition, the age gradient may be seen to satisfy a property that is of use for graphical representations. First, we rewrite the equation governing the (mean) age, that is, relation (9), as follows: Thus, on the boundary under consideration, where concentration C is prescribed to be zero, this equation simplifies to ∇C · K · ∇a = 0. Since the concentration is zero on the boundary, its gradient must be normal to it and, hence, must be parallel to unit normal vector n, leading to [n · K · a] x∈Γ = 0. Then, it is readily seen that this expression is equivalent to a zero normal diffusive age flux boundary condition which, as pointed out in Section 3.1, simplifies to [∇a · n] x∈Γ = 0 if one of the principal axes of the diffusivity tensor is parallel to n. Clearly, (42) does not contradict the hypothesis that the age has a finite value on Γ. Departure and arrival boundary conditions (15), (16) and (38)- (40) have been derived without consideration of the direction the velocity on the boundary. However, it is likely that a diagnostic strategy will be built in such a way that a departure (resp. arrival) boundary will be an incoming (resp. outgoing) boundary, that is, a boundary on which v · n ≤ 0 (resp. v · n ≥ 0 ).
For illustration purposes, we revisit the one-dimensional problem dealt with in Section 3.2. We keep the departure boundary with Dirichlet boundary conditions at x = 0 and introduce an arrival boundary at x = L so that the domain of interest now has a finite length (0 ≤ x ≤ L) (Figure 4). For the sake of simplicity, we focus on steady-state solutions. Accordingly, concentration C(x) and age concentration α(x) obey and Assuming that the concentration is equal to unity on the incoming boundary entails no loss of generality. This is readily understood. Furthermore, dealing with similar idealised, one-dimensional problems have been found to be of use in previous studies [38,62].  The concentration and age concentration are ( Figure 5): and where dimensionless parameter Pe = UL/K is the Peclet number, that is, the ratio of the timescale characterising diffusion (L 2 /K) and that associated with advection (L/U). In the vicinity of the departure boundary (x = 0), the age, a(x) = α(x)/C(x), admits asymptotic expansion which, unsurprisingly, simplifies to a(x) ∼ x/U in the limit Pe → ∞. As for the arrival boundary (x = L), the age tends, as expected, to a finite value with a zero gradient, with a(L) → L/U as Pe → ∞. The larger the Peclet number, the closer the solutions are to their zero diffusion counterparts, that is, a unit value of concentration, with the age concentration and age equal to x/U. Such solutions cannot satisfy the boundary conditions prescribed at x = L. This is why the correct concentration and age concentration exhibit a boundary layer adjacent to the outgoing boundary.  (46), and (c) the associated age. Panel (d) depicts the inconsistent age that is obtained as the ratio of inconsistent age concentration (49) and correct concentration (49). Dimensionless variables are represented. They are identified by asterisks and are x * = x/L, α * = Uα/L, a * = Ua/L and a * inc = Ua * inc /L.
To document the impact of inconsistent boundary conditions, we replace for a moment the Dirichlet boundary condition for the age concentration at x = L by Neumann condition [−K d α inc /dx] x=L = 0, where subscript "inc" refers to the inconsistent solution. This relation is the simplest type of radiation condition, which may be found in Table 1 of Bendsten et al. [65] in conjunction with a Dirichlet boundary condition for the concentration. Obviously, the concentration is not changed. The inconsistent age concentration reads The modified age is slightly different from the correct one in the neighbourhood of the incoming boundary, but has no finite limit on the outgoing boundary, In other words, the age resulting from the imposition of inconsistent boundary conditions on the downstream boundary is such that a inc → ∞ as x → L, which is unjustifiable in a finite-sized domain of interest with an open boundary meant to allow particles to leave the domain. What is worse, the inconsistent boundary conditions impact a significant fraction of the domain of interest ( Figure 5). Undoubtedly, such behaviour is unacceptable.

Arrival Boundary: An Alternative Approach
Imposing, as suggested in the previous Section, that the concentration be zero on the arrival boundary implies that the advective mass flux through this boundary is zero. Thus, the outgoing mass flux crossing the boundary is purely diffusive. An alternative approach consists in prescribing that the diffusive flux through the boundary be zero so that the outgoing flux is entirely of an advective nature. This requires boundary conditions equivalent to (11), (13), and (14) to be enforced on the arrival boundary.
To illustrate the impact of these boundary conditions, the idealised problem of the previous Section is revisited. The only modification pertains to the outgoing boundary, where the diffusive fluxes are prescribed to be zero, yielding This age is larger than that obtained by imposing Dirichlet boundary conditions at the departure boundary ( Figure 6). This can be proven rigorously for the present one-dimensional flow. Remarkably, a similar inequality also holds valid for a much wider class of problems [99,100].   Figure 4. Panel (a) depicts the age ensuing from Dirichlet boundary conditions imposed at the departure boundary (x = L) and, hence, is the same age as that represented in panel (c) of Figure 5. The age obtained by prescribing Neumann boundary conditions at x = L, that is, age (52), is illustrated in panel (b). Dimensionless variables are represented. They are identified by asterisks and are defined as follows: x * = x/L and a * = aU/L.

Gas Exchanges through the Water-air Interface
The gas flux at the water-air interface is usually parameterised by having recourse to the concept of piston velocity [101][102][103][104][105]. Accordingly, the net outgoing (that is, from the water body to the atmosphere) mass flux (kg m −2 s −1 ) through surface Γ of a gas whose concentration in the water is where is the piston velocity (ms −1 ), whilst C s is the saturation concentration, that is, the water surface concentration in equilibrium with the atmosphere. In general, and C s depend on time and position. If the surface concentration is greater (resp. smaller) than the saturation concentration, the gas flux is directed from the water body to the atmosphere (resp. from the air to the water). This is why the boundary condition (53) may be viewed as a relaxation boundary condition-the air-water mass flux tends to nudge the surface concentration toward C s . Formula (53) provides an estimate of the net flux at the water-air interface. No other information about the water-air interface is needed in order to model the concentration in the water (or in the atmosphere) of the gas under consideration. For age calculations, however, it is necessary to realise that the right-hand side of (53) actually represents the difference between the upward mass flux (ρφ ↑ C = [ρ C] x∈Γ ) and the downward one (ρφ ↓ C = [ρ C s ] x∈Γ ) (Figure 7). This piece of information is essential for building the boundary condition for the age distribution function. The outgoing (resp. incoming) mass flux of the gas particles whose age lies in the interval where c a is the age distribution in the atmosphere. Thus, the boundary condition for the age distribution function reads Under the assumption that the integral of c a over the age τ is C s , the boundary conditions for the concentration and age concentration are obtained by taking the 0th-and first-order moment of (54) ( Table 1): Obviously, relations (53) and (55) are equivalent. Equations (54)-(56) are usually referred to as Robin boundary conditions. It is often assumed that all the gas particles entering the water column have the same age, a a (t, x), that is, the gas age at the lower boundary of the atmosphere. As a result, their age distribution function is c a = C s δ(τ − a a ). Therefore, in the water body, the age of the dissolved gas is the sum of atmospheric age a a (t, x) and the time spent in the water body since entering it through the water-air. For ventilation studies, it may be appropriate to assume that the atmospheric age is zero (Table 1). Table 1. Outgoing and incoming specific fluxes (ratio of a flux to the water density) at the water-air interface for the age distribution function (φ ↑,↓ c ), the concentration (φ ↑,↓ C ), and the age concentration (φ ↑,↓ α ). As for the downward flux, the general expression and two simplified ones are taken into consideration, which consists in assuming that all the incoming gas particles have the same age, a a (fourth column), and that this age is zero (fifth column).

Variable
Outgoing (Upward) Incoming (Downward) Specific Flux Specific Flux General Expression c a = C s δ(τ − a a ) c a = C s δ(τ) If the piston velocity is small, then the boundary conditions derived above tend to simplify to no-flux expressions relevant to an insulating boundary. If, on the other hand, the piston velocity is large, (54)-(56) tend to degenerate into Dirichlet boundary conditions. To assess the impact of the piston velocity, a dimensionless parameter should be derived. This can be achieved with the help of a steady-state water column model (Figure 8a). The domain of interest is defined by inequalities −h ≤ z ≤ 0, where z is the vertical coordinate and h is the height of the water column, whose water-air boundary is located at z = 0. As is customary in water column modelling, the lower boundary of the domain is considered to be impermeable. For the sake of simplicity, it is assumed that vertical diffusion, represented by means of constant diffusivity K, is the only process to be taken into account. Accordingly, the concentration and age concentration obey the following differential problems: and Unsurprisingly, the concentration is a constant: C(z) = C s . The age concentration and age are readily seen to be and a(z) = a a + − (2h + z)z 2h where is the sought-after dimensionless parameter (Figure 8). The latter is the ratio of the timescale associated with the piston velocity, h/ , and the classical diffusive timescale, h 2 /K, that is, The larger the piston velocity, the smaller this dimensionless parameter. In the limit → 0, the age tends to the function of z that would be obtained under Dirichlet boundary conditions. The age is the sum of the atmospheric age and a 0 (z), which is the age ensuing from the assumption that the atmospheric age a a is zero. Interestingly, this age shift property is satisfied in a wide class of multi-dimensional, time-dependent age calculation problems [106].
In reality, as opposed to what is represented in the above highly idealised water column model, the impact of vertical turbulent diffusion due to the surface forcing (chiefly surface wind stress) does not always extend to the bottom. Therefore, the relevant vertical length scale is not necessarily the height of the water column. It could be significantly smaller. A plausible option consists in selecting the thickness of the surface mixed layer, when such a hydro-dydnamic feature can be identified. Furthermore, the vertical eddy diffusivity and the piston velocity are likely to be timeand position-dependent, implying that their typical order of magnitude should be evaluated and subsequently introduced into (61). Accordingly, the final formulation of dimensionless parameter (61) is where K, H, and W denote a typical value of the vertical eddy diffusivity near the water-air interface, the relevant vertical length scale, and the order of magnitude of the piston velocity, respectively. For advective and diffusive transport problems, Haine [104] studied the relationship between solutions obtained under Dirichlet boundary conditions and those ensuing from Robin conditions. The obtained theoretical results are rather general, but are beyond the scope of the present study.
There is no denying that the developments related to gas exchanges through the water-air interface are somewhat cumbersome. If obtaining diagnostic quantities for such phenomena were to prove too laborious, we may wonder if it would be worth the effort, suggesting that we should perhaps restrict ourselves to the simplest approach, which consists in assuming that the atmospheric age is zero, as laid out in the rightmost column of Table 1.

A Simple Ventilation Assessment Problem
Each of the illustrations included in the preceding Section was meant to help comprehend the role of a single type of boundary condition. To gain further insight into the role of boundary conditions in age calculations, it may be desirable to consider a slightly more sophisticated situation. Seeking inspiration in [65,107], we will tackle a relatively simple, two-dimensional, "horizontal-vertical" ventilation assessment study with constant hydro-dynamic parameters.
Let x and z denote the horizontal coordinate and the vertical one, respectively. The domain of interest is defined by inequalities 0 ≤ x ≤ L and −h ≤ z ≤ 0 (Figure 9). Water flows in the direction of increasing x with constant (horizontal) velocity U. Horizontal and vertical diffusion is taken into account with the help of constant diffusivities K x and K z . The bottom (z = −h) is an insulating boundary. The vertical boundaries located at x = 0 and x = L are open. The latter is an arrival boundary and the former is a departure one, and so is the water-air interface. Therefore, the age is the time elapsed since leaving the incoming boundary (x = 0) or the water-air interface (z = 0). For the sake of simplicity, only steady-state solutions will be considered. Accordingly, the concentration of the ventilation tracer, C(x, z), is the solution of the following partial differential problem: Then, the associated age concentration, α(x, z), obeys Finally, the age is a(x, z) = α(x, z)/C(x, z).
To the best of our knowledge, there exists no simple analytical solution to the differential problem (63)- (64). This is why we built a numerical solution for it by having recourse to the tracer transport module of the finite-element, discontinuous Galerkin model, SLIM (www.slim-ocean. be) [108][109][110]. The computations are based on a fine mesh consisting of 39,204 rectangular elements. The resolution is increased near the domain boundaries in such a way that the discrete solution is believed to be very close to the exact one.
There are two crucial dimensionless numbers associated with the present transport problem. The horizontal (resp. vertical) Peclet number is the ratio of the horizontal (resp. vertical) diffusion timescale L 2 /K x (resp. h 2 /K z ) to the advective timescale L/U, yielding Pe x = UL/K x (resp. Pe z = h 2 U/(LK z )). Since the vertical velocity is zero, introducing a vertical Peclet number may seem to be somewhat questionable. It is argued in Appendix C that the aforementioned vertical Peclet number is, roughly speaking, in line with common practice.
In shallow domains, such as rivers, estuaries, or coastal regions, vertical (turbulent) diffusion plays an important role. Therefore, assuming that the vertical Peclet number ranges from 1 to 100 with a typical order of magnitude of 10 would presumably be acceptable. In the aforementioned domains of interest, the aspect ratio (that is, the ratio of the vertical length scale to the horizontal one) is significantly smaller than unity, which is why the horizontal Peclet number should probably be somewhat larger than the vertical one. The concentration, age concentration, and age are displayed in Figure 10 for (Pe x , Pe z ) = (10, 10) and (Pe x , Pe z ) = (100, 10). The boundary conditions at the water-air interface mostly impact the solution in the upper part of the domain, so that the maximum of the age is always attained at (x, z) = (L, −h), that is, at the lower-right corner of the domain of interest. The smaller the horizontal Peclet number, the larger the impact of the boundary conditions prescribed at the upper boundary of the domain. Clearly, the present age distribution is more complex than that obtained in Section 3.4, though there, the horizontal processes taken into account are similar. Unsurprisingly, in the vicinity of the lower boundary of the domain (z = −h), with relatively large values of Pe z , the age tends to be similar to that derived from (45)- (46). Pe x = 10 and Pe z = 10 Pe x = 100 and Pe z = 10 Concentration Age concentration Age Non consistent age Figure 10. Illustration of the concentration, age concentration, and age from the solution of the partial differential problem (63) and (64) for (Pe x , Pe z ) = (10, 10) and (Pe x , Pe z ) = (100, 10). Dimensionless variables are represented. They are identified by asterisks and are defined as follows: x * = x/L, z * = z/h, α * = Uα/L, and a * = Ua/L. The non-consistent age ensuing from inappropriate boundary conditions for the age concentration on the departure boundary is represented in the lowermost row of the graph.
Finally, on the departure boundary (x = L), we introduce an inconsistent boundary condition for the age concentration similar to that of Section 3.4, that is, [−K ∂ α inc /∂x] x=L = 0. This causes the age to be infinite on the departure boundary. Nonetheless, as opposed to what was observed in the one-dimensional solution of the abovementioned Section, in the two-dimensional problem, the error does not affect a large fraction of the domain, because of the impact of the correct boundary conditions at the water-air interface. This is why we suspect that the studies that relied on non-consistent boundary conditions of the type referred to here are plagued by errors that, though non-negligible, are not catastrophic.
As illustrated by Figure 11, the error due to the use of non-consistent boundary conditions occurs in a region adjacent to the arrival boundary, whose width increases as the distance to the water-air interface increases. For the purpose of a sensitivity analysis, it is convenient to evaluate the maximum width of the region where the error is significant. A simple measure (Λ) is defined as follows: a inc (Λ, −h) = a(L, −h). This width, which is represented in Figure 12, is a decreasing (resp. increasing) function of the horizontal (resp. vertical) Peclet number. This is in agreement with elementary physical intuition. Figure 11. Difference between the non-consistent age and the correct one for (Pe x , Pe z ) = (10, 10). Dimensionless variables similar to those of Figure 10 are represented. On the departure boundary, the age difference is infinite, but the colour is saturated at a value of 10.

Discussion and Conclusions
The (mean) age of a constituent of the water, or a group of constituents, including the water itself (that is, the aggregate of all the constituents), is a diagnostic timescale depending on time and position that can be obtained from the solution of a system of partial differential equations. The general form of them has been well-known since the turn of the century [59,69]. Over the past two decades, relatively little attention has been devoted to the construction of the initial and boundary conditions under which the age-related equations are to be solved. This is, however, a critical ingredient of an age-based diagnostic strategy-for the solutions of age-related (or any other) parabolic partial differential equations to be unambiguously determined, the initial and boundary conditions must be precisely defined. In this regard, casualness must be ruled out, as has been exemplified above by documenting the detrimental impact of some inconsistent boundary conditions. While initial conditions are rather easily built, which is why we have not tackled them explicitly, constructing appropriate boundary conditions is less straightforward, as has been shown in the present study. Hopefully, the latter will help clarify the methodology to set up an age-based diagnostic approach. The first steps of it should be as follows: 1.
Set out the reasons why the age, rather than other timescales (or diagnoses of another nature), is likely to be of use to help interpret the aquatic processes under consideration; 2.
Select the constituent whose (mean) age is to be evaluated and explain the rationale of this choice; 3.
Define the age, especially where and when the age of a particle of the constituent under study is to set or reset to zero, as well as where, when, and how this particle will cease to be taken into consideration; 4.
Build the boundary conditions for the age distribution function in accordance with the outcome of the previous three steps; 5.
Derive consistent boundary conditions for the concentration and age concentration using the methodology developed in this article (see also Appendix D).
Obviously, the following steps will consist in solving, analytically or numerically, the relevant partial differential problems and discuss the obtained results, which cannot be achieved in a fruitful manner if the foundations of the diagnostic approach are shaky.
Steps 1 to 3 above seem to be rather straightforward, if not trivial. However, as was seen by Delhez et al. [111], an ill-conceived diagnostic strategy may lead to the evaluation of timescales contradicting their very definition, eventually leading to dubious interpretations. Clearly, we should bear in mind the wise piece of advice of Bolin and Rodhe [55]: "To avoid misunderstandings and even erroneous conclusions, it is important to introduce precise definitions and to use them with care". It is because this word of caution has been overlooked time and again that Viero and Defina [41] referred to the field of diagnostic timescales as a modern Tower of Babel. Hopefully, the present study will be considered as a modest contribution to the deconstruction of this edifice. We strongly believe that all the developments made above are also relevant to partial ages, a recently developed generalisation of the concept of age [76,112]. This is because there is no fundamental difference between the concept of age and that of partial age: similar lines of argument should apply to both types of diagnoses.
It is impossible to address all the existing types of boundary conditions in a single paper. However, the approach advocated herein (that is, deriving the boundary conditions for the concentration and age concentration from those relevant to the age distribution function) can probably be applied to open boundary conditions other than those dealt with above, in particular, tracer-adapted versions of radiation conditions [113][114][115][116], sponge layers [117,118] and other techniques [119]. This has yet to be convincingly substantiated.
Clearly, dealing with realistic case studies is beyond the scope of the present article, for its key objective is the development of the theory to build consistent age-related boundary conditions. However, the boundary conditions used in the idealised ventilation rate assessment of Section 4 are most likely to be similar to those that would be implemented in realistic ventilation studies. This is especially true for the surface and bottom boundary conditions. Many ventilation studies [57,58] did not have to cope with lateral open boundaries. Their lateral boundaries were insulating ones, leading to trivial boundary conditions as may be seen in Section 3.1. In general, the open boundaries are believed to be the most problematic ones. To estimate water renewal of semi-enclosed domains, many authors resorted to Dirichlet boundary conditions [62,120] on open boundaries. This approach is undoubtedly the easiest one. However, we are convinced that, in the near future, more subtle boundary conditions will be worked out, involving fluxes rather than prescribed values, which will be related, in one way or another, to the boundary conditions for the momentum equations. Such boundary conditions are being developed and will be described in forthcoming articles.
Although the present study focuses on Eulerian developments, it must be realised that, in principle, all the age calculations, for idealised or realistic flow processes, can be achieved by means of Lagrangian methods, as well as Eulerian ones-as explained in van Sebille et al. [68] and some of the references therein, both approaches should converge to similar solutions. In the Lagrangian framework, it is not necessary to explicitly evaluate the concentration and the age concentration, which is why all of the above developments about boundary conditions are likely to be irrelevant to Lagrangian techniques. However, Lagrangian calculations have disadvantages-the representation of diffusive processes by means of stochastic terms is not straightforward [121,122] and the necessary computational resources are unlikely to be smaller than those required for carrying out Eulerian computations. In addition, most hydrodynamic models are Eulerian, making it somewhat easier to implement diagnoses rooted in the same framework. Finally, properties of diagnostic timescales are generally simpler to derive in the Eulerian framework than in the Lagrangian one [62,69,89,91,97,100]. All this being said, literal interpretations are much easier to produce in the Lagrangian framework, which is the reason why, in the present article, as well as in many publications of the diagnostic timescale literature, the theoretical developments and calculations are Eulerian, whereas the explanations and interpretations resort to a vocabulary rooted in the Lagrangian formalism.
The present study focused on the calculation of the age of a tracer as a means to help understand complex aquatic fluid flows. However, as underscored by Wunsch [123], there are many more aspects in tracer and timescale methods, the significance of which cannot be overestimated (e.g., time dependency, the number of space dimensions, Green's function theory, inverse methods, stochastic boundary conditions, and timescales other than those of CART or similar ones). In addition, the aforementioned article introduces a number of analytical solutions that should no longer be overlooked. In this respect, the reference to Carslaw and Jaeger [124] is a very important one. Clearly, much more work is needed in this rich field of research.

Conflicts of Interest:
The authors declare no conflict of interest.
Wh/K z . Since the Boussinesq approximation is assumed to hold valid, the velocity is divergence-free, implying that the order of magnitude of the vertical velocity satisfies W/h = U/L. Combining this expression with the vertical Peclet number defined above, we would obtain Pe z = h 2 U/(LK z ). Thus, though it is not immediately self-evident, the Peclet number introduced in Section 4 is in line with the widely-used definition of this dimensionless number.

Appendix D
All the boundary conditions for the age distribution function set out above are special cases of the following generic form: where b(t, x) is an appropriate vector, whilst η(t, x) and µ(t, x, τ) are scalar functions, whose specific form depends on the nature of the boundary condition to be built. To derive the boundary conditions for the concentration and age concentration, we integrate Equation (A7) over the age, yielding b · ∇C + ηC + If the age is prescribed to be zero on the boundary (a Γ = 0), then (A11) simplifies to [α(t, x)] x∈Γ = 0. Thus, we have obtained a pair of Dirichlet boundary conditions that can be applied on a departure boundary (see Section 3.1). The simplest arrival boundary conditions, that is, Equations (39) and (40), may be derived from (A7)-(A9) by setting b(t, x) = 0, η(t, x) = 1 and µ(t, x, τ) = 0 (see Section 3.4).