GOAL-ORIENTED A POSTERIORI ERROR ESTIMATION FOR THE TRAVEL TIME FUNCTIONAL IN POROUS MEDIA FLOWS

Abstract. In this article we consider the a posteriori error estimation and adaptive mesh refinement for the numerical approximation of the travel time functional arising in porous media flows. The key application of this work is in the safety assessment of radioactive waste facilities; in this setting, the travel time functional measures the time taken for a non-sorbing radioactive solute, transported by groundwater, to travel from a potential site deep underground to the biosphere. To ensure the computability of the travel time functional, we employ a mixed formulation of Darcy’s law and conservation of mass, together with Raviart-Thomas Hpdiv,Ωq-conforming finite elements. The proposed a posteriori error bound is derived based on a variant of the standard Dual-WeightedResidual approximation, which takes into account the lack of smoothness of the underlying functional of interest. The proposed adaptive refinement strategy is tested on both a simple academic test case and a problem based on the geological units found at the Sellafield site in the UK.


Introduction.
In recent decades the use of numerical simulations in hydrogeological studies has become commonplace across a range of applications.Amongst these, modelling the post-closure safety performance of deep geological storage of radioactive waste is of particular interest for a posteriori error estimation.Efficient and reliable simulations are required in order to assess the suitability of a specific location for siting a waste repository.Furthermore, there is a critical need to verify any computational results with rigorous error bounds as the effects of an inaccurate simulation could be extremely costly.In the safety assessment of radioactive waste facilities, one of the key quantities of interest is the time taken for a non-sorbing radioactive solute, transported by groundwater, to travel from a potential site deep underground to the biosphere [28,41,44].Additionally, accurate computation of the travel time has applications in streamline methods for modelling other subsurface flows; for instance in oil and gas reservoir management [36].
The suitability of finite element methods has been demonstrated for many of the complex geometries and physical effects that are associated with the numerical approximation of groundwater flow and contaminant transport problems [16,17,20,21,37,41,42].There are, however, problems associated with the use of nodal-based elements, such as lack of mass conservation at an elemental level and unphysical streamlines, as noted in [21,24,42,43].These problems are not observed when using a (conforming) mixed finite element method, or finite difference method, in which pressure and velocity are computed simultaneously.Indeed, when employing the latter class of methods on triangular meshes, the tracing of streamlines is straightforward and has been shown to yield physical results, even on very coarse meshes, when nodalbased approximations may fail [24,37,42].While these works use only the lowest ˚School of Mathematical Sciences, University of Nottingham, University Park, Nottingham NG7 2RD, UK, email: andrew.cliffe@nottingham.ac.uk.
order approximation for the Darcy velocity, it is also possible to compute streamlines to high-order, for divergence-free flows, using a mixed finite element method [36].
A posteriori estimation is a key tool for the control of the discretization error arising in numerical simulations involving the finite element method [1,3,4].Traditionally, a posteriori error analysis has focussed on estimating the error in global quantities of the solution.Goal-oriented a posteriori error estimation seeks to determine whether some physical quantity of interest, calculated using the numerical solution, is within a given tolerance of the true value.Goal-oriented techniques for a posteriori error estimation were introduced in [6,25]; see also [7,8,26,29,30,31,32,33,34,38,39,45], and the references cited therein.The formulation in [6] constructs an adjoint-weighted error bound, whilst the alternative formulation in [25] constructs an unweighted error bound based on exploiting stability estimates.Adjoint-weighted residual techniques have been applied to a wide range of physical problems including, amongst others, the Poisson problem [6], elasticity theory [45], incompressible viscous flows [7,Chp. 8], radiative transfer [7,Chp. 11], and non-linear hyperbolic conservation laws [30].To date, most work has focussed on single target functionals, where, typically, the functional is either linear or non-linear, but Fréchet differentiable.Simultaneously, a posteriori estimation of the error in a given set of functionals has been undertaken in the articles [33,29].Finally, we point out that, in the context of our current article, some limited work has been undertaken in the context of a posteriori error estimation for the accurate computation of streamlines and the travel time.Indeed, the article [17] exploits an energy norm error estimate to adaptively refine the computational mesh in order to increase the accuracy of the computation of streamlines using a stream function approach.However, whilst there is literature on the application of adjoint-weighted techniques for contaminant transport [8,39], there is, to the authors' knowledge, no literature in which the travel time is directly considered as the quantity of interest.
The main contribution of this article is the introduction of a goal-oriented a posteriori error estimate for the travel time in groundwater flows.This is a key quantity of interest in many physical simulations; in particular, post-closure assessment of deep geological radioactive waste repositories, which has received little attention in the error estimation literature.Furthermore, there are a limited number of results in the adjoint-weighted goal-oriented a posteriori error estimation literature concerning non-Fréchet differentiable functionals.This article introduces a goal-oriented a posteriori error estimate that does not explicitly require Fréchet differentiability of the quantity of interest.Moreover, we apply this general approach to compute the travel time for a simplified model based on the geological units found at the Sellafield site in the UK.
This article is organized as follows.In section 2 we introduce the classical mixed formulation of Darcy flow and define the travel time functional.In section 3 we pose the weak formulation of the Darcy flow from section 2 and introduce a suitable mixed finite element approximation.In section 4 we establish goal-oriented a posteriori error estimates for linear and non-linear Fréchet differentiable functionals and adapt these ideas to derive an estimate for the error in the travel time functional.In section 5 we present a series of numerical experiments to demonstrate the efficiency of the proposed a posteriori error estimator.Finally, in section 6 we make some concluding remarks and highlight ongoing and future work.
2. Problem Definition.In this section we introduce the first-order system of partial differential equations (PDEs) that describe saturated groundwater flow.We then define the travel time functional for a non-sorbing solute in a contaminated groundwater plume; an a posteriori estimate is derived in section 4 for this quantity of interest.

