Postprocessing of standard finite element velocity fields for accurate particle tracking applied to groundwater flow

Particle tracking is a computationally advantageous and fast scheme to determine travel times and trajectories in subsurface hydrology. Accurate particle tracking requires element-wise mass-conservative, conforming velocity fields. This condition is not fulfilled by the standard linear Galerkin finite element method (FEM). We present a projection, which maps a non-conforming, element-wise given velocity field, computed on triangles and tetrahedra, onto a conforming velocity field in lowest-order Raviart-Thomas-Nédélec (RTN0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {RTN}_{0}$\end{document}) space, which meets the requirements of accurate particle tracking. The projection is based on minimizing the difference in the hydraulic gradients at the element centroids between the standard FEM solution and the hydraulic gradients consistent with the RTN0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {RTN}_{0}$\end{document} velocity field imposing element-wise mass conservation. Using the conforming velocity field in RTN0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {RTN}_{0}$\end{document} space on triangles and tetrahedra, we present semi-analytical particle tracking methods for divergent and non-divergent flow. We compare the results with those obtained by a cell-centered finite volume method defined for the same elements, and a test case considering hydraulic anisotropy to an analytical solution. The velocity fields and associated particle trajectories based on the projection of the standard FEM solution are comparable to those resulting from the finite volume method, but the projected fields are smoother within zones of piecewise uniform hydraulic conductivity. While the RTN0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {RTN}_{0}$\end{document}-projected standard FEM solution is thus more accurate, the computational costs of the cell-centered finite volume approach are considerably smaller.


Introduction
Groundwater flow is commonly simulated by substituting Darcy's law into the continuity equation, resulting in an elliptic (steady-state problem) or parabolic (transient problem), second-order differential equation of the hydraulic head [4,29,32]. Conservative solute transport is traditionally described by the advection-dispersion equation (ADE) employing the velocity field originating from the solution Olaf A. Cirpka olaf.cirpka@uni-tuebingen.de Philipp Selzer philipp.selzer@uni-tuebingen.de quadrilaterals and deformed hexahedra, split into tetrahedra, for non-divergent flow by Cordes and Kinzelbach [11]. Pollock's [34] method has been classified as semi-analytical because it uses an analytical solution for particle trajectories based on a numerical approximation of the velocity field. It computes the exit point from a particle's entry point on an element boundary by evaluating the fluxes over the volume boundaries, assuming that the normal component of the velocity vector on each volume-face is constant and each velocity component varies linearly in its own direction within the cell [29]. This assumption about the velocity field corresponds to the lowest-order admissible velocity approximation [23,29]. Such a velocity field is in the lowest-order Raviart-Thomas-Nédélec (RT N 0 ) space, which is also used for the velocity approximation of lowestorder mixed (hybrid) finite element methods [3, 29-31, 37, 38]. Higher-order schemes are possible but seldom used in the simulation of groundwater flow and solute transport [7,8,23,29].
By construction, finite volume methods lead to conforming velocity fields. Such a field is characterized by a continuous normal component of the velocity field on element boundaries and local mass conservation in each element. While continuous Galerkin finite element methods approximate the unknown head field as continuous functions, they do not yield a conforming velocity field [32,37]. If a particle is tracked on the basis of an element-wise approximated velocity field originating from P 1 Galerkin FEM, the non-conforming property of the velocity field leads to severe numerical artifacts, possibly including particle stagnation because the approximated normal velocity component points into opposite directions at the two sites of an element interface [37]. Therefore, the consistent application of particle tracking methods to standard FEM flow fields remains a problem.
Several authors have emphasized that the P 1 Galerkin FEM is locally mass conservative if appropriate grids and control volumes, which differ from the elements, or nodal fluxes are considered [15,21,37]. Especially, locally mass-conservative sub-control volumes can be assigned to the element-wise approximation of the Darcy velocity field obtained by P 1 Galerkin FEM for some sub-controlvolumes on triangular elements carrying heterogeneous material coefficients [37]. According to the latter authors, such regions can be obtained by Voronoi tessellation or by patches around the nodes bounded by direct connections between the midpoints of edges of triangular elements sharing the respective node. They call these control volumes "internal A-cells." A delineation of control volumes in analogy to vertex-centered finite volumes for triangular grids was proposed by Durlofsky [15]. In three dimensions, however, comparable control volumes are not definable with similar ease if hydraulic conductivity varies in space [36,37]. Also, a conforming velocity interpolation within polygonal control volumes with a large number of faces, which may result from the reconstruction of Durlofsky [15] or from vertex-centered Finite Volumes, can be nontrivial. That is, while local mass conservation is a necessary condition for reliable particle tracking, it may not be sufficient. In any case, considering the elements of a P 1 Galerkin FEM grid itself, fluid mass is not conserved, and the normal component of the velocity is discontinuous, experiencing a jump on element interfaces, while the primary unknown is continuous [37,41].
Cordes and Kinzelbach [11] introduced a scheme to reconstruct a conforming velocity field on linear triangular and bilinear quadrilateral Galerkin finite elements. It is based on the mass conservation property of the internal A-cells [11,37]. Considering the total fluxes on the boundaries of the internal A-cell, a system of equations can be set up for the total fluxes on all its inner boundaries [11]. This system has N equations with N variables, in which N is the number of inner edges of the internal Acell. However, the circular structure causes the system of equations to be underdetermined. For regularization, Cordes and Kinzelbach [11] set the constraint that the rotation of the hydraulic gradient within the elements must be zero. A direct extension of the method of Cordes and Kinzelbach [11] to three dimensions is impossible, as there are 1.5 times more unknowns than equations related to internal nodes for tetrahedra, resulting in a highly underdetermined system of equations.
Larson and Niklasson [27] presented an approach to construct conforming velocity fields by postprocessing of P 1 Galerkin FEM solutions that was originally defined for the case of a constant isotropic hydraulic conductivity. Sun and Wheeler [43] introduced a similar, element-wise approach of flux correction to obtain local mass conservation and continuous normal components of fluxes on element boundaries for non-conforming velocity approximations, originating from continuous Galerkin finite elements [35,41]. The approach of Larson and Niklasson [27] and Sun and Wheeler [43] is based on a discontinuous enrichment of the velocity field by applying a piecewise constant correction term which is added to the non-conservative, element-wise velocity approximation of the FEM on the element interfaces [33]. An adapted formulation applied to the Richards equation for variably saturated flow is given by Kees et al. [24] and further employed and investigated by [41].
An advantage of the flux correction proposed by Larson and Niklasson [27] and Sun and Wheeler [43] is that it can be formulated as a local variant, which is parallelizable [27,41]. However, numerical experiments conducted by Schiavazzi [40] indicate that numerical errors are nonnegligible for this scheme, if hydraulic conductivity strongly varies in space even in the isotropic case. Schiavazzi [40] shows that the absolute values of element-wise constant Darcy velocities can be orders of magnitude difference in the postprocessed solution compared to the initial estimate of a velocity space based on a P 1 Galerkin FEM. Considering the particle tracks presented by Schiavazzi [40], it appears that also some artificial rotation is numerically introduced in the postprocessed solution. This might also be the reason why neither most of the particle tracks evaluated by Schiavazzi [40] cross significantly zones of low permeability nor concentrations evaluated by Scudeler et al. [41] enter zones of lower permeability to a significant extent, although Schiavazzi [40] reports that the absolute magnitude of the elementwise postprocessed velocities can be orders of magnitude higher than those estimated on the basis of a P 1 Galerkin FEM. Furthermore, it has to be noted that, depending on the actual computational scheme, numerically induced rotation is known to become a problem in the computation of the RT N 0 mixed (hybrid) finite element method [20].
Recently Odsaeter et al. [33] presented an adapted scheme, in which the residual between an elementwise given non-conforming velocity approximation and a conforming velocity solution is minimized with respect to a weighted L 2 -norm imposing local mass conservation. Odsaeter et al. [33] introduce weighting factors, which are chosen such that the normal component of the gradient rather than the Darcy velocity is considered on element boundaries, leading to very good results that apparently avoid the problems of the flux corrections applied and evaluated by Schiavazzi [40] and Scudeler et al. [41].
In this paper, we present a novel formulation of an RT N 0 projection employing concepts similar to those described by Odsaeter et al. [33]. Our RT N 0 projector maps a non-conforming, element-wise approximated Darcy velocity field of an FEM solution, representing the nonmass-conservative flow, onto a conforming velocity field in lowest-order Raviart-Thomas-Nédélec space. The target velocity space ensures continuity of the normal velocity component on element boundaries, zero divergence of the element-wise velocity field for non-divergent flow, and element-wise mass conservation. Implemented as a postprocessing code, the RT N 0 projector is formulated such that it can be coupled to any FEM code, but it comes at the costs that a new global system of equations has to be solved, which is of a size comparable to that of mixed finite elements. The key of the approach is that we minimize differences in the hydraulic gradient at the element centroids between the P 1 Galerkin FEM solution of hydraulic heads and the target RT N 0 velocity field, subject to the constraint that the RT N 0 velocity field is element-wise mass conservative. We consider the residuals of the hydraulic gradient rather than that of the Darcy velocity because the hydraulic gradient varies much less than the Darcy velocity if the hydraulic conductivity field is heterogeneous. Moreover, such a minimization procedure should guarantee that the postprocessed velocity field has the same order of convergence as the original velocity approximation. In this regard, our postprocessing scheme is equivalent to the projection presented by Odsaeter et al. [33], who give a thorough analysis of the convergence behavior applicable also to our scheme. We found that a projection, in which residuals of the Darcy velocity were minimized, comparable to flux corrections formerly applied by Schiavazzi [40] and Scudeler et al. [41], introduced severe numerical artifacts because meeting the Darcy velocities in low-conductivity regions was less important in the optimization procedure than meeting those in high-conductivity regions, resulting in velocities of too high magnitude, and erroneous rotation and even reversal of the flow direction in zones of low hydraulic conductivity. Our RT N 0 projector is reasonably applicable to any finite element method yielding non-conforming element-wise velocity fields. It can, in principle, be applied to various element types also in higher dimensions. The theory is given for simplices in two and three spatial dimensions, i.e., triangles and tetrahedra. For demonstration purposes, we evaluate the numerical results on triangular and tetrahedral grids in detail. For comparison purposes, we define a cell-centered finite volume method on simplicial centroids against which the results of the RT N 0 projector are evaluated. On this basis, we present semi-analytical particle tracking methods for divergent and non-divergent flow, similar to those of Pollock [34].