Primal PDE Problem.
Let Ω be an open, bounded, Lipschitz domain in R 2 with polygonal boundary BΩ.Let BΩ be partitioned into two open subsets Γ D and Γ N , such that ΓD Y ΓN " BΩ and Γ D X Γ N " H.We model groundwater flow in Ω using Darcy's law and conservation of mass, cf.[5,23], to yield a Poisson-type boundary value problem for the Darcy velocity u and hydraulic head H, given by u `K∇H " 0 where n denotes the unit outward normal to BΩ and K is the hydraulic conductivity.
We specify that K is a 2 ˆ2 positive definite matrix and its smallest eigenvalue is bounded away from zero uniformly in Ω with respect to x, i.e., there exists α 0 ą 0 such that α 0 |ξ| 2 ď ξ Kξ @ξ P R 2 , @x P Ω, which implies that K is invertible almost everywhere in Ω.However, K may contain discontinuous coefficients, which are potentially highly anisotropic.Furthermore, we make the following assumptions regarding the model data: pΓ D q, and g N P L 2 pΓ N q.
For the numerical experiments presented in this article, we consider problems which possess homogeneous Neumann boundary conditions for the Darcy velocity on certain portions of the boundary, due to the geology and tomography of the region, and inhomogeneous Dirichlet boundary conditions for the hydraulic head on the remainder of the boundary.As such, we impose the condition that g N " 0 throughout this article.However, we do this without loss of generality, since the numerical treatment of the homogeneous case does not differ from that of the inhomogeneous one, cf.[12, Chp.IV.1].

Quantity of Interest.
The physical quantity of interest we consider in depth in this article is the travel time for a non-sorbing solute from some release point x 0 P Ω to the boundary BΩ of the computational domain.In the context of deep geological storage of radioactive waste, this corresponds to the time taken for leaked radioactive material to be transported from the repository to the biosphere.
As the fluid is flowing only through the connected pores in the solid matrix, the Darcy velocity does not represent the fluid velocity.This transport velocity is instead given by where φ is the porosity of the rock.
If we make the assumption that there are no dispersion or sorption effects, then, since the Darcy velocity is not time-dependent, the motion of a particle simply corresponds to streamlines of the flow.This implies that the particle position X, through the domain, is governed by the differential equation subject to the initial condition Xp0q " x 0 . (2.3) If we denote the path taken through the domain by the particle as P, then the travel time T p , is given by where } ¨}2 is the Euclidean norm in R 2 .This then allows us to define the travel time functional J p¨q, by where Ppvq is the path obtained for the transport velocity given by v T " 1 φ v. Remark 2.1.The dependence on v in the computation of the path over which the integral is evaluated, in addition to the term in the integrand, leads to significant challenges in demonstrating Fréchet differentiability of this functional with respect to v.In fact, it is clear that this functional is not globally continuous and, therefore, not globally Fréchet differentiable.However, it remains unclear as to whether or not, for certain data, the functional is Fréchet differentiable, or Gâteaux differentiable, on some open neighbourhood of the solution, in a suitable Hilbert space.With this in mind, in section 4.3 we adapt the standard goal-oriented a posteriori error estimates for Fréchet differentiable functionals to overcome these difficulties, subject to certain restrictions that are discussed later.