Governing equation
We consider a bounded model domain ⊂ R d with dimensionality d ∈ {2, 3}. The domain boundary, ∂ = , is subdivided into a Dirichlet boundary, D , and a Neumann boundary, N , such that ∂ = D ∪ N and D ∩ N = ∅.
We consider the elliptic steady-state groundwater flow equation: subject to the following boundary conditions: is a volumetric source/sink term per unit volume, q N is a specified volumetric flux density, and n is the unit normal vector pointing outwards. If the material properties are locally isotropic, the hydraulic conductivity tensor reduces to K = kI, where k ∈ R [LT −1 ] is the directional independent scalar hydraulic conductivity and I is the identity matrix of order d. Furthermore, the specific discharge, or Darcy velocity, q ∈ R d [LT −1 ], is defined by Darcy's law: and related to the average linear velocity is the porosity assumed constant in the present analysis, even though the extension to element-wise constant porosities is straightforward.

Spatial discretization of the domain
Let T be a topological triangulation of into triangles (d = 2), or tetrahedra (d = 3). Furthermore, we restrict ourselves to conforming grids such that boundaries of neighboring elements perfectly fit. Elements within the triangulation are denoted E i ∈ T , for i = 1, 2, ..., N, in which N is the total number of elements. E I is the set of inner elements, E D is the set of elements sharing a Dirichlet boundary, and E N is the set of elements sharing a Neumann boundary.
The set of nodes of T is denoted N , while N E is the set of nodes of the element E. Note that all sets of faces are pairwise disjoint such that any face F is an element of only one of the sets, F I , F D , or F N . Furthermore, let F E be the set of faces of the element E, such that F ∈ F E i is an individual face of the element E i . On every face, there is a unit normal vector, ν F , following a sign convention. Furthermore, on every face. there is an outer unit normal vector, n E , pointing out of E.

Further notation
We denote measures by | · |, in particular |E| denotes the area (d = 2), or the volume (d = 3) of an element, E ∈ T , and |F | denotes the length (d = 2), or the area (d = 3) of a face, F ∈ F. The Euclidean norm of a vector is denoted by · 2 , while ·,· is the standard inner product. Furthermore, [·] F denotes the jump of a function on a face.

Piecewise polynomial spaces
For clarification purposes, we briefly introduce some relevant, fundamental function spaces skipping their subspaces refined by boundary conditions [3,32,33]. In the following sections, we will only refer to the respective function spaces, if we believe this is absolutely necessary for proper understanding of computations. Let P 1 be the space of continuous piecewise linear polynomials such that where C is the space of continuous functions and Q d 1 denotes the tensor product of linear polynomial spaces in each spatial direction [33]. The space P 1 (T ) is defined such that it contains the approximated solution of the hydraulic head h computed by the standard Galerkin FEM. Therefore, it follows that Galerkin solution. That is, the hydraulic head is conforming and continuous on elemental boundaries, while the normal component of the element-wise evaluated specific discharge can be discontinuous, experiencing a jump on element boundaries, leading to non-conforming specific discharges.
In lowest-order Raviart-Thomas-Nédélec space for simplices [3,30,31,38,39], a vector-valued function in R d , in our case the specific discharge, is considered element-wise such that in which L 2 ( , R d ) is the usual Lebesgue space of squareintegrable functions in R d [3]. Considering the definition of the RT N 0 space in equation 5, it is obvious that not only the normal component of a vector-valued function is continuous on the faces, but also its divergence is elementwise constant and linear dependent on b, i.e., the divergence is element-wise zero for non-divergent flow. In this study, we intend to project an element-wise approximated specific discharge field originating in our application from a P 1 Galerkin FEM onto a conforming velocity field in RT N 0 space. This can be expressed by the following map:

P 1 Galerkin FEM
The standard Galerkin FEM is a common numerical method to solve the groundwater flow equation (Eq. 1 [4,22,25]). P 1 Galerkin FEM yields a continuous approximation of the hydraulic head field and can easily handle unstructured grids and anisotropy. Considering the space P 1 (T ), the linear shape functions χ i corresponding to the nodes i can be defined element wise such that they sum up to unity at all points in all elements: in which χ i are the shape functions according to the space P 1 (T ) equaling 1 on node i and 0 on any other node j = i. We replace the hydraulic head field h(x) and the source/sink term field f (x) by the approximate, piecewise linear functionsĥ andf , respectively, such that in which h i is the hydraulic head value at node i. Applying the weak form of the groundwater flow equation (Eq. 1) with the weighting functions being identical to the shape functions, the standard Galerkin FEM discretization reads element wise as in which χ E are the elemental shape functions, K E is the element-wise given hydraulic conductivity tensor,ĥ are the nodal hydraulic head values of element E, andf are the nodal sources/sink strengths. A nodal formulation of Eq. 9 in analogy to Forsyth [17] is also possible. Based onĥ, a unique, element-wise constant gradient can be computed on simplices, leading to an element-wise constant, nonconforming specific discharge vector. This unique gradient can be obtained by evaluating the gradient of the shape functions multiplied with the vector of nodal hydraulic head values. Then, the element-wise approximated specific discharge vector is approximated by

Cell-centered FVM on simplices
In a cell-centered finite volume method (FVM) applied to the groundwater flow equation (Eq. 1), the conservation of fluid volume is enforced for each element E of a given spatial subdivision, or topological triangulation T : Then, the divergence theorem is employed such that the normal component of the specific discharge n E · q is assumed constant on each element face, F = ∂E i ∩ ∂E j , for two neighboring elements E i and E j of a topological triangulation T leading to conforming fluxes in RT N 0 space: in which q F is the specific discharge on F (Fig. 1). Because the normal flux ν F · q F is assumed to be constant on any F , and the source/sink strength f is assumed to be constant on E, the integrals over F and E reduce to a multiplication with |F | and |E|, respectively, assuming an RT N 0 space.
Higher-order approximations of the specific discharge on F are possible also in FVM [7,8].
A cell-centered FVM discretization of Eq. 1 results in a linear system of equations: in which A is the mobility matrix,ĥ is the vector of all approximated cell-related hydraulic head values, and r is a vector containing sources, and sinks integrated over the volume, as well as boundary conditions. In order to construct the mobility matrix A, the normal component of the hydraulic gradient ν F · ∇h at each face needs to be approximated from the head values in the cells. Towards this end, the standard procedure in FVM is to formally allocate the cell-related hydraulic head value to the cell center, compute the sum of distances of both cell centers to the face, and divide the difference of the cellrelated head values by this distance. In our FVM approach, we consider the centroids of triangles, or tetrahedra, as the cell centers rather than the circumcenters because the centroid is in general more representative for an element Fig. 1 Definitions exemplified for two adjacent triangles. The elements E i and E j share a face F on which a unit normal vector ν F following a sign convention is defined. The triangular face itself is defined by the line between the two nodes x n1 and x n2 . The coordinates of the nodes opposite of F are denoted P i for the node in E i and P j for the node in E j , respectively than the circumcenter, at least if hydraulic conductivity is isotropic. Moreover, the centroid always lies within the respective element, no matter how deformed it is. The disadvantage of taking the centroid as the characteristic point of an element in a finite volume formulation is that the line connecting the centroids of two neighboring elements is in general not orthogonal to their shared face, making an orthogonal projection of the distance vector onto the face necessary when computing the associated coefficients of the mobility matrix. Cordes and Kinzelbach [12] presented a finite volume method on triangles for isotropic hydraulic conductivity employing circumcenters. Such an approach is easier to compute than employing centroids, as no orthogonal projection is required, but the results are less accurate, and a valid application is restricted to triangles strictly meeting the Delaunay criterion. Otherwise, circumcenters may lie on an element face for right-angled triangles, which could introduce divisions by zero, if the corresponding face is also opposite to a right angle of the neighboring volume. The circumcenters of elements with angles larger than 90°lie outside of the elements, which can result in coefficients of the mobility matrix with the wrong sign.
For an accurate representation of Dirichlet boundaries and Neumann boundaries with a fixed, non-zero flux, we introduce ghost nodes by orthogonally projecting the centroid of an element next to the boundary, E ∈ E D ∪E N,± , onto the boundary. That is, we expand the discretization by introducing surface elements of dimension d − 1. We denote this expanded discretization by T D,N = T ∪ F D ∪ F N,± , where F D is the set of faces discretizing the Dirichlet boundary D , and F N,± is the union of the sets of faces discretizing the inflow and the outflow boundary, respectively.
The mobility coefficients related to a face of two neighboring elements of the expanded discretization with indices i and j is given by in which k is the element-wise isotropic hydraulic conductivity and c ∈ R is an arbitrary scalar with c > 0. Finite volume methods employing full material tensors have been developed for node-centered dual grids based on a primal grid of triangles, and deformed quadrilaterals ( e.g., [16]), and for general grids in two dimensions including cellcentered triangular grids [18]. All these methods include transformation of coordinates, which are beyond the scope of our present consideration. Furthermore, if l is the distance between two centroids, l i is the part of l in the element i, and l j is the part of l in the neighboring element j , and s ⊥ i,j is a scaling factor such that only the orthogonal part of a mobility coefficient normal to an internal face is considered: in which l is the distance vector between the centroids of two neighboring elements E i and E j , and f is the vector lying on F = ∂E i ∩ ∂E j . Then, the element A ik of the mobility matrix A ∈ R N×N with two indices i and j ranging over the elements in T D,N is given by and the right-hand side vector r ∈ R N is given by in which f i is the volumetric source/sink strength per unit volume for an element E i ,ĥ i,D is the assigned Dirichlet boundary condition on the Dirichlet face E i ∈ F D , and q f ix,i is the specific discharge on a Neumann face with non-zero flux,