Finite Element Discretization.
In this section we define the mixed weak formulation of the first-order system of PDEs introduced in section 2. Given this weak formulation, we then define a suitable finite element subspace on which the discrete weak formulation yields a convergent numerical approximation.
3.1.Weak Formulation.We first introduce the function spaces that will be used in the weak formulation and definition of the finite element discretization.To this end, we define the Sobolev space Hpdiv, Ωq by Hpdiv, Ωq " cf. [12].Let x¨, ¨yBΩ denote the duality pairing between H 1 2 pBΩq and H ´1 2 pBΩq.Then the Green's formula allows us to define the normal trace of a function v P Hpdiv, Ωq on BΩ [12, Chp.III.1].We now proceed by defining the subspace H 0,N pdiv, Ωq Ă Hpdiv, Ωq of functions whose normal trace vanishes on Γ N by H 0,N pdiv, Ωq " v P Hpdiv, Ωq : xv ¨n, Φy BΩ " 0 @Φ P H 1 0,D pΩq where Given the definition of H 0,N pdiv, Ωq in equation (3.2), we now state the weak formulation of the first-order system of PDEs stated in equation (2.1).Multiplying the first and second equations in (2.1) by test functions in H 0,N pdiv, Ωq and L 2 pΩq, respectively, and applying the Green's formula given in equation (3.1) to the first equation in (2.1) yields the weak formulation: find u P H 0,N pdiv, Ωq and H P L 2 pΩq such that apu, vq `bpv, Hq " Gpvq @v P H 0,N pdiv, Ωq, bpu, qq " F pqq @q P L 2 pΩq.
The bilinear forms ap¨, ¨q and bp¨, ¨q are defined as apv, wq " ż Ω K ´1v ¨w dx, v, w P Hpdiv, Ωq, and bpv, qq " ´żΩ q div v dx, v P Hpdiv, Ωq, q P L 2 pΩq, respectively, and the linear functionals Gp¨q and F p¨q are defined as Gpvq " ´xv ¨n, g D y BD , v P Hpdiv, Ωq, and F pqq " ´żΩ f q dx, q P L 2 pΩq, respectively.
It is well known that the weak formulation given in equation (3.3), subject to the restrictions on the data outlined in section 2, possesses a unique solution [12, Chp.II.1].

Mixed Finite Element
Formulation.We consider shape-regular (conforming) meshes T h that partition Ω Ă R 2 into disjoint open triangular domains κ, such that Ω " Y κPT h κ.Here, we implicitly assume that T h respects the decomposition BΩ = ΓD Y ΓN of the boundary, in the sense that an edge of a boundary element κ is solely contained within ΓD or ΓN .Given κ P T h , we write Bκ to denote the boundary of κ; the outward unit normal vector to Bκ is given by n κ .By h κ we denote the element diameter of κ P T h and introduce the mesh function h, defined by h " maxth κ : κ P T h u.For k P N 0 and κ P T h , we write P k pκq to denote the space of polynomials of degree at most k on κ.With this notation, we introduce the Raviart-Thomas Element [46] given by It can be shown, cf.[12], for example, that v P RT k pκq may be fully characterized by the moments of up to order k of v ¨nκ on the edges of κ and the moments of up to order k ´1 of v on the interior of κ.Given the definition of RT k pκq, and the characterization of its degrees of freedom, we may now build a finite-dimensional subspace of Hpdiv, Ωq from the spaces RT k pκq, κ P T h .To this end, we define the following finite dimensional subspace of Hpdiv, Ωq by RT k pΩ, T h q " tv P Hpdiv, Ωq : v| κ P RT k pκq @κ P T h u , cf. [12].Thereby, setting V 0 h,k " tv P RT k pΩ, T h q : v P H 0,N pdiv, Ωqu and Q h,k " p P L 2 pΩq : p| κ P P k pκq @κ P T h ( , the discrete mixed formulation is defined by: find It can be shown, cf.[10,12], that for a given k ě 0 the choice of the finite element spaces V 0 h,k and Q h yields a convergent approximation defined by the discrete mixed formulation given in equation (3.4).
We stress that the exploitation of the mixed formulation introduced in section 3.1 for the first-order system of PDEs given in equation (2.1) and the subsequent discretization based on employing Hpdiv, Ωq-conforming elements is crucial for the computation of streamlines and travel times.Indeed, if the system is solved as a single second-order Poisson-type PDE using standard H 1 pΩq-conforming or discontinuous Galerkin finite element methods, for example, then there is no guarantee that the normal trace of the numerical Darcy velocity has consistent sign with respect to adjacent elements; this may then lead to non-physical streamlines and non-computable travel times.The continuity of normal traces is a key property of Hpdiv, Ωq-conforming finite elements that is demonstrated in the following lemma.Lemma 3.2.Let Y pΩ; T h q Ă `L2 pΩq ˘2 be the broken Sobolev space defined by Y pΩ; T h q " !v P `L2 pΩq ˘2 : v P Hpdiv, κq @κ P T h ) .
Proof.See [12, Chp.III.1] Remark 3.3.We point out that alternative choices of Hpdiv, Ωq-conforming elements, such as BDM [11], may also be employed.Indeed, the proceeding a posteriori error estimation naturally holds for any appropriate choice of stable Hpdiv, Ωqconforming finite element spaces.