Flux-postprocessing: the RT N 0 projection
Our flux correction method aims to project a nonconforming, non-mass-conservative velocity field onto a mass-conservative and conforming one. This is done by minimizing the element-wise residual between the hydraulic gradient of the P 1 Galerkin FEM solution and a solution corresponding to the velocity approximation in RT N 0 space, in which the normal component of the velocity is constant on and continuous across every face, and the velocity components vary linearly within an element. The RT N 0 basis functions, ψ F (x), can be set globally for the topological triangulation T for all elements sharing face F . Any face is either shared by two neighboring elements, F = ∂E i ∩∂E j , or it is a face of only one element, F ∈ F D ∪ F N . For each face F = ∂E i ∩ ∂E j , there is one vertex in both elements E i and E j each that is opposite of the face F , or F is a face of only one element, F ⊂ ∂E i . P i denotes the coordinates of the node opposite to F in element E i , and P j those of the node opposite to F in E j , if there is an element E j . Then, the basis function ψ F (x) is Then, the element-wise RT N 0 velocity field in global coordinates for a simplicial element E i is given by in which q ⊥ F = ν F · q F is the orthogonal component of the specific discharge vector on F , which is assumed to be constant on the face F ; a j is an individual scalar for each spatial direction x j ; and b is a scalar independent of the spatial direction equaling the divergence of the velocity field integrated over E i times (d|E i |) −1 ; therefore, b = 0 holds for non-divergent flow. It should be noted that also in cell-centered FVM fluxes normal to the faces, q ⊥ F are directly computable, such that the RT N 0 basis functions of Eq. 19 can be used to obtain the velocity field within the cells of an FVM. For a divergent velocity field, the RT N 0 velocity approximation varies within the element, and the same holds for an FVM velocity field with internal sources or sinks. For illustration, see Fig. 2 showing a divergent and a non-divergent velocity field. To make specific values comparable on an element by element basis, we evaluate the velocity at the element centroid, and denote such an element-wise velocity originating from the is evaluated at the centroid x c of an element E i , an element-wise residual vector q,E i between the RT N 0 -based velocity q RT N 0 E i and the element-wise constant P 1 Galerkin FEM velocity q F E E i can be determined: However, it is advantageous to consider the residual of the hydraulic gradient at the centroid instead of the specific discharge within the proposed global optimization procedure. If the specific discharge was considered directly, residuals in zones of higher hydraulic conductivity, which are small in relative terms, would have a higher contribution to the objective function of the optimization procedure than residuals in zones of lower hydraulic conductivity, which are high in relative terms. We propose to minimize the L 2 norm of the hydraulic gradient at the centroids x c , which is a slight adaption of the approach of Odsaeter et al. [33], who considered normal components of hydraulic gradients at the faces. Our residual vector at the centroid of the single element E i is in which ∇χ E i is the gradient of the linear shape functions χ E i in element E i . Considering Eq. 21, a global system of equations can be set up for the squared element-wise errors: in which the index i ranges over the elements and m is an index ranging over the faces, both according to a global numbering convention. Furthermore, T ∈ R Nd is the vector of element-wise errors in each spatial direction, where Nd is the number of elements times the number of dimensions and (q ⊥ F m ) F m ∈F ∈ R M is the vector of normal components of the specific discharge following global indexing of faces, where no-flow boundaries are excluded. In the following, ((∇h) F E E i ) E i ∈T ∈ R Nd is the vector of element-wise hydraulic gradients with (∇h) F E E i = ∇χ E iĥ obtained by a P 1 Galerkin FEM solution, and N ∈ R Nd×M is the matrix, in which M is the number of faces minus the number of faces belonging to Neumann boundaries containing all RT N 0 basis functions ψ F (x c ), evaluated at the centroids, scaled by K −1 E i according to a global numbering scheme of elements and faces: in which ψ F m (x c ) is the RT N 0 basis function of the face F m evaluated at the respective component of the centroid. N im has d elements, so that N im q ⊥ F m are the contributions of the normal flux at face F m to the hydraulic gradient vector at the centroid of element E i . The full matrix N ∈ R Nd×M is assembled from the individual contributions N im .
The squared errors of Eq. 22 are minimized under the constraint that mass is conserved element wise, leading to the following auxiliary condition: in which i ranges over all elements and m ranges over all faces following a global numbering scheme and f E i is the divergence of the specific discharge in element E i which equals the source/sink strength. M ∈ R N×M is a matrix with the following elements: Minimizing the squared residuals of Eq. 22 subject to the constraint of Eq. 24 is done by the method of Lagrangian multipliers, in which the following objective function is minimized: in which λ ∈ R N is the vector of element-wise Lagrangian multipliers. To obtain the minimum of the objective function, the derivatives with respect to (q ⊥ F m ) F m ∈F and λ must be zero: leading to the following system of equations in block matrix form: Equation 29 describes the system of equations of the RT N 0 projection, which is square and has the order M + N independent of the spatial dimensionality d. For incorporating a fixed, non-zero flux boundary condition with a normal flux component q ⊥ F m , we replace the corresponding entry on the main diagonal of the matrix by a one, set all other entries of the respective row to zero, and specify the normal component of the flux in the entry m of right-hand side vector.

Particle tracking methods
In advective particle tracking, the trajectories of ideal, point-like particles are computed via integrating the particle velocity over the travel time. A trajectory X p (t) = (x j,p (t)) j =1,...,d in x ∈ for a particle p is given by where t 0 is the starting time, t is a discrete time increment, v is a velocity field, and τ is the travel time. In the following, we present element-wise analytical solutions of Eq. 30 to compute the exit location of a particle on an element face from the start location of that particle within the element, or from the entry location on another element face, for nondivergent and divergent flow on arbitrarily shaped triangles and tetrahedra in global coordinates. With this, we can track a particle from one element face to the next throughout the whole domain.

Particle tracking in non-divergent flows on triangles
We first consider non-divergent flow on triangles, which implies an element-wise constant velocity vector. Given an entry point of a particle on the starting face F s , the element-wise constant velocity vector defines a straight line, starting at the entry point. The point of the first intersection of this line with one of the element faces not being the starting face is the exit point. It is computed by intersection of hyperplanes, i.e., lines (d = 2), or planes (d = 3), containing the element faces with the trajectory line. Let x p = (x p , y p ) be the starting point of a particle in an element E i , or on a face F S ∈ F E i . Then, the other points lying on the trajectory, given the velocity vector, where r > 0 is a positive scalar. Every face F of an element E i can be described by two nodes x n1 = (x n1 , y n1 ) and x n2 = (x n2 , y n2 ). Then, a possible exit point, x e,p = (x e,p , y e,p ), on the face F is given by line-line intersection: in which F s is the face on which the particle starts, if the particle starts on a face. Let (x j ) j =1,..,d = x d be the distance vector x e,p − x p , then the travel times leading to possible intersection points are simply given by The actual exit point x e is the possible exit point with the corresponding smallest positive travel time: in which i is the corresponding index of the possible exit point with the corresponding minimum travel time τ i .

Particle tracking in divergent flows on triangles
For divergent flow on triangles in global coordinates, we again consider the intersection of the particle trajectory within an element with the element faces. Considering the space (5) and Eqs. 18 and 19, we know that for every spatial dimension x j with x = (x, y), the velocity varies linearly within the element: Furthermore, we can compute the acceleration of the particle along its trajectory in the j th direction by in which t is time and p denotes the position of the particle. The derivatives of the right-hand side of Eq. 35 are Combining Eqs. 35 and 36 and given two distinct particle positions p 1 and p 2 on a trajectory within an element separated by a travel time t yield On every face, F ∈ F E i , of an element E i a vectorized line equation for every point, x l = (x l , y l ) on F , can be defined independent of the orientation of the line including the respective face: in which s is a scaling factor, and t l = (t x,l , t y,l ) is a tangential unit normal vector parallel to the respective face, and starting from one facial node, x n1 = (x n1 , y n1 ), pointing towards the second facial node, x n2 = (x n2 , y n2 ). Furthermore, we know that for two points p 1 and p 2 on a trajectory, the travel time t in Eq. 37 is the same no matter which spatial dimension x j is considered, such that t = t x = t y , and the factor b in Eq. 37 is element-wise constant. Considering these findings and Eq. 37, we observe that a x + bx e,p a x + bx p = a y + by e,p a y + by p , (39) in which x e,p = (x e,p , y e,p ) is a possible exit point of a particle on a face and x p = (x p , y p ) is an entry point of that particle. Then, a possible exit point on a face is given by substituting Eq. 38 into Eq. 39 yielding the scalar factor s such that Eq. 38 points to the possible exit point x e,p on a face: s = a y +by n1 a y + by p − a x + bx n1 a x + bx p · bt x,l a x + bx p − bt y,l a y +by p −1 .
Then, the vector of travel times τ = ( t) F ∈F E i \{F s } for divergent flow on triangles is given by evaluating the travel times in an arbitrary spatial dimension: Finally, the actual exit point is determined according to Eq. 33.

Particle tracking for non-divergent flow on tetrahedra
On tetrahedra, the velocity is again element-wise constant for non-divergent flow and a possible exit point of a particle on a face is given by simple line-plane intersection. Again, let x p = (x p , y p , z p ) be the starting point of a particle in an element E i , or on a face F s ∈ F E i , and let v n,E i = (v x,n , v y,n , v z,n ) be a vector of unit length pointing in the same direction as the element-wise velocity vector, v E i , and let n e = (n x,e , n y,e , n z,e ) be a normal vector on a possible exit face, , d y , d z ) be a distance vector between the entry point x p and an arbitrary point, p F = x ∈ F , on the respective face. Given the following scalar factor then the possible exit point on a face is given by The vector of travel times is in which x d = x e,p − x p . The actual exit point is again determined with Eq. 33.

Particle tracking for divergent flow on tetrahedra
For divergent flow, a particle tracking scheme on tetrahedra can be derived in analogy to those for triangles. On tetrahedra, the element-wise velocity field, with x = (x, y, z), is given by and every face is embedded in a plane equation such that where x l = (x l , y l , z l ) is the vector of coordinates of an arbitrary point on a face F , x n1 = (x n1 , y n1 , z n1 ) is the vector of coordinates of one of the nodes defining F , and p 1 = (p x1 , p y1 , p z1 ) and p 2 = (p x2 , p y2 , p z2 ) are the direction vectors of unit length parallel to F , which have to be linearly independent. For every point of a trajectory of a particle p, we again know that the travel time is the same no matter which spatial direction we consider t = t x = t y = t z . Reconsidering (Eq. 37), we observe that a x + bx e,p a x + bx p = a y + by e,p a y + by p = a z + bz e,p a z + bz p , in which x e,p = (x e,p , y e,p , z e,p ) is again a possible exit point of a particle on a face and x p = (x p , y p , z p ) is an entry point of a particle. Obviously, there are three possibilities to compute s for d = 3. We opt for s 1 such that t x = t y which yields furthermore, s 2 is set such that t x = t z which is then given by The factor t is computed by s 1 = s 2 yielding Possible exit points x e,p on planes including the faces F ∈ F E i \ {F s } for an element E i , are then uniquely determined by inserting t (Eq. 50) and either s 1 (Eq. 48) or s 2 (Eq. 49) for s in the plane (Eq. 46) for a face F . It is obvious that s 1 (Eq. 48), s 2 (Eq. 49), and t (Eq. 50) have, in any case, to be defined such that gaps in definition are avoided. That is, it has always to be checked, if the direction vectors t l for triangles and p 1 and p 2 for tetrahedra are chosen such that divisions by zero are avoided.
The vector of travel times, τ = ( t) F ∈F E i \{F s } , for divergent flow on tetrahedra is again given by evaluating the travel times in an arbitrary spatial dimension: Finally, the actual exit point is determined according to Eq. 33.

Particle tracking through the entire domain
In our particle tracking code, we track particles from one Dirichlet or Neumann boundary to the next. Given a starting location of a particle, we first search for the starting element. If particles are tracked only from Dirichlet or Neumann boundary conditions, we can also restrict the search to the respective elements sharing a boundary. If the start location is at an interface or a distinct node, we pragmatically choose the first possible starting element found in the search. Details on the search employing barycentric coordinates are given in the appendix. If the starting element is determined, all the following elements are determined, by tracking a particle from face to face. This scheme is in analogy to that of Pollock [34]. According to some criteria reflecting numerical accuracy, our particle tracking algorithm automatically switches between a tracking for non-divergent or divergent flow, depending on the properties of the element-wise flow field. The tracking procedure is continued until a particle reaches a Dirichlet or Neumann boundary.