A Posteriori Error Estimation.
In this section we follow the standard theory of adjoint-weighted goal-oriented a posteriori error estimation to derive error estimates for both linear and non-linear Fréchet differentiable functionals.We then present the a posteriori error estimate for the travel time functional based on employing a one-sided difference approximation of the Gâteaux derivative of the functional.Finally, we discuss the limitations on the applicability of this estimate with regards to round off error and the truncation error in a generalized Taylor expansion.
Similarly, the (primal) finite element approximation may be written in the equivalent form: find , respectively, the following Galerkin orthogonality property holds We may now proceed as in [4,7] by defining the following adjoint problem, corresponding to the quantity of interest J, as follows: find rz, rs P H 0,N pdiv, Ωq ˆL2 pΩq such that A `rv, qs, rz, rs ˘" J `rv, qs ˘@rv, qs P H 0,N pdiv, Ωq ˆL2 pΩq. (4.3) Exploiting the linearity of J and the Galerkin orthogonality property (4.2), we may derive the following error representation formula for the error in the linear functional J: for all rz h,k , r h,k s P V 0 h,k ˆQh,k .Here, Rp¨, ¨q is referred to as the weighted residual, defined by for α, β P H 0,N pdiv, Ωq ˆL2 pΩq.The error representation formula (4.4) requires knowledge of the adjoint solution rz, rs of the continuous problem (4.3).In practice, this will not be available and we are therefore required to compute a numerical approximation rẑ ĥ, k, rĥ , ks to rz, rs, based on employing a mesh of granularity ĥ with polynomials of degree k.Notice that the Galerkin orthogonality property necessitates that rẑ ĥ, k, rĥ , ks R V 0 h,k ˆQh,k , otherwise the error estimate will be identically zero.Therefore, we adopt the approach of computing a finite element approximation to rz, rs with T ĥ " T h and k " k `1, i.e., rẑ ĥ, k, rĥ , ks " rẑ h,k`1 , , rh,k`1 s P V 0 h,k`1 ˆQh,k`1 , cf. [4].Replacing rz, rs with rẑ h,k`1 , rh,k`1 s in (4.4) then yields the (approximate) error estimate Rewriting the above error estimate as a summation of local error indicators η κ , over the elements κ of the mesh T h , we arrive at the following (approximate) a posteriori bound.
Proposition 4.1.Under the foregoing notation, the following (approximate) a posteriori bound holds: where η κ " R `ru h,k , H h,k s, rẑ h,k`1 ´zh,k , rh,k`1 ´rh,k s ˘|κ is the weighted residual evaluated on a single element κ P T h .

4.2.
Non-Linear Functionals.Suppose now that the quantity of interest corresponds to a non-linear functional Jp¨q, that is Fréchet differentiable with derivative at rv, qs denoted by J 1 " rv, qs ‰ p¨q.We define the mean value linearization of J by J" ru, Hs ‰" In this case, the formal adjoint problem corresponding to the quantity of interest J is then given by : find rz, rs P H 0,N pdiv, Ωq ˆL2 pΩq such that A `rv, qs, rz, rs ˘" J" ru, Hs ‰" ru h,k , H h,k s ‰`r v, qs ˘@rv, qs P H 0,N pdiv, Ωq ˆL2 pΩq.(4.6)On the basis of this adjoint problem, together with the mean value linearization of J, cf.above, we deduce the following error representation formula: In practice, however, as the solution to the weak formulation given in equation (3.3) is unknown, we make the assumption that the underlying linearization error is small and employ the (approximate) adjoint problem: find rz, rs P H 0,N pdiv, Ωq L2 pΩq such that A `rv, qs, rz, rs ˘" J 1 " ru h,k , H h,k s ‰`r v, qs ˘@rv, qs P H 0,N pdiv, Ωq ˆL2 pΩq; (4.7) for simplicity of notation, we use rz, rs to denote the solution to both adjoint problems (4.6) and (4.7), though we stress that these are clearly not identical, in general.Writing rẑ h,k`1 , rh,k`1 s P V 0 h,k`1 ˆQh,k`1 to denote the numerical approximation to (4.7) computed on the mesh T h , with polynomials of degree k `1, proceeding as above, we deduce that where η κ is defined in analogous manner to the case when J is a linear functional.