Measures of differences between the discretization methods
We define difference measures comparing the specific discharge values at all element centroids of the RT N 0 projection with a reference solution being either analytical or obtained by the cell-centered FVM. The first difference measure abs is the ratio between the absolute values: in which q is the specific discharge evaluated at the centroid of element E i , obtained by the RT N 0 projection of the P 1 Galerkin FEM solution, and q Ref E i (x c ) is the specific discharge computed by the reference solution. The second difference measure quantifies the discrepancy in direction by evaluating the scaled angle between the specific discharge vectors: in which ∠(·, ·) denotes the angle in degrees between two vectors. For evaluating the overall behavior of a model, the element-wise measures abs and dir are assembled to vectors listing all element-wise values.

Implementation
The grids for the two-dimensional test cases are generated by the algorithm "triangle" [42], accessed via MeshPy

Non-divergent flow on triangles
We first compare the RT N 0 projection to the cell-centered FVM, and the original For the case of non-divergent flow, the latter implies that the calculated particle trajectories are streamlines, and any two neighboring trajectories are streamtubes with the same discharge. Figure 3a shows particle trajectories using the elementwise velocity approximation q F E E based on the P 1 Galerkin FEM. As discussed, specific discharge and average linear velocity fields, q(x) and v(x), respectively, are nonconforming, causing jumps of the normal components on faces. The erroneous physical interpretations of these jumps are virtual sources and sinks on the faces. These numerical sources and sinks are preserved in the tracking patterns. In the most extreme case, the element-wise velocity vectors of neighboring elements both point towards the same face. In such cases, we could not track the particles any further and we made the trajectories end on the respective face. Figure 3b shows the particle trajectories of the conforming cell-centered FVM solution of the groundwater flow (1) on triangles, whereas Fig. 3c shows the trajectories resulting from the RT N 0 projection of the P 1 Galerkin FEM velocity approximation. The trajectories of the cell-centered FVM solution are somewhat more angular with small rapid changes of the direction at cell interfaces. This is due to the strong constraints of the FVM discretization and has also been observed for mixed finite element solutions in RT N 0 [20,23,29,37]. In contrast, our RT N 0 projection leads to smoother trajectories, which we assess as being more correct in sub-domains of uniform hydraulic conductivity. We believe that the smoothness of the RT N 0 projection stems from the smoothness of the hydraulic gradients in the P 1 Galerkin FEM solution. While the latter needs to be corrected in order to obtain a mass-conservative RT N 0 solution of the velocity field, the correction is minimized within the constraints of a conforming, mass conservative velocity approximation. Of course, the velocity fields of both the FVM solution and our projection are in RT N 0 space, which are conforming and element-wise mass conservative. While the RT N 0 projection of the P 1 Galerkin FEM velocity approximation is smoother, the computational effort to obtain the cell-centered FVM solution is considerably lower.
We quantify the difference between the two conforming velocity fields depicted in Fig. 3b and c by the ratio of absolute velocities according to Eq. 52 and the angle between the two velocities according to Eq. 53. Figure 3d and e show the spatial distribution of these two difference measures, and Table 1 lists norms over the entire domain. As abs is a ratio, we give the respective infinity norm as ||(|1 − abs |)|| ∞ From this, we conclude that the velocities originating from our RT N 0 projection are very similar to those directly obtained from cell-centered FVM, with respect to both the absolute values and directions of the velocities throughout the domain, including the zones within or nearby the nearly impervious walls of very low hydraulic conductivity.
In particular, the RT N 0 projection does not introduce a noticeable numerical rotation, which might occur in other flux correction schemes [40], and occasionally also within the normal framework of mixed finite elements in RT N 0 space [20]. Note that even though both solutions are in RT N 0 space, one would expect different results as they are computed differently. As both solutions are numerical approximations, we cannot say which solution is more accurate. In general, the difference measures depicted in Table 1 are all very low and also depend on the numerical accuracy of the linear solvers used in the computations.

Divergent flow on triangles
In our next benchmark, we consider a 100 × 100 m domain with uniform isotropic hydraulic conductivity of k = 10 −5 (m/s), in which groundwater recharge takes place in a 30 × 30 m squared subdomain in the center of the domain, indicated by a blue-gray shaded square in Fig. 4. The recharge flux is q in = 200 (mm/a)≈ 6.3376 · 10 −9 (m/s). The domain is discretized by 1252 triangular elements and 667 nodes. The left and right boundaries of the domain are Dirichlet boundaries with a head difference of 1.0 (m), whereas no flow is allowed across the upper and lower boundaries. Figure 4a and b show the corresponding particle trajectories for the FVM velocity field and the RT N 0 projection of the P 1 Galerkin FEM velocity field, respectively. Both trajectory patterns appear reasonable. The trajectories diverge in the recharge area depicted by the blue-gray square. The trajectories of the FVM solution are again more angular, while the RT N 0 -projected solution preserves the smoothness of the P 1 Galerkin FEM solution.
Like in the non-divergent case, we quantify the difference of the two conforming velocity fields in the divergent case depicted in Fig. 4a and b by the ratio of absolute velocities according to Eq. 52 and the angle between the two velocities according to Eq. 53. Figure 4c and d show the spatial distributions of these two difference measures, and Table 2 lists norms over the entire domain. These metrics confirm that the two velocity fields are quite similar. All measures listed in Table 2 are very low. We also computed the infinity norm of the relative error in mass conservation, which was on the order of machine precision for double-precision real variables.

Empirical consistency and convergence tests for anisotropic hydraulic conductivity on triangles
For evaluating the performance of the RT N 0 projection in cases with anisotropic hydraulic conductivity, we consider a unit square = [0, 1] 2 in which a fixed non-zero flux is defined as Neumann boundary condition over the central sections (y ∈ [0.25, 0.75]) of the left and right boundaries, all other boundary sections are no-flow boundaries. The domain is discretized by 812 elements and 439 nodes, the maximum diameter of the elements, h T , is chosen such that 1/h T = 16. We consider a uniform hydraulic conductivity with the component k xx = 10 −4 (m/s) remaining identical in all cases, whereas the k yy component is varied such that k yy /k xx = 1, k yy /k xx = 0.1, and k yy /k xx = 0.01. The off-diagonal entries remain zero. The normal flux at the inflow and outflow boundaries is uniformly set to n · q N,+ = −n · q N,− = −10 −4 (m/s). As in the preceding case, the starting positions of particles are chosen such that the spacing reflects an equal cumulative flux across the boundary, so that the particle trajectories represent streamlines. For illustration, see Fig. 5, which illustrates the case of k yy /k xx = 0.1.
We show consistency by comparing the results of the RT N 0 projection to the corresponding analytical solution As abs is a ratio, we give the respective infinity norm as ||(|1 − abs |)|| ∞ rather than a numerical result based on cell-centered FVM. The analytical solution for our test case is in which Q (m 2 /s) is the total volumetric flux and L y (m) denotes the total vertical length of the domain, which equals unity in our case. Furthermore, q in,r k xx α r , and and in which L w is the length of the inflow/outflow boundary and L x [m] denotes the total horizontal length of the domain, which equals unity in our case. Figure 5a and b show the corresponding particle trajectories based on the analytical solution and the RT N 0 projection of the P 1 Galerkin FEM velocity field, respectively, for the case of an anisotropy ratio of 0.1. Figure 5c and d show the ratio of absolute velocities of the RT N 0 projection over the analytical solution and the difference in the direction, respectively. The trajectory pattern gives the visual impression that the result of the RT N 0 projection of the P 1 Galerkin FEM velocity field is very close to the analytical solution. The patterns of the two other anisotropy ratios show the comparable similarities, but the trajectories are more centered in the middle section for k yy /k xx = 0.01 and more evenly distributed in the vertical direction for k yy /k xx = 1. Table 3 lists the norms of difference measures, confirming the good agreement of the RT N 0 projection with the analytical solution for the isotropic case and for the case with the anisotropy ratio of k yy /k xx = 0.1. For an anisotropy ratio of k yy /k xx = 0.01, the RT N 0 projection of the P 1 Galerkin FEM velocity introduces numerical artifacts, which mainly occur in low-flow regions near the no-flow boundaries. Already the case of k yy /k xx = 0.1, Fig. 5c and d indicate these tendencies. For k yy /k xx = 0.01, however, the orientation of the velocities near the As abs is a ratio, we give the respective infinity norm as ||(|1 − abs |)|| ∞ top and the bottom is mostly reverted. In these regions, the absolute velocity is significantly smaller than that in the central section, particularly in cases of high anisotropy.
Because of the small values, getting the velocity right in these regions has a relatively low significance in the minimization procedure, which induces the described numerical artifacts.
To analyze the performance of the RT N 0 projection more rigorously, we extend this test by an empirical convergence analysis comparing our obtained results with those of the base P 1 Galerkin FEM solution. We employ the same models as exemplarily shown in Fig. 5 using iteratively refined grids. The base grid discussed above has a maximum diameter of the elements, h T , such that 1/h T = 16. For the empirical convergence analysis, we use this grid among other refinement levels, whereas keeping all adjustments for the grid generation identical on all levels.  The empirical order of convergence is given in parentheses. h T is the maximum diameter of the triangles, i.e., the maximum face length, and N is the number of elements. Furthermore, q RT N 0 and q Ana are the velocity fields of the RT N 0 projection and the analytical solution, respectively According to Table 4, the velocity field approximated by the RT N 0 projection essentially shows the same order of convergence as the original velocity field based on the P 1 Galerkin FEM. It is also worth noting that, besides of the coarsest grid for an anisotropy ratio of unity, all other norms indicate that the conforming velocity field of the RT N 0 projection is overall closer to the analytical solution than the original flux approximation of the P 1 Galerkin FEM. Nevertheless, the numerical artifacts of reverted velocity fields and wrong magnitudes of the velocities induced by the RT N 0 projection for an anisotropy ratio of k yy /k xx = 0.01 are also preserved in the reduced empirical convergence order shown in the first refinement step in Table 4.
We conclude the test case by reporting on the computation time needed to solve the equation systems resulting from the RT N 0 projection for the isotropic case, for which we can also perform the cell-centered FVM simulations as a comparison. The systems of equations are solved by the direct solver UMFPACK [14], accessed via the backslash operator of Matlab.
The CPU times shown in Table 5 clearly indicate that the final system of equations resulting from the cell-centered FVM can be solved quicker for larger systems than the system resulting from the RT N 0 projection. This is so because the system of equations resulting from the cellcentered FVM discretization is symmetric, positive definite, and of order N, whereas the system of equations resulting from the RT N 0 projection is, though symmetric, an indefinite saddle point problem of order N + M comparable to an RT N 0 MFEM discretization.

Three-dimensional test case
Our three-dimensional benchmark model is a cube with an edge length of 10.0 (m). The background isotropic hydraulic conductivity is k = 10 −4 (m/s). Two cuboid inclusions ( x ∈ [2,4], y ∈ [0, 5], z ∈ [2, 4] and x ∈ [6,8], y ∈ [5,10], z ∈ [6,8], respectively) have a strongly reduced isotropic hydraulic conductivity of k = 10 −10 (m/s). In Fig. 6a, the inclusions are highlighted as blue blocks. The uniform porosity is = 0.4. The domain is discretized into cubes of 1 (m) edge length, which are split into six tetrahedra such that effects of grid orientation in the cellcentered FVM solution is minimized. This leads to 6000   Figure 6a shows a few particle trajectories starting at the inlet boundary of the domain, which have been computed using the RT N 0 projection of the specific discharge field originating from the P 1 Galerkin FEM solution. The particle trajectories circumvent the inclusions in a physically reasonable manner. In contrast to the twodimensional test cases, we do not show particle trajectories for the FVM solution, because the visual impression does not indicate any obvious differences to trajectories obtained from the RT N 0 projected velocity field of the standard FEM solution. Figure 6b and c visualize the measures of differences in the velocity field between the FVM solution and the velocity approximation by our RT N 0 projection. Figure 6b shows the ratio of absolute velocities, and Fig. 6c the normalized angle between the two velocity vectors at the cell centroids. Table 6 lists the associated norms for the entire domain. The magnitude of the differences are visualized in Fig. 6b and c by both the size and color of spheres at the cell centroids. In the largest part of the domain, these differences are negligible. Only at certain points next to the inclusions the difference measures exceed average values. The difference Table 6 Measures for the differences between the discretization methods in the 3-D application with two nearly impervious cuboid inclusions Mean(·) Median(·) ||(| · |)|| ∞ abs 9.89 · 10 −1 9.74 · 10 −1 8.34 dir 5.83 · 10 −2 5.86 · 10 −2 1 As abs is a ratio, we give the respective infinity norm as ||(|1 − abs |)|| ∞ in the direction is mainly caused by a slight shift of points where the particle trajectories separate to circumvent the inclusions. By coincidence, two cells show a larger difference in the absolute value of the velocity. As indicated by the difference metrics listed in Table 6, the results are very similar. Like in the 2-D test cases, the trajectories based on the FVM solution on simplices are more angular than those based on the RT N 0 projection of the P 1 Galerkin FEM solution, because the projection tries to preserve the smoothness of the standard FEM solution. We would deem the projected results to be physically more accurate, but the FVM solution is achieved at considerably lower computational costs.