Travel Time Functional.
For notational convenience, we rewrite the travel time functional defined in (2.4) in the following equivalent manner: As noted in Remark 2.1, the smoothness of the travel time functional is unclear; the key technical difficulty stems from integrating the quantity of interest on a path which is unknown a priori.Given the difficulty in demonstrating Fréchet differentiability of the travel time functional in some neighbourhood of the solution ru, Hs, we adopt a different approach to that presented in section 4.2 for Fréchet differentiable functionals.With this in mind, we exploit a numerical approximation of the Gâteaux derivative in the definition of the corresponding adjoint problem.The idea of numerically approximating these derivatives has been used in a range of applications, such as inexact Newton methods [14,18], the numerical solution of stiff systems of ODEs [13] and minimization techniques [27].
Employing a generalized Taylor expansion and neglecting terms that are Op q, we may approximate the Gâteaux derivative of J , evaluated at ru, Hs, in the direction rv, qs, based on employing the one sided difference formula where 0 ă ! 1.With this approximation, we formally consider the following adjoint problem corresponding to the travel time functional: find rz, rs P H 0,N pdiv, ΩqˆL 2 pΩq such that The adjoint solution may be interpreted as a generalized Green's function relating the partial differential equation given in (2.1) and (the linearization of) the travel time functional in equation (4.3).In the current setting, the linearization of the travel time functional is approximated by (4.8).However, we stress that this approximation leads to the formulation of an approximate adjoint problem (4.9), whose right-hand side functional, i.e., J 1 " ru h,k , H h,k s ‰`r v, qs ˘is nonlinear with respect to the test function rv, qs; thereby, we may not apply a standard Galerkin finite element method to approximate the adjoint solution rz, rs defined by (4.9).To overcome this issue, we introduce a further approximate discrete adjoint problem which is defined based on employing a linear functional approximation of J 1 " ru h,k , H h,k s ‰`r v, qs ˘.To this end, we define the unique linear functional J : V 0 h,k`1 ˆQh,k`1 Ñ R which precisely coincides with the nonlinear functional J 1 when evaluated/sampled at each of the basis functions which form the spanning set for the finite element space V 0 h,k`1 ˆQh,k`1 .Thereby, we proceed by defining the approximate discrete adjoint problem by: find 10) The study of the well-posedness of the approximate adjoint problem (4.9) and its approximation by the discrete (approximate) counterpart (4.10) are beyond the scope of this article; for the purposes of the proceeding analysis, we assume that the former is indeed the case and that (4.10) provides a suitably accurate approximation of the analytical solution of (4.9).In particular, we demonstrate numerically that the resulting (approximate) error representation formula, cf.(4.11) below, leads to highly efficient estimation of the discretization error in the travel time functional on optimized computational meshes.
Under the assumption that the error committed in the approximation of the Gâteaux derivative of J is small compared to the discretization error in the computation of the primal solution, measured in terms of the travel time functional, in practice we compute the following approximate error representation formula where rẑ h,k`1 , rh,k`1 s satisfies (4.10); as before, we exploit the same mesh T h employed for the computation of the primal approximation ru h,k , H h,k s, using the enriched polynomial degree k `1.The local (weighted) elementwise error indicators η κ are defined in a similar manner as in the previous two subsections.
For a heuristic choice of , we follow the arguments presented in [22].When the functional is evaluated on a computer with finite precision, the computed value of the functional is only evaluated to a relative accuracy of m ; for double precision IEEE arithmetic, m « 10 ´16 .Writing J ‹ to denote the numerical evaluation of J , then for some δ " Op m q, we note that where rv, qs P H 0,N pdiv, Ωq ˆL2 pΩq.Thereby, the numerator in the approximation of J 1 , cf. the right-hand side of (4.8), may be computed to a relative accuracy of Op m q.
Hence, the approximation of the derivative of the travel time functional involves two error terms: one corresponding to the round off error in the evaluation of J , and the other corresponding to the truncation of the Taylor series expansion of the derivative.Taking both of these error contributions into account, we deduce that To control the error in the evaluation of J 1 , and thereby also J , should be chosen in order to balance the two error terms arising in (4.13).With this in mind, we require that « ?m ; for double precision arithmetic, a suitable choice is to select « 10 ´8.This implies that this approach will not be able to provide a good representation of the discretization error to full precision and will become polluted by either truncation or round off error once the discretization error is of the same order.However, for the engineering applications of interest in this article, this level of precision is more than sufficient.

Numerical Experiments.
In this section we present a series of numerical experiments to demonstrate the quality of the computed error representation formula (4.11) within an automatic adaptive mesh refinement algorithm.Here, the elements are marked for refinement/derefinement based on the size of the local error indicators |η κ | using the fixed fraction refinement strategy, with refinement and derefinement fractions REF% and DEREF%, respectively.The computations presented in this section have been undertaken based on employing the AptoFEM finite element package [2].Here, the primal finite element solutions are computed using RT 0 elements for the Darcy velocity, with a piecewise-constant approximation of the hydraulic head, i.e., k " 0, while the corresponding adjoint finite element solutions are computed with k " 1; i.e., RT 1 elements are employed for the adjoint Darcy velocity, with a discontinuous piecewise linear approximation to the adjoint hydraulic head.We point out that with the restriction to piecewise-constant pressures, we can link the mixed finite element method employed here to commonly used finite volume methods [15,49,48].
Finally, we discuss the numerical evaluation of the travel time functional J rv h,k , q h,k s, for a given rv h,k , q h,k s P V 0 h,k ˆQh,k , k ě 0. Firstly, we note that, for divergence-free flows, the exploitation of RT 0 elements, i.e., k " 0, implies that the numerical approximation of the Darcy velocity is constant elementwise and, as such, the computation of the path Ppv h,0 q and travel time are trivial.In examples in which the velocity is not divergence-free, the method proposed in [37] may be employed for the computation of the path involved in the evaluation of the travel time functional.For the assembly of the discrete adjoint solution, we must compute the approximation to the derivative of J given in (4.8).This requires the computation of the path and travel time for each degree of freedom corresponding to an element that intersects the primal path, i.e., Ppu h,0 q.In this case, it is necessary to compute the path and travel time functional for a Darcy velocity that lies in RT 1 and, as such, the method employed for RT 0 Darcy velocities is no longer applicable.For the divergence-free case it should be possible to use the streamline method proposed in [36] to compute an approximation to the path.For simplicity, however, we employ a forward Euler discretization with sufficiently small time steps that any discretization error is negligible and may be ignored.It should be noted, however, that this will occur only for the elements corresponding to the degree of freedom for which the numerical approximation of the derivative is being evaluated; this will be at most two elements per path, with the remainder of the path being either pre-computed or computed as required using the method in [37].
Example 1.In this first example, we consider a problem that has a known analytical solution, which is sufficiently simple that we may compute an exact value for the travel time functional for a given release point x 0 .To this end, we let Ω " p0, 1q 2 , Γ D " BΩ, Γ N " H, K " I, f " ´6px `1q, and g D " px `1q 3 `y `3; thereby, the analytical solution to the system of PDEs (2.1) is given by u " ´ˆ3px `1q 2 1 ˙and H " px `1q 3 `y `3. (5.1) Setting φ " 1 and x 0 " p0.9, 0.2q, the exact value of the travel time functional may be evaluated as follows.The x-component of the path satisfies the ordinary differential equation 2), (2.3).This implies that Assuming that the path intersects the boundary of Ω at x " 0, gives t " 1 {3 p1 ´1{p0.9`1qq" 3 {19; i.e., J `ru, Hs ˘" 3 {19.Repeating this process for y demonstrates that X p 3 {19q " p0, 23 {95q, confirming our assumption that the path intersects the boundary of Ω at x " 0. We point out that this problem is not of physical relevance, but serves as a useful testcase to demonstrate the practical performance of the proposed a posteriori error indicator on adaptively refined computational meshes.With this in mind, we define the effectivity index θ, as follows: θ " ř κ η κ J `ru, Hs ˘´J `ru h,0 , H h,0 s ˘. (5.2) For this example, we set " 10 ´8, and refinement and derefinement percentages to REF " 20% and DEREF " 10%, respectively.In Figure 5.1 we show the finite element solutions ru h,0 , H h,0 s and rz h,1 , r h,1 s computed on the initial mesh, denoted by T p1q h .From the plots of the adjoint solution rz h,1 , r h,1 s, we can observe that there is an apparent discontinuity in r h,1 along the path Ppu h,0 q.The adjoint velocity z h,1 is zero almost everywhere in the domain and only non-zero in the vicinity of Ppu h,0 q.Interpreting the adjoint solution as a generalized Green's function for the travel time functional and the system of PDEs given in equation (2.1), we observe that the computed numerical approximation corresponds to a δ´function type source, or sink, along Ppu h,0 q.This is, in a sense, qualitatively similar in character to the adjoint solutions computed for first-order hyperbolic conservation laws, when the functional of interest is a point evaluation of the primal solution, cf.[30,32], for example.However, in this setting, we do observe some growth in the height of the δ´type adjoint solution along the path of interest, as we move from the release point to the boundary of the computational domain Ω.
In Figure 5.2 we plot the third mesh T p3q h , and the corresponding path, obtained using the proposed adaptive mesh refinement strategy.Here, we observe that the adaptive algorithm has refined elements in the region surrounding the path, with more refinement taking place around the release point x 0 , and the point at which the path intersects BΩ.The induced mesh refinement is due to the presence of the weighting terms involving the difference between the (approximated) adjoint solution Finally, in Table 5.1 we demonstrate the performance of the proposed adaptive strategy; here, we show the number of degrees of freedom in underlying finite element space V 0 h,0 ˆQh,0 , the true error in the functional J `ru, Hs ˘´J `ru h,0 , H h,0 s ˘, the computed error representation formula ř κPT h η κ , and the effectivity index θ.Here, we see that the quality of the computed error representation formula is extremely good, in the sense that the effectivity indices are very close to unity on all of the meshes employed.
Example 2. In this example, we consider a more complex testcase in order to demonstrate the applicability of the proposed error estimation techniques for the travel time functional for a more realistic groundwater flow example.To this end, we model the flow in a two-dimensional domain consisting of six rock strata, each with greatly differing hydrogeological properties.The geological units considered are based on those found at the Sellafield site in the UK; for details, we refer to [41,44].Note that, since we only intend to demonstrate the applicability of the proposed a posteriori estimate, we consider a greatly simplified geometry; in particular, this does not represent the geometry considered in any hydrogeological assessment.Moreover, we neglect any faults or other complex geological or topographical features.With this in mind, the results presented here are not of direct relevance to the results presented in [41,44].
Here, we let the computational domain Ω be as shown in Figure 5.3.We partition this domain into a series of sub-domains Ω i , i " 1, . . ., 6, that correspond to each of the geological units shown, and are enumerated from the top of the domain to the bottom, i.e., Ω 1 corresponds to Calder sandstone, and so on.In order to prescribe boundary conditions for this problem we make the following modelling assumptions: i) The rock below the Borrowdale Volcanic Group (BVG) stratum is of a much lower permeability than the other geological units; ii) There is a flow divide on both the left and right edges of the domain; and iii) The atmospheric pressure along the top of the domain p atm , is 1.013 ˆ10 5 P a.
Given these assumptions, we set Γ D be the top of the boundary and Γ N be the remainder.We define the Dirichlet boundary data g D by where z is the distance above the vertical datum shown at 0 in Figure 5.3, ρ is the density of water, and g is acceleration due to gravity.In addition, the right-hand side forcing function f is set equal to zero.The data for permeability k and porosity φ are taken from [44] and are summa- rized in Table 5.2, where the permeability has units of m 2 .The data corresponds to the mean mp¨q and standard deviation σp¨q of log normally distributed random variables that represent the upscaled permeabilities and porosities of the geological units shown in Figure 5.3.As we are considering deterministic flow in this article, we simply take the mean value as the permeability.This yields piecewise-constant     permeability and porosity throughout Ω.We make the simplifying assumption that the permeability is isotropic.The hydraulic conductivity K is then related to the permeability k by where µ is the kinematic viscosity of the groundwater.The values for these quantities used in the numerical experiments is contained in Table 5.3.
Whilst the hydraulic conductivity is piecewise constant, it is highly discontinuous at the boundaries of the geological units, varying across all geological units by three orders of magnitude.Throughout this section, we ensure that the computational mesh T h respects the boundaries of these geological units.Additionally, we choose appropriate units such that the resultant travel time is Op1q, enabling the choice of " 10 ´8 in the computation of the approximate adjoint solution.Finally, we choose refinement and derefinement percentages to be REF " 25% and DEREF " 15%, respectively, and select x 0 " p5.0, ´2.0q.
Figure 5.4 shows the finite element solutions ru h,0 , H h,0 s and rz h,1 , r h,1 s computed on the initial mesh T p1q h , yielding a finite element space V 0 h,k ˆQh,k of dimension 14641 degrees of freedom.The geometry of the rock strata is superimposed to demonstrate the differences in the solutions between the different regions.Considering the adjoint solution rz h,1 , r h,1 s, we observe that there is an apparent discontinuity in r h,1 along the path Ppu h,0 q, as was the case in Example 1, in the strata of BVG, Carboniferous Limestone and Brockram Breccia.Again, the adjoint velocity z h,1 is zero in the majority of the domain and only non-zero close to Ppu h,0 q, cf.Example 1.However, in this example, there are interesting features in the adjoint velocity close to the release point x 0 and close to Ppu h,0 q in the stratum of BVG; here, the velocity along the path is in the opposite direction to the primal velocity, however, the velocity    surrounding the path is in the same direction.Figure 5.5 shows the first four meshes generated by the proposed adaptive algorithm, denoted by T p1q h

´T p4q
h , and the paths computed based on exploiting the primal solution computed on this mesh.Here, we can clearly observe that the adaptive algo- rithm has refined elements in the region surrounding the path, with more refinement taking place around the release point x 0 , as well as in the strata of Carboniferous Limestone, Brockram Breccia, and St. Bees evaporites.From the plot of the path we can see there are rapid changes in direction of the path and, as such, it is reasonable to expect that these could be a significant potential cause of discretization error.In Table 5.4 we present the performance of the proposed adaptive algorithm for the estimation of the travel time functional.Here, the effectivity indices are computed based on evaluating the travel time functional on a very fine adaptively refined computational mesh, which yields J `ru, Hs ˘« 0.49 ˆ10 5 years.As for the previous example, we again observe that the effectivity indices computed on all of the meshes generated by our adaptive refinement strategy are very close to one, which indicates that reliable and efficient estimation of the error in the computed travel time functional has been attained for this physically motivated example.
To test the robustness of the proposed goal-oriented a posteriori error indicator, we now consider the application of the underlying adaptive algorithm starting from a very coarse mesh.To this end, we construct a mesh so that the resulting finite element space V 0 h,k ˆQh,k consists of only 1253 degrees of freedom, cf. Figure 5.6(a).The final three meshes generated by the our adaptive algorithm, together with the computed path, are also depicted in Figure 5.6.Once again, we observe that the mesh is refined in the vicinity of the computed path.The performance of the approximate error representation formula computed on this sequence of meshes is given in Table 5.5.Here, we observe that while there is some degradation in the quality of the computed effectivity indices, compared with those presented in Table 5.4, which were computed based on starting from a finer initial mesh, the computed error representation formula still provides a highly accurate approximation to the error in the computed travel time functional, even on much coarser finite element meshes.
Finally, in Figure 5.7 we compare the performance of the goal-oriented a posteriori error indicator proposed in this article with a standard energy norm error estimator; for the latter, we refer to [9].Here we clearly observe the superiority of the goaloriented a posteriori error indicator; on the final mesh, the true error in the travel time functional is over two orders of magnitude smaller than J `ru, Hs ˘´J `ru h,0 , H h,0 s computed on the sequence of meshes produced using the energy indicator.Moreover, we note that the true error in the computed travel time functional does not converge to zero when the latter refinement indicator is employed; indeed, the energy indicator concentrates the mesh in the regions of Brockram Breccia and St. Bees Evaporates, where the inverse of the hydraulic conductivity is smallest, cf. Figure 5.8.
Example 3. In this final example, we consider a more challenging version of the problem considered in Example 2. To this end, we consider the same geometry and boundary conditions as set out in the previous example, but introduce variability in k within each stratum by modelling the permeability as a realization of a spatially correlated log-normal random field, cf.[35,47].On the initial finite element mesh T p1q h , we compute a realization of k, based on evaluating the Karhunen-Loève (KL) decomposition of the random permeability field in each stratum, cf.[40].To this end, for a given stratum l, 1 ď l ď 6, let T p1q h,l denote the finite element submesh consisting of elements from T p1q h which lie within the current stratum.For ease of presentation, we assume that the elements within T p1q h,l are numbered consecutively from 1 to M l .Writing x i , i " 1, . . ., M l , to denote the centroids of each element in the submesh T p1q h,l , consider the exponential covariance function in stratum l given by where λ " 1 denotes the correlation length and σ l is the standard deviation of the logarithm of the permeability of the l th stratum.We decompose the M l ˆMl covariance matrix for the element centroids, pC l q ij " C l px i , x j q, i, j " 1, . . ., M l , into where Φ l is the matrix of eigenvectors of C l and Λ l is the (ordered) matrix of eigenvalues.Then in the l th stratum the value of k " k l in the i th element κ i P T p1q h,l is given by k l | κi " exppZ i q, (5.7)  i " 1, . . ., M l , where Z i is defined by where ξ j are iid random variables with ξ j " N p0, 1q, j " 1, . . ., M l , and m l is the mean of the logarithm of the permeability of the l th stratum.The computed realization, evaluated on the initial finite element mesh, is depicted in Figure 5.9; here, we observe that the minimum and maximum of the permeability vary by 6 orders of magnitude.The performance of the proposed goal-oriented adaptive refinement algorithm is presented in Table 5.6.As for the previous examples, we again observe that the quality of the computed error representation formula is extremely good, even for this more complex and demanding problem; indeed, the effectivity indices are very close to unity on all of the meshes generated by our adaptive algorithm.Figure 5.10 shows the first four meshes generated by the proposed adaptive algorithm, denoted by T p1q h ´T p4q h , and the paths computed based on exploiting the primal solution computed on this mesh.Here, we can clearly observe that the adaptive algorithm has refined elements in the region surrounding the path, with more refinement taking place around the release point x 0 .In particular, we see that the trajectory of the path has changed dramatically, compared with the path computed in Example 2.Moreover, here we observe some small 'wiggles' in the path generated by the variation in the permeability field.

Conclusions.
In this article we have introduced a goal-oriented a posteriori error estimate for the travel time functional for a non-sorbing solute transported by groundwater from some release point deep underground to the biosphere.On the basis     on this error bound, we have designed and implemented the corresponding adaptive finite element mesh refinement strategy.This general approach has been validated for both a simple problem with a known analytical solution, together with a more realistic example based on the hydrogeology of the Sellafield site in the UK, cf.[41,44].
There are several natural extensions to the work considered in this article; indeed, the robustness of the error estimate under more demanding conditions is crucial for more realistic applications.This verification could include, for instance, testing the error estimate in a domain that includes a fracture network or other complex topographical features or using a standard permeability-porosity data set.Application of this method to random porous media will be considered in our forthcoming article [19]; of course, the consideration of three-dimensional problems is also of significant interest.Finally, there are also important questions related to the analysis of the underlying adjoint problem; indeed, the well-posedness of both the continuous weak formulation and its discrete counterpart needs to be addressed, together with the regularity of the adjoint solution.

Fig. 5 . 4 :
Fig. 5.4: Example 2. The primal and adjoint data obtained on the initial mesh with the of the rock strata superimposed.
The initial mesh, T p1q h , and resultant path.
The second mesh, T p2q h , and resultant path.
The third mesh, T p3q h , and resultant path.
The fourth mesh, T p4q h , and resultant path.

Fig. 5 . 5 :
Fig. 5.5: Example 2. The initial mesh and the meshes obtained after the first three applications of the adaptive algorithm together with the path obtained using the primal solution corresponding to the mesh.
The fourth mesh, T p4q h , and resultant path.
The fifth mesh, T p5q h , and resultant path.
The sixth mesh, T p6q h , and resultant path.

Fig. 5 . 6 :
Fig. 5.6: Example 2.The initial mesh and the fourth, fifth, and sixth meshes generated after the application of the adaptive algorithm together with the path obtained using the primal solution corresponding to the mesh.
The initial mesh, T p1q h , and resultant path.
The second mesh, T p2q h , and resultant path.
The third mesh, T p3q h , and resultant path.
The fourth mesh, T p4q h , and resultant path.

Fig. 5 . 10 :
Fig. 5.10: Example 3. The initial mesh and the meshes obtained after the first three applications of the adaptive algorithm together with the path obtained using the primal solution corresponding to the mesh.

Table 5 . 1 :
DOFS J `ru, Hs ˘´J `ru h,0 , H h,0 s ˘řκPT h η κ θ Example 1. Error data obtained from the adaptive algorithm.rẑ h,k`1 , rh,k`1 s and rz h,k , r h,k s, which multiply the computable residual terms involving the numerical solution ru h,k , H h,k s in the definition of the local error indicators |η κ |.These weights represent the sensitivity of the relevant error quantity with respect to variations of the local element residuals.

Table 5 . 3 :
Example 2.Values for physical constants that are used to compute the hydraulic conductivity.

Table 5 . 5 :
Example 2.Error data obtained from the adaptive algorithm, based on starting from a coarse initial mesh.