Conclusions
We have proposed an RT N 0 projection of velocities obtained from solving the groundwater flow equation by P 1 Galerkin FEM. Our projection yields physically reasonable flow fields preserving the smoothness of the P 1 Galerkin FEM solution for isotropic test cases and mild anisotropy. For comparison, we have presented an analytical solution and a cell-centered finite volume formulation on simplices in which we account for Dirichlet and Neumann boundary conditions by ghost nodes. The results of our projection are similar to those obtained by an analytical solution and the cell-centered finite volumes for two-and three-dimensional isotropic and mildly anisotropic flow fields. The velocity approximations of our RT N 0 projection and of the FVM are in the same function space, and the velocity fields are similar in magnitude and the direction for the evaluated cases.
We claim that our flux correction method, like that of Odsaeter et al. [33], gives better results than the postprocessing methods of Schiavazzi [40] and Scudeler et al. [41], in which velocity differences between the P 1 Galerkin solution and the RT N 0 projection showed much larger relative differences in low-flow regions. The latter can be avoided by minimizing differences in hydraulic gradients rather than velocities between the P 1 Galerkin FEM and the RT N 0 solutions. While Odsaeter et al. [33] considered hydraulic gradients on element faces, we minimize the difference in hydraulic gradients at the element centroids. Our approach yields physically reasonable velocity fields regarding the magnitude of the velocity vectors and does not introduce erroneous numerical rotation, neither for two-dimensional nor for three-dimensional flow fields with isotropic or mildly anisotropic hydraulic conductivity, which seems to occur in other flux correction schemes even in the isotropic case [40], and have occasionally been observed in RT N 0 -based mixed (hybrid) FEM schemes [20]. However, if stronger anisotropy is considered, our RT N 0 projection may numerically induce rotation and wrong magnitudes of the velocities leading to partial failure of the scheme in low-flow regions. Similar problems were already observed occasionally for former flux corrections, but here even in isotropic cases, which were evaluated by Schiavazzi [40] and applied by Scudeler et al. [41].
By construction, the velocity fields of our RT N 0 projection are conforming and mass conservative. In comparison with velocity fields obtained from cell-centered FVM, the projected fields are smoother and trajectories are less angular so that we assess the projected velocity fields to be physically more accurate than the fields obtained by the FVM scheme, at least in regions of uniform hydraulic conductivity. However, the computational costs of a cellcentered FVM scheme are much lower than those of our RT N 0 projection. For an FVM solution, one system of equations of the order of elements has to be solved, in which the matrix is symmetric and positive definite. By contrast, the projection requires first solving a system of equations of the order of the number of nodes involving a symmetric, positive-definite matrix to obtain the P 1 Galerkin FEM solution, followed by solving a saddle point problem in the RT N 0 projection step, in which the matrix is symmetric, non-definite, and of the order of the number of elements plus the number of faces without no-flow boundary conditions. Both velocity fields are in RT N 0 space. For large-scale applications, the computational advantage of the FVM scheme may outweigh the smoothness of the RT N 0 projection. Conversely, two variable velocity fields may cause erroneous twisting of streamlines in heterogeneous formations, which can have significant effects on the simulation of transverse mixing [6,9].
We have presented semi-analytical particle tracking methods for conforming RT N 0 velocity fields on triangles and tetrahedra. Extending the RT N 0 projection, FVM discretization, and particle tracking scheme to other element types (e.g., triangular prisms, pyramids, hexaeder) is possible, but addressing deformed elements would require formulating all schemes in local coordinates applying a contravariant Piola transformation [39]. Projecting velocity fields obtained by other finite element methods like higher-order Galerkin FEM onto a higher-order conforming velocity field would be possible as well, although one might opt to minimize weighted differences in the hydraulic gradient at more than a single point per element in such cases, but these considerations are beyond the scope of the present analysis.
We have developed the presented RT N 0 projection for steady-state, potentially divergent flow including mild anisotropy. An extension to transient flow is probably possible by treating the change of storage like a source/sink term in steady-state flow. In order to do this, however, the nodal loads computed by P 1 Galerkin FEM must be distributed to the elements associated with the nodes. Likewise, extensions to unsaturated flow or other nonlinear flow laws would require modifications, although the complexity of the final equation system to solve would remain the same. Especially for nonlinear flow laws, the application of an RT N 0 projection to postprocess an existent non-conforming flow solution might still be computationally cheaper than computing a conforming flow solution from scratch using an RT N 0 MFEM or comparable methods. As a postprocessor, the RT N 0 projection has to be performed only at time points at which a velocity field is wanted, and the nonlinear dependence of hydraulic conductivity on the simulated head and saturation fields would be evaluated within the primary solution step, deeming the projection itself linear.
Acknowledgment We like to thank Ole Klein from the Interdisciplinary Center for Scientific Computing at Heidelberg University for a very insightful discussion on finite element methods and the numerics behind a possible RT N 0 projection. Furthermore, we are grateful to Max Allmendinger for proofreading the mathematical notation and several joyful discussions on mathematical topics. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/. Appendix: Using barycentric coordinates to find the starting element of a particle trajectory Every pointx s within a simplex can be described by a linear combination of the coordinates of the nodes,x i , and a barycentric, nodal weight, β i ≥ 0, respectively, such that

Funding information
with i∈N E β i = 1, and β j = 1 − in which j is the index of a node and β j is its associated weight, such that every point in the simplex can be described by the coordinates of its nodes and d + 1 weights. Combining (57) and (58) leads to in which T is a transformation matrix only depending on the coordinates of the nodes; T −1 can easily be evaluated analytically; β is the vector of the d independent, barycentric, nodal weights; andx j is the vector of coordinates of node j of the element. We exploit the concept of barycentric coordinates in the search for the element in which the starting point of our particle tracking scheme resides. To do so, we set x p =x s and solve equation (60) for every element. If a particle lies within an element, all weights β = (β i ) i=1,...,d+1 are within the interval 0 ≤ β i ≤ 1 and sum up to unity. We stop the search at the first instance at which both criteria are met.