A combined GDM--ELLAM--MMOC scheme for advection dominated PDEs

We propose a combination of the Eulerian Lagrangian Localised Adjoint Method (ELLAM) and the Modified Method of Characteristics (MMOC) for time-dependent advection-domina\-ted PDEs. The combined scheme, so-called GEM scheme, takes advantages of both ELLAM scheme (mass conservation) and MMOC scheme (easier computations), while at the same time avoids their disadvantages (respectively, harder tracking around the injection regions, and loss of mass). We present a precise analysis of mass conservation properties for these three schemes, and after achieving global mass balance, an adjustment yielding local volume conservation is then proposed. Numerical results for all three schemes are then compared, illustrating the advantages of the GEM scheme. A convergence result of the MMOC scheme, motivated by our previous work on the convergence of ELLAM schemes, is provided, which can be extended to obtain the convergence of GEM scheme.


Models and assumptions
1.1. Introduction. In this paper, we consider a time-dependent advection-dominated PDE (1), and study some numerical schemes for this equation that are based on characteristic methods. These types of PDEs are encountered in many important fields, such as mathematical models in porous medium flow (e.g. reservoir simulation), and fluid dynamics (e.g. Navier--Stokes equations). A short summary, which includes most of the commonly used numerical schemes for (1), together with their advantages and disadvantages, have been presented in [17].
In particular, our work focuses on two types of numerical schemes based on characteristic methods which are popularly used, namely the Eulerian Lagrangian Localised Adjoint Method (ELLAM) and the Modified Method of Characteristics (MMOC). The advantages of these schemes lie on the fact that they are based on characteristic methods, and thus capture the advective component of the PDE better than upwinding schemes. Moreover, these schemes are not limited by CFL constraints, and hence large time steps can be taken for numerical simulations. These are usually combined with finite difference (FD), finite element (FE) or finite volume (FV) discretisations, in order to provide a complete numerical scheme for (1). To cite a few examples, the FE-MMOC [13], FE-ELLAM [5], and FV-ELLAM [19], have been used to discretise (1). Other variants of the ELLAM, as well as a summary of its properties, have also been presented in [24]. More recent variants of the ELLAM involve the volume corrected characteristics mixed method (VCCMM) [1,4]. Aside from the global mass conservation property of ELLAM, these ensure that local volume conservation is achieved. On the other hand, more recent studies of the MMOC involves MMOC with adjusted advection (MMOCAA) [11]. Compared to the MMOC, MMOCAA has better global mass conservation properties, which is usually required for an accurate numerical simulation of models that are related to engineering problems. On the other hand, of particular difficulty in the implementation of ELLAM is an accurate evaluation of integrals involving steep back-tracked functions [26]. An inaccurate evaluation of these integrals will yield a loss in mass conservation; a fix to simplify the evaluation of these integrals, which will preserve the mass conservation property, has been proposed recently in [4].
Here, we present a detailed analysis in terms of the mass conservation properties, and provide an idea of combining the ELLAM and MMOC schemes, so that we can capitalise on the advantage of both schemes, while at the same time avoid their disadvantages. Having achieved global mass balance, a novel, less expensive adjustment yielding local volume conservation is then proposed. The complete coupled scheme then consists of a characteristic component (the combined ELLAM-MMOC), accompanied by a discretisation of the diffusive terms using the Gradient Discretisation Method (GDM) framework [14]. This framework contains many classical methods (finite elements, finite volumes, etc.) for diffusion equations. The complete coupled scheme, named GEM (for GDM-ELLAM-MMOC), therefore presents in one form several possible discretisations of the advection-diffusion model.
The main contributions of this work are • a new adjustment algorithm to achieve local volume conservation (Section 2.5), • precise analysis of mass balance errors for MMOC scheme (Section 3.4), • combined ELLAM-MMOC scheme for the advection-reaction equation (2) (Section 4), • application of the ELLAM-MMOC scheme and local volume adjustments for the miscible flow model (GEM scheme, Section 5).
The paper is made up of two major components. The first component focuses on the advection-reaction equation (2) and the characteristic based schemes used to discretise this equation. We start here by studying the ELLAM scheme and the MMOC scheme in Sections 2 and 3 respectively. In particular, we look into the global mass balance properties of both schemes, and then discuss, for the ELLAM, how to achieve local volume conservation after achieving global mass balance. We then propose in Section 4 a combined ELLAM-MMOC scheme, and discuss its advantages over both the ELLAM and the MMOC. The second component then focuses on the application of the proposed ELLAM-MMOC combination to a miscible flow model (28) in porous medium. This model involves diffusion terms, that are discretised using the generic GDM framework. The local volume adjustments performed in Section 2.5 are then adapted for the numerical schemes for model (28). Numerical results presented in Section 5.4 illustrate the clear advantages of this GEM scheme. Finally, we provide in Section 6 a convergence result for the MMOC scheme. This can then be extended, together with the convergence result for the ELLAM scheme [9], to obtain the convergence of the ELLAM-MMOC scheme for (1).

Models.
Our objective is to design a robust, characteristic-based numerical scheme for a model of miscible displacement in a porous medium. This model, described in Section 5, involves an elliptic equation for the pressure, and an advectionreaction-diffusion equation for the concentration of the invading fluid. For simplicity, we describe the characteristic-based scheme for the concentration equation without explicitly referring to the pressure equation. We therefore consider the scalar model in which T > 0, Ω is an open bounded domain of R d (d ≥ 1), the porosity φ, the diffusion tensor Λ and the velocity u are given, u · n = 0 on ∂Ω, and f (c) = f (c, x, t) is a function R × Q T → R. The unknown c(x, t) represents the amount of material (a fraction) present at (x, t). The characteristic method only deals with the advective part of the model, and will therefore be described on the advectionreaction equation (corresponding to Λ ≡ 0): Note that the boundary is non-characteristic due to the assumption u · n = 0 on ∂Ω, and thus no boundary conditions need to be enforced in (2).

1.3.
Assumptions on the data, and numerical setting. Throughout the article we assume the following properties: φ ∈ L ∞ (Ω) and there exists φ * > 0 s.t. φ ≥ φ * a.e. on Ω (3b) u ∈ L ∞ (0, T ; L 2 (Ω) d ) and divu ∈ L ∞ (Q T ). (3c) Our objective in this paper is to describe numerical methods in a general setting, to ensure that our design and analysis of ELLAM-MMOC schemes applies at once to various possible spatial discretisations (e.g. finite-element or finite-volume based). To achieve this, we use the Gradient Discretisation Method (GDM), a generic framework for numerical methods for diffusion equations [14]. Although most of our work will be done here on the advective-reactive parts of (1), we will demonstrate that the GDM also provides all the required tools to describe ELLAM and MMOC schemes.
The GDM consists in replacing, in weak formulations of the models, the continuous (infinite-dimensional) spaces and corresponding operators by a discrete (finitedimensional) space and reconstructions of functions and gradients. A space-time Gradient Discretisation (GD) is C = (X C , Π C , ∇ C , I C , (t (n) ) n=0,...,N ), where • X C is a finite-dimensional real space, describing the unknowns of the chosen scheme, • Π C : X C → L ∞ (Ω) is a linear operator that reconstructs a function from the unknowns, • ∇ C : X C → L ∞ (Ω) d is a linear operator that reconstructs a gradient from the unknowns, • I C c ini is a rule to interpolate c ini onto an element of X C , • (t (n) ) n=0,...,N are the time steps, and we let δt (n+ 1 2 ) = t (n+1) − t (n) . Different choices of C lead to different schemes, e.g. finite elements or finite volumes [14].
Finally, we assume in the following that u is approximated on each time interval (t (n) , t (n+1) ) by a function Precise assumptions regarding this approximation will be described in Section 6.
In the rest of the paper, the variables are only made explicit in the integrals when there is a risk of confusion. Otherwise we simply write, e.g., Ω qdx.

ELLAM scheme for the advection-reaction equation
2.1. Motivation. For any sufficiently smooth function ϕ, the product rule yields Hence, (2) gives, for any time interval (t (n) , t (n+1) ), To simplify the second term on the left hand side of the above equation, the ELLAM requires that test functions ϕ satisfy with ϕ(·, t (n+1) ) given. The advection equation (2) then leads to the relation 2.2. ELLAM scheme. The ELLAM scheme consists in exploiting the motivation above, in the discrete context of the GDM in which trial and test functions are replaced by reconstructions Π C applied to trial and test vectors in X C .
Definition 2.1 (ELLAM scheme). Given a Gradient Discretisation C and using a weighted trapezoid rule with weight w ∈ [0, 1] for the time-integration of the source term, the ELLAM scheme for (2) reads as: find (c (n) ) n=0,...,N ∈ X N +1 C such that c (0) = I C c ini and, for all n = 0, . . . , N − 1, c (n+1) satisfies where v z is the solution to Here and in the rest of the paper, we let f k := f (Π C c (k) , ·, t (k) ) (or with a suitable average over (t (k) , t (k+1) ) if f is not continuous in time).
Remark 2.2 (About the time integration). The velocity field u was approximated by its value at time t (n+1) , given by u (n+1) . Other choices for the approximation of u, such as a centred approximation 1 2 (u (n) + u (n+1) ), may be made, but we noticed in our tests that this does not noticeably change the numerical solution. A weighted trapezoid rule is applied for time integration in Definition 2.1 for the purpose of achieving mass conservation, as discussed in [1].
Define the flow F t : Ω → Ω such that, for a.e. x ∈ Ω, Under Assumption (4), the existence of this flow is proved in [9,Lemma 5.1]. The solution to (7) is then understood in the sense: for t ∈ (t (n) , t (n+1) ] and a.e.
For any functions f and g, defining the vector functions f (n,w) and g F by the time-stepping (6) can be rewritten in the condensed form 2.3. Physical Interpretation. We provide a simple physical interpretation of the ELLAM, by supposing to simplify that Π C is a piecewise-constant reconstruction on a given mesh M. We also assume that for each cell K ∈ M, there is M 1 M and taking z K as test function, (11) and (9) give which reduces to where |E| φ = E φ is the available porous volume in a set E ⊂ R d . The first term on the right hand side of (12) tells us that the amount of material c (n+1) K present in a particular cell K ∈ M at time t (n+1) is obtained by locating where the material in cell K comes from, hence back-tracking the cell K to F −δt (n+1/2) (K), measuring how much of the material c (n) M is taken from each M ∈ M, and deposing this material into the cell K. These are accompanied by the contribution of the source term f in the particular cell K, which is given by the second term. We note here that this second term has a very similar treatment as the first term, i.e. the contribution that comes from f at time t (n) is determined by the traceback region associated to cell K.

2.4.
Mass balance properties. One desirable property for numerical schemes is conservation of mass. Essentially, we want a discrete form of the following equation, obtained by integrating (2) over Ω and which tells us that the change in c is dictated by the amount of inflow/outflow given by the source term.
Hence, it is important to define a measure of the mass balance error. Setting e := (1, 1), the (discrete) mass balance error is defined by A mass balance-preserving method is one for which e mass = 0. The ELLAM scheme (6) satisfies this property. Indeed, taking z 1 = K∈M z K , which satisfies Π C z 1 = 1 over Ω, as a test function in (11) gives e mass = 0.

Remark 2.3 (Steep back-tracked functions).
The natural physical interpretation of ELLAM, together with its mass conservation property, seem to indicate that the ELLAM scheme should be preferred over other numerical schemes for the advection equation (2). However, for Darcy velocities typically encountered in reservoir simulations, the streamlines of the flow F t concentrate around injection wells, and the functions v z defined by (9) are then extremely steep in these regions. An accurate approximation of the integral of these functions in cells close to the injection well then requires to track a lot of quadrature points, which is very costly [26]. In some instances, even tracking several points along these regions would not give an accurate depiction of the integral. This is one of the main issues with ELLAM implementations. Fixes have been proposed, but they consist is resorting to a different approach, near the injection wells, than the ELLAM process [4]. We aim at designing a numerical scheme that readily behaves well, without having to implement specific fixes in certain regions. The MMOC will be instrumental to that objective.

2.5.
Local volume conservation. The main difficulty of implementing an EL-LAM type scheme is the evaluation of the integral K φΠ C c (n) dx for each cell K, where K = F −δt (n+ 1 2 ) (K). In general, the region K (see Figure 1, left) cannot be exactly described, and hence we approximate it by polygons obtained from backtracking the vertices, together with a number of points along the edges of the cell K. Figure 1 (right) gives an illustration of the approximate traceback region K obtained by tracking the vertices, together with the edge midpoints of the cell K.
In general, | K| = |F −δt (n+ 1 2 ) (K)|. However, the equality of these volumes is essential in numerical simulations to preserve the accuracy of numerical solutions. To illustrate this point, consider the simple case of a divergence free velocity field in (2), with φ = 1, f = 0 and c ini = 1. In this test case, the exact solution is given by c(x, t) = 1. In theory, upon implementing a numerical scheme with piecewise constant approximations for the unknown c, we should have the following simplified form of (12) at the first time step: However, due to the approximation of the traceback region, we only have in general. This example shows that an inaccurate approximation of the volume of the tracked cell renders the numerical scheme unable to recover constant solutions. Hence, we need to perform some adjustments on the polygonal region K in order to yield | K| = |F −δt (n+ 1 2 ) (K)|, which we shall define as the local volume constraint for K. In [1,2], the local volume constraint is satisfied by moving through the tracked cells one by one, and adjusting the locations of the tracked midpoints along the direction of the characteristics. This was illustrated to work for square cells; however, there is no guarantee that such adjustments will yield a valid mesh configuration for irregular cells. A more generic approach is taken in [10], where the local volume constraint is satisfied by solving an optimisation problem of finding a mesh which is closest to the tracked mesh formed by K∈M K, satisfying the local volume constraints. Common to these algorithms is an explicit expression for the adjusted traceback regions. However, as can be seen in (12), this is not necessary for piecewise constant approximations, standard in reservoir simulations based on finite volume methods. The important aspect is the computation of the quantities |M ∩ K|.
We propose an algorithm which adjusts |M ∩ K| for each cell K, in the sense that these adjusted volumes would be something we expect to recover from a mesh obtained by adjusting the tracked points. The proposed algorithm works on any type of cells, but for simplicity of exposure we illustrate it in Figure 2 using square cells, and traceback regions K i approximated by tracking the vertices and edge midpoints of K i . In Figure 2, the blue lines denote the trajectory taken by the velocity field, and the red squares form part of the original mesh cells K ∈ M, whereas the black cells are their traceback regions: shaded points are the tracked vertices, and the hollow points are the tracked edge midpoints. In practice, after performing the tracking, aside from the final location of these vertices and midpoints, we also store the cell that they finally reside in. The algorithm is implemented cell-wise, starting from the first cell K 1 , and proceeds as follows: Consider a cell K 1 with neighbors K 2 , K 3 , etc. This leads to a tracked cell K 1 with neighbors K 2 , K 3 , etc. Suppose that K 1 intersects the residing cells M 1 , M 2 , M 3 and M 4 (see Figure 2, left). i) We start by measuring the error e K1 : The relation e K1 > 0 (resp. e K1 < 0) means that we overestimate (resp. underestimate) the volume of the tracked region. ii) We then compute the magnitude |u| of u at the tracked midpoints and also check whether u points into K 1 or not. If the velocity points towards the same direction for two consecutive midpoints, then we also compute the magnitude of u for the vertex in between them. iii) We now illustrate how to adjust the volumes of the regions. If e K1 < 0, then it means that the volume | K 1 | should be increased. Based on the velocity field u in Figure 2, this should be done by increasing along the direction of K 2 and K 3 . Now, the velocity field along the edge midpoints that are located at K 1 ∩M 2 and K 1 ∩M 3 points outward of K 1 and hence the vertex at K 1 ∩ M 4 should also be included. Denote then by |u 1,i | the magnitude of the velocity field evaluated at these tracked points in K 1 ∩ M i , i = 2, 3, 4. We will then adjust each of the volumes by subtracting, to each | K 1 ∩ M i |, the quantity |u1,i| 4 j=2 |u1,j | e K1 . This will then make K 1 satisfy the local volume constraint. iv) To make sure that this quantity really represents something that would have come from a perturbed mesh (see Figure 2, right), the quantities |u1,2| 4 j=2 |u1,j | e K1 should be added onto | K 2 ∩ M 2 |, and the corresponding quantities for the other edges).
We then proceed to adjust the other K i in the same manner so that they satisfy their respective local volume constraints.
As a remark, we note that one set of adjustments may not suffice to satisfy the local volume constraints for all cells. For example, if, after the adjustment of K 1 , we have that e K2 > 0, then we need to decrease the volume of K 2 . This can only be done (respecting the direction of the velocity field) by moving along the direction of K 1 . This will lead to e K1 > 0. Hence, after all the adjustments in the first stage, we need to re-evaluate e Ki and re-adjust the volumes. Of course, from a computational point of view, not all cells would be able to satisfy e Ki = 0 exactly, so we terminate our algorithm once |e Ki | is below a desired tolerance value for all cells. Another potential issue that may be encountered is when K 1 is a boundary cell (see Figure  3). If K 1 lies on the boundary and e K1 > 0, then, we need adjustments which will decrease the volume of K 1 . Thus, adjusting along the direction of the velocity u will only worsen the problem, since it will further increase the volume of K 1 and hence the value of e K1 . In this case, the idea is to consider −u instead. This translates to pushing inwards the points that would have been pushed outwards if e K1 < 0. The re-adjustment of the other cells follow accordingly, and convergence is still expected. Making the volume adjustments in the direction of u (resp. −u) corresponds to saying that we have tracked backward (resp. forward) too much, and hence to fix these, we must track forward (resp. backward) a little bit further. As has been mentioned in [8], when dealing with irregular cells, we need to track more than just edge midpoints so that K is close to F −δt (n+ 1 2 ) (K). Hence, for irregular cells, a slight modification for our algorithm should be made: we may opt to take smaller time steps (to make sure that the errors e Ki are small to start with, and adjustments can be made in a similar manner as with square cells, by only placing markers on the tracked midpoints and adjusting accordingly), or we may mark more than just the tracked midpoints (in order to have a better idea of the geometry of the tracked cell, and provide a more comprehensive list of the cell volumes to adjust). For our tests in Section 5.4, when such a modification was required, we chose to take smaller time steps.

MMOC scheme for the advection-reaction equation
3.1. Motivation. We use the product rule to write div(uc) = cdiv(u) + u · ∇c. By treating φ ∂c ∂t + u · ∇c as a directional derivative, and denoting by τ the associated characteristic direction, we rewrite (2) as follows where ζ = (φ 2 + |u| 2 ) 1 2 . The MMOC then considers the following approximation: Here 1 2 ) and the flow F t is defined by (8).
Definition 3.1 (MMOC scheme). Given a Gradient Discretisation C and using a weighted trapezoid rule with weight w ∈ [0, 1] for the time-integration of the source term, the MMOC scheme for (2) reads as: find (c (n) ) n=0,...,N ∈ X N +1 C such that c (0) = I C c ini and, for all n = 0, . . . , where we recall that e := (1, 1), and where we have set (by slight abuse of the notation (10)) 3.3. Physical Interpretation. As with the ELLAM, an interpretation of the MMOC will be provided for the simple case wherein we have a piecewise constant approximation of c. Fixing K ∈ M and taking in (18) the test vector z K , and thus The first term on the right hand side of the equation tells us that the amount of material c (n+1) K present in a particular cell K ∈ M at time t (n+1) is obtained by taking all cells M ∈ M, advecting material from each of these cells (by computing the trace-forward regions F δt (n+1/2) (M )), and determining which portion of each cell flows into K. The second term simply represents the change that comes from the source term f . We note that, unlike in the ELLAM, the contribution of the source term f for the MMOC is taken exactly to be from cell K. By itself, this term tells us that, if the source term f is nonconstant over regions close to one another, the MMOC will give either an excess or miss some amount that has flowed into the region K. The third term in (20) represents taking away a percentage of the net inflow/outflow in cell K and, in some sense, attempts to balance out the excessive or missing amount resulting from the second term.
Physically, the equivalence is expected, as we are now just comparing the first terms of equations (12) and (20), which both compute the amount of substance that has flowed into cell K. This can be done in two ways: Either we first locate the regions from which the substance has come from (ELLAM), or we let the substances in all cells flow, and determine which ones enter the cell K (MMOC). Mathematically, this equivalence can be established by performing a change of variables in |F δt (n+ 1 2 ) (M ) ∩ K| φ and using the property φ(F (22)).

3.4.
Analysis of mass balance error. Consider the MMOC scheme (18). By taking the test function z 1 = K∈M z K , we have Π C z 1 = 1 in Ω and we obtain thus the discrete mass balance equation From this, we see that one of the disadvantages of MMOC schemes over ELLAM schemes is that, in general, MMOC schemes do not preserve mass. The mass balance error e mass for MMOC is estimated by using (21) to substitute Ω φΠ C c (n+1) dx in (14). Performing a change of variables, we obtain It follows from [9, Lemma 5.2] that Hence, By the triangle inequality and recalling the definition (19) of (Π C c) (n,w) , we infer This estimate shows that the mass balance error e mass is minimal when δt (n+ 1 2 ) tends to 0 (as F −t → Id as t → 0) or the approximate amount of substance c near the non divergence-free regions, denoted by U , is almost constant. More precisely, Due to this, the MMOC is in most cases expected to give a relatively large mass balance error compared to ELLAM. Hence, it is not practical to perform the local adjustments described in Section 2.5 for the MMOC scheme.

Remark 3.3 (Conservation of mass for the MMOC)
. Estimate (24) shows that if the velocity field is divergence free then the MMOC scheme conserves mass, which is consistent with Remark 3.2.

Remark 3.4 (Forward tracking and cost near the injection cells).
Contrary to the ELLAM, the MMOC requires to forward-track test functions (see (20) in the case of piecewise constant approximations). Hence, in the MMOC, functions whose support is near the injection wells are not backward tracked into the injection cells, which makes them very steep and difficult to integrate (see Remark 2.3), but forward tracked far from these cells into non-steep functions that are easier to integrate.

A combined ELLAM-MMOC scheme for the advection-reaction equation
Here, we propose a combined ELLAM-MMOC scheme, to benefit from the mass balance property of the ELLAM and mitigate its costly implementation near the injection wells by using the MMOC method, much less expensive in these regions.
We start by applying a pure ELLAM scheme over the first few time steps, until c is almost constant in areas near the non divergence-free regions. After which, we do a split ELLAM-MMOC scheme, where we apply MMOC over these areas, and ELLAM elsewhere. The interest of such a scheme is twofold. First, the computational cost is reduced compared to a pure ELLAM scheme as we no longer have to compute integrals of steep functions. Second, upon using MMOC only in regions where divu = 0 or c is already almost constant, no mass balance error occurs. This combined scheme takes out the main disadvantages of both methods.
Discretise this by applying ELLAM on the first part αc (and αf ) and MMOC on the second part (1 − α)c (and (1 − α)f ). Defining Π C (αc) (n) as αΠ D c (n) , one time step of this leads to Remark 4.1 (Interpretation of the combined ELLAM-MMOC). An interpretation can be given by considering c 1 = αc and c 2 = (1 − α)c as two miscible fluids (that are also miscible in their surroundings) that do not react with each other, and are advected by the velocity u. We can consider the combination of these two as one single fluid with concentration c, that is advected at velocity u (one can also consider that c 1 is made of red molecules, c 2 of green molecules, in which case the combination c is yellow; to advect this yellow fluid, one can advect the red molecules with u and the green ones with u too). The presentations (2) or (25) correspond to one or the other of these interpretations: do we want to consider both fluids together, or do we treat them separately. For the numerical method, it consists in applying ELLAM on one and MMOC on the other.
The following definition summarises the combined ELLAM-MMOC scheme. satisfies 4.2. Analysis of mass balance error. Taking z 1 = K∈M z K (so that Π C z 1 = 1 in Ω) in (26) and plugging into (14), the mass balance error e mass of the ELLAM-MMOC scheme is estimated as follows: By using (22) and doing a change of variable F −t as in (23), we obtain Hence, the mass balance error e mass of the ELLAM-MMOC scheme is minimal when δt (n+ 1 2 ) tends to 0 or, setting U = {x ∈ Ω : divu (n+1) (x) = 0}, if (U ) (that is, pure ELLAM is used on the traceback of non-divergence free regions), or (that is, if MMOC is used -partially or entirely-on a domain D that is inside the traceback of non-divergence free regions, then the approximate concentration should almost be constant and stationary on D, and α should also be constant on D).
Assume that α is piecewise constant on M and only takes the values 0 and 1. Each cell M ∈ M can then be classified as M ELLAM (corresponding to α = 1) or M MMOC (corresponding to α = 0). The above relation is then re-written If K ∈ M ELLAM , then we only need to compute the integral of the first term on the right hand side of (27) since the latter terms will be zero. Otherwise, only the second and third terms are computed. These are approximated by taking the average value of f and divu (n+1) on the respective cells. We will demonstrate in Section 5 that, with a proper choice of α, the combined ELLAM-MMOC scheme can be implemented with an equivalent or cheaper computational cost than ELLAM, does not degrade much (or at all) the global mass conservation properties (contrary to MMOC), and has a much better local mass conservation compared to ELLAM on non-Cartesian meshes.

4.4.
Comparison with the MMOCAA. Of particular interest is a comparison with the MMOC scheme with adjusted advection (MMOCAA), first introduced in [11]. The MMOCAA is a modification of MMOC designed to conserve the global discrete mass. Simply stated, the modification consists of perturbing the foot of the characteristic F −δt (n+ 1 2 ) (x) by a term of order O((∆t) 2 ). Fixing η ∈ (0, 1), set For simplicity of notation, we denote the difference in mass for the MMOC scheme d mass to be We then define otherwise.
In order to enforce a discrete conservation of mass, the term Π C c(F , is equal to 0. For a more detailed presentation and implementation of the MMO-CAA, we refer the readers to [11,12,21]. It should be noticed that, contrary to the underlying principles of the ELLAM-MMOC scheme, there is no physical justification for using this parameter γ to enforce the mass conservation (that is, d mass = 0). Moreover, in some instances, at the first few time steps of a simulation, Ω φ(x)(Π C c (n) )(F −δt (n+ 1 2 ) (x))Π C z(x)dx = Ω φ(x) Π C c (n) (x)Π C z(x)dx and thus mass conservation cannot be achieved for any γ; see [12]. Also, in order to be able to determine the proper value for γ, one needs to evaluate both However, even though γ ensures that the MMOCAA achieves global mass balance, local mass conservation is not necessarily satisfied. Hence, after computing γ, adjustments, such as those at the core of our ELLAM-MMOC approach (see Section 2.5), need to be performed in order for the MMOCAA to provide numerical results that preserve local mass balance.

5.
Application: miscible flow model 5.1. Presentation of the model. We consider the model of miscible flow in porous medium represented by the following coupled system of elliptic and parabolic PDEs [22,16]: The unknowns are p(x, t), u(x, t), and c(x, t) which denote the pressure of the mixture, the Darcy velocity, and the concentration of the injected solvent, respectively. The functions q + and q − represent the injection and production wells respectively, and D(x, u) is the diffusion-dispersion tensor Here  (1) is the mobility ratio of the two fluids. As usually considered in numerical tests, we take no-flow boundary conditions: The concentration equation is completed by an initial condition, and the pressure equation by an average condition:

GDM-ELLAM-MMOC (GEM) scheme.
The idea is to implement a time-marching algorithm, wherein gradient discretisations (as described in Section 1.3) are used to approximate the diffusive terms for both (28a) and (28b). Some examples for which different GDs are applied for each equation in (28) have been presented in [9]. Note that the GDs used for the pressure equation do not need to involve the time components (time steps and interpolant of the initial condition). As highlighted in Section 4, combining the ELLAM and MMOC for the treatment of the advective terms removes the main disadvantages of each of the schemes. Hence, we propose to use the combined ELLAM-MMOC scheme for the advective component. We will refer to the combination of the GDM with the ELLAM-MMOC scheme for the complete coupled model (28) as the GDM-ELLAM-MMOC (GEM) scheme. Writing A = K µ , the following definition of the GEM scheme is inspired by the construction of the GDM-ELLAM scheme in [9,8] and by the design of the ELLAM-MMOC scheme for the advection-reaction model (Definition 4.2).
Definition 5.1 (GEM scheme). Let P = (X P , Π P , ∇ P ) be a space GD for the pressure, and C = (X C , Π C , ∇ C , I D , (t (n) ) n=0,...,N ) be a time-space GD for the concentration. Let α : Ω → [0, 1]. The GEM scheme for (28) reads as: find (p (n) ) n=1,...,N ∈ X N P and (c (n) ) n=0,...,N ∈ X N +1 C such that c (0) = I C c ini and, for all n = 0, . . . , N −1, q ± (·, s)ds (or, alternatively, q ± n = q ± (t (n) ) if q ± are continuous in time). ii) A Darcy velocity u (n+1) P is reconstructed from p (n+1) and, to account for the advection term in the concentration equation, the following advection equation is considered; it defines space-time test functions from chosen final values: iii) Using a weighted trapezoid rule with weight w ∈ [0, 1] for the time-integration of the source term and setting U where we recall the notations (10), and we set q ± N = q ± N −1 if these quantities are defined by averages on time intervals (there is no time interval (t (N ) , t (N +1) )).

Remark 5.2 (About the gradient discretisation and the velocity reconstruction).
The gradient discretisation P for the pressure equation (29) can be chosen so that the scheme is locally conservative. Some examples of such schemes are the mixed finite element (MFEM) and hybrid mimetic mixed (HMM) schemes. If the gradient discretisation gives a velocity field A(x, Π C c (n) )∇ P p (n+1) that is already in H(div, Ω), which is the case if the GD corresponds to an MFEM scheme, then the Darcy velocity u (n+1) P is taken equal to that field. For some other gradient discretisations, such as the ones corresponding to the HMM scheme, the velocity field is not in H(div, Ω) and a specific reconstruction, based on the numerical fluxes, must be made (see e.g. [8,Section 2.3]).
Key to an efficient and accurate implementation of the GEM scheme is a proper choice of α so that mass conservation is achieved without having to deal with the steep source terms encountered in ELLAM. Remarks 2.3 and 4.3 give us an idea of how to define the function α. In the context of the complete coupled model (28), the non-divergence free regions are the injection and production cells. Moreover, it is expected that, for an injection cell C + , (C − ). This indicates that for an efficient application of the GEM scheme, the MMOC component should be implemented on regions near the injection cells once the concentration c is almost constant in these regions. This happens after some time T + when the injection cells C + are almost filled up, i.e. c ≈ 1 in C + . Before this, we should implement a pure ELLAM scheme. Hence, we start by defining α = 1 over Ω for all n such that t (n) ≤ T + , the time where the injection cells are filled up; typically, T + ≈ 1 to 1.5 years. Note that T + can actually be found during the simulation, by checking if the concentration is almost constant in and around the injection cells or not. To be specific, T + is determined to be the time such that |Π C c (n) − Π C c (n+1) | < in the cells surrounding the injection well(s) and the well(s) themselves. For our numerical tests, we take = 10 −4 . For n such that t (n) > T + , and assuming for simplicity one injection cell C + and one production cell C − , a possible choice is Here, |x − C + | and |x − C − | denote the distance between x and the center of the cells C + and C − , respectively. This tells us to use ELLAM for regions far from the injection well, and MMOC otherwise. In particular, the choice of α = 1 near the injection cell and α = 0 away from the injection cell ensures that the second condition for mass conservation in Remark 4.3 is satisfied around the traceback region of C + . Moreover, for regions that are not in the traceback of C + , we have that divu = 0, and hence mass conservation is attained by MMOC c.f. Remarks 3.2 and 3.3. In case of multiple injection and production wells, the same rule can be applied by taking α(x) = 1 if the closest well to x is a production well, and α(x) = 0 if the closest well to x is an injection well.

5.3.
Adaptation of local mass conserving adjustments to the miscible flow model. In this section, we write the algorithm proposed in Section 2.5 in the context of the miscible flow model, for the GDM-ELLAM and GEM scheme. Since the GDM-MMOC does not achieve a global mass balance, it does not make sense to post-process the data and enforce local mass balance for each cell in the mesh. For the miscible flow model (28), the velocity field is divergence free in most regions, except for the injection and production wells (see (28a)). Hence, for most cells K ∈ M, we know that |F 2 ) (K)| cannot be determined exactly are: the production cells, and the cells K = C + that track into the injection cells (fully or partially). We thus formulate an approximate local volume constraint for the cells K = C + that track into the injection cells. Due to Liouville's Theorem, the exact volume of the traceback region of an injection cell C + is deduced to be |F We compute an approximate trace-forward region C + of C + by forming a polygon with vertices and edge points of C + tracked forward in time. In practice, compared to the other cells K in the mesh, we track more points along the edges of C + in order to obtain a good enough approximation of C + . In particular, if, on average, n points are tracked along the edges of a cell K ∈ M, then we found that 4n + 1 is an appropriate number of points to track along the edges of C + . We then set the approximate local volume constraint for the cells K = C + that track into C + to be Essentially, (33) may be interpreted in the following manner: the region K ∩ C + is the part of K that is expected to track back into the injection cell C + . Hence, since the other part of K stays in the divergence free region, its volume |K| φ − |K ∩ C + | φ remains unchanged. The changed volume is then approximated as a ratio of the part in C + that will be tracked out of itself (i.e. (1 − e −β )|C + | φ ). We can then adapt the algorithm proposed in Section 2.5 to the miscible flow model. We note that the local approximation (33) is valid only if we have a good approximation of the trace-forward region of the injection cells, and hence we must track more points on the injection cells. As an additional fix, since we expect the injection cells to be eventually filled with the injected fluid, we set c = 1 on these cells. Assuming that (33) gives a good enough approximation, the local volume constraint for the production cells should then be satisfied approximately (since we have global mass conservation). 2 ) (C + )| φ should be very close to their actual values. Here, it is acceptable to track more points than what we track on average for the other cells K ∈ M to ensure that such an accuracy is obtained, since the global impact in terms of cost is minimal. In particular, one clear advantage of the GEM scheme comes from the absence of the approximations (33), which tells us that GEM is expected to give a better local mass conservation property compared to the GDM-ELLAM.

Remark 5.3 (Convergence tests).
For this miscible displacement flow model, very few analytical solutions exist, and those that are known have a velocity field u such that div(u) is non-zero almost everywhere in the domain (see, e.g., [25, Section 3.1]), which means that the wells (corresponding to q + and q − ) are also completely diffused in the domain. Also, when implementing an ELLAM scheme for these test cases, we do not encounter steep back-tracked functions as in Remark 2.3. As a consequence, the most efficient implementation of the GEM scheme boils down to taking α = 1 everywhere in (31), corresponding to a pure ELLAM scheme. Due to this, we do not provide convergence tests for the GEM scheme, as it will yield the same results as ELLAM, and would not give any interesting insight on the efficiency of the GEM scheme in preserving global and local mass balances -its very raison d'être.

Numerical Results.
In this section, we compare numerical results obtained from the GEM scheme to those obtained from GDM-ELLAM. Both of these schemes employ the post-processing technique outlined in Section 5.3 to achieve local mass balance. For completeness, we will also present a comparison with GDM-MMOC, but without the local adjustments. As with the GEM scheme, for the GDM-MMOC, an ELLAM scheme is first implemented for the first few time steps, when t (n) ≤ T + , after which, a pure MMOC scheme is implemented, i.e. for t (n) > T + , we take α = 0 over Ω in (31). Here, we use the HMM scheme [15] for discretising the diffusive terms. This HMM corresponds to a certain choice of the gradient discretisation C and P, see [14,Chapter 13]. The numerical simulations are performed under the following standard data (see, e.g., [27]): (1) Ω = (0, 1000) × (0, 1000) ft 2 , (2) injection well at (1000, 1000) and production well at (0, 0), both with flow rate of 30ft 2 /day,   Figure 5. This is followed by a solution obtained by a HMM-GEM scheme (GEM scheme using the HMM gradient discretisation, with advection components computed as described in (27)) in Figure  6. Also, for this test case, we present a numerical result obtained from an MFEM-ELLAM scheme. This numerical output from MFEM-ELLAM was obtained by a straight application of the MFEM-ELLAM algorithm as presented in [27], with several hundred of quadrature points per cell around the injection well, without any post-processing or tweaks. One of the authors of [27] has informed us that a post-processing was implemented near the injection wells in order to obtain a better concentration profile than those we find in Figure 6, right. However, due to issues on intellectual property rights, this post-processing technique was not revealed, and hence our tests were not able to account for this. These are accompanied by Table  1, which presents some important features, such as the number of points tracked along each edge, overshoots/undershoots, e (N ) mass , and the approximate amount of oil 1 |Ω| φ Ω φΠ C c (N ) recovered after 10 years. Here, e (N ) mass refers to the accumulated mass balance error (percentage) obtained over all the time steps, i.e. .
In practice, we have e   Upon comparing the concentration profiles, we see that for all of the HMM based schemes, the overshoot is very low, with the maximum overshoot being less than 0.2%. However, it can be noted in Table 1 that ELLAM's 0.18% overshoot is much larger, by a factor of almost 20, than those of the MMOC and GEM, even on a very simple square mesh. For the MFEM-ELLAM, the overshoots were at 14.17%, which is much larger than those observed from the other schemes. Also, we note here the presence of severe undershoots, of around 21%, near the diagonal. Aside from the overshoots, there are no noticeable differences between the concentration profiles for the HMM-ELLAM and the HMM-GEM. On the other hand, we note that the HMM-MMOC scheme introduces some artificial diffusion along the diagonal, which slightly smears the expected fingering effect. Next, upon comparing the approximate amount of oil recovered after 10 years, the 68.44% to 69.14% obtained for the HMM-GEM scheme is comparable to the amount from the HMM-ELLAM, which ranges from 69.76% to 70.09%. The HMM-MMOC, on the other hand, provides an overestimate of the oil recovered when the tracked cells are approximated only by vertices and edge midpoints, due to the excess diffusion it introduces along the diagonal becoming more prominent. When 3 points are tracked along each edge, the amount of oil recovered for HMM-MMOC is almost the same as that for HMM-ELLAM and HMM-GEM. We note that the approximate amount of oil recovered by MFEM-ELLAM is only at 57.39%, which is possibly caused by the undershoots on the region near the diagonal.
Lastly, we compare the mass balance errors. Since the recovery for MFEM-ELLAM is only at 57.39%, we expect the mass balance errors to be quite large. Now, we focus on the mass balance errors for the HMM based schemes, obtained once we track 3 or more points along each edge. The error obtained from the GEM (0.85%) is much better than the one from MMOC (2.80%), and close to that obtained from ELLAM (0.21%). These results agree with the analysis provided in Sections 3.4 and 4.2, due to the fact that the MMOC will fail to conserve mass as soon as the fluid starts invading the production well (which translates to

Hexahedral meshes.
We then compare the numerical results on hexahedral meshes. Unlike the regular square cells for the Cartesian type meshes, the cells for the hexahedral meshes are irregular. To implement the HMM-ELLAM scheme on these types of cells, the proper amount of points to track along each edge is determined by the cell regularity parameter [8], defined for each cell K ∈ M to be log 2 (m Kreg ) points are then tracked along each edge of cell K. For hexahedral and Kershaw meshes, we have max K∈M ( log 2 (m Kreg ) ) = 3 and 6, respectively. We note however, that due to the adjustments described in Section 5.3 which need to be implemented in order to achieve local mass conservation, 2 log 2 (m Kreg ) + 1 points would need to be tracked along the edge of each cell for hexahedral type meshes. In the following figures, three tests are performed for the HMM-ELLAM: the first of which does not involve any adjustment to achieve local mass conservation, followed by an adjustment only on the cells that are not involved with either the injection or production cells (i.e. the algorithm in Section 5.3 without (33)), and finally an adjustment based on the full algorithm in Section 5.3.  Upon looking at the results in Figures 7-9 and Table 2, it is noticeable that the solution from HMM-ELLAM has a large discrepancy and overshoot near the injection well. This is due to the fact that the algorithm in Section 5.3 fails to converge for hexahedral meshes. We note here that such a failure of the adjustment algorithm was not noticed on Cartesian meshes, whether in the tests conducted above or (with a different approach to the adjustment) in [4]. Our tests on hexahedral meshes demonstrate here the difficulty of designing a robust algorithm to locally adjust the mass balance for ELLAM. As a matter of fact, it has been pointed out in [10] that there is no guarantee that adjustments for ELLAM type schemes, such as those in [4], will terminate or yield a valid mesh configuration. The same issue  could happen to our proposed adjustment applied on the GEM method but, as seen in the tests above, this adjustment seems to be more robust. There are two possible explanations for why the algorithm in Section 5.3 fails to converge for hexahedral meshes: Firstly, compared to the Cartesian mesh, the volume of the injection cell, and each of the cells around it (around 700 to 1000 square units), is much smaller than the volume of the other cells in the mesh (on average, 2000 square units). Since these cells are already small to begin with, tracking them backwards will lead to traceback regions which are much smaller, and hence will be more prone to errors. Taking ∆t = 36 days, the expected volume of the traceback region of the injection cell e −α |C + | ≈ 1e − 04 square units, which is only around 1e-08% of the total volume. This gives us an idea that taking a smaller time step might be able to mitigate these errors; however, even taking a much smaller time step of ∆t = 1 day, there is still no convergence. Due to this, we conclude that it is caused by the second possibility, i.e. for hexahedral meshes, (33) does not give a good approximation of the local volume constraint for the cells tracked back into the injection cell. This can be seen more clearly upon comparing Figure 8 (left and right). Without enforcing (33), the solution seems to behave better. Actually, upon implementation of (33), the local volume constraint is satisfied on the cells which are not tracked back either into injection or into production cells. Having eliminated the local mass balance errors from these cells, they accumulate onto the cells near the injection cell. Since (33) does not give a good approximation, the accumulated error around this region does not spread and cancel out properly, and hence severely distorts the quality of the numerical solution. Contrary to the HMM-ELLAM, due to the absence of (33), this problem is not as severe with the GEM scheme, and can be solved by taking a smaller time step of ∆t = 18.
With the exception of the case for which ELLAM is adjusted with the local volume constraint (33), the amount of oil recovered from all three schemes are comparable, as they are within 2% of each other. Upon comparing the mass balance errors, we note that the HMM-GEM (0.74%) outperforms the HMM-MMOC (1.82%), and is in the same range as the best HMM-ELLAM implementation (0.31%). This example shows that, on non-Cartesian meshes, due to the absence of (33), the HMM-GEM is able to provide a better-looking solution, with reduced overshoots and acceptable mass conservation properties, compared to the HMM-ELLAM method. Moreover, the HMM-GEM conserves mass better than the HMM-MMOC.

5.4.3.
Kershaw meshes. Finally, we compare the numerical results on a much more challenging mesh, the Kershaw mesh.  As with the Cartesian mesh test case, no significant difference can be observed between the numerical solutions obtained from HMM-ELLAM and HMM-GEM. Also, as expected, the mass balance error for the HMM-MMOC is quite large. Notably, the numerical solution on Kershaw type meshes is skewed towards the lower right corner. This is expected due to the fact that the numerical fluxes for HMM schemes on this type of mesh are prone to grid effects, as explained in [8].
To summarise the previous tests, the HMM-GEM exhibits a slightly better volume conservation property than the HMM-ELLAM. This is due to the local volume constraint for the HMM-ELLAM being inexact in the sense that it depended on (33) to give a good approximation; whereas for the HMM-GEM scheme, exact local volume constraints were imposed. In particular, on non-Cartesian meshes with small or mildly distorted cells, HMM-GEM can control local volume constraints more easily and more accurately than HMM-ELLAM, as was seen in the test case on hexahedral meshes. Hence, with the choice of α driven by the discussion in Section 4.2, GEM achieves both a good preservation of the physical bounds on c, and of mass conservation.

Convergence result for GDM-MMOC scheme
In this section, we present a convergence result for the GDM-MMOC scheme. To simplify the exposition, this convergence is stated for the scalar advectionreaction-diffusion model (1); at the expense of heavier notations, a convergence result of GDM-MMOC could be stated for the coupled model (28), as in [9] for the GDM-ELLAM method. This convergence result is established using only weak regularity assumptions on the solution (which are satisfied in practical applications), and not relying on L ∞ bounds (which are impossible to ensure at the discrete level given the anisotropic diffusion tensors and the general grids used in applications). Hence, the convergence is established under compactness arguments, and error estimates are not available. As a matter of fact, the convergence analysis of numerical approximations of (28) under weak regularity assumptions has recently received an increasing interest; see, e.g., [6,7] for finite volume methods and [23,18] for discontinuous Galerkin methods. We also note that this convergence analysis assumes a perfect computation of the tracked regions. An analysis which accounts for the approximation in tracked regions, and adjustment strategies, has been performed in [3], under the assumption that u ∈ C 1 (Q T ) d . Future work will aim to perform a convergence analysis, without this strong regularity assumption.

Λ(x, s)ds
If c = (c (n) ) n=0,...,N is a solution to the GDM-MMOC scheme, we define Π C c ∈ L ∞ (Q T ) and ∇ C c ∈ L ∞ (Q T ) d by ∀n = 0, . . . , N − 1 , ∀t ∈ (t (n) , t (n+1) ) , for a.e. x ∈ Ω, The convergence of the GDM-MMOC can be established under most of the assumptions made for the GDM-ELLAM (i.e. [9, Assumptions (A1),(A3)-(A4)]). Assumption (B1) makes precise the sense in which the velocities used to construct the flows over each time step (see (8)) approximate the given velocity u.  ) and all n = 0, . . . , N m − 1, we have, as m → ∞, u m → u weakly in L 2 (Q T ) d . We now define the notion of weak solution to the advection-reaction-diffusion model, and state the convergence result for the GDM-MMOC. The proof of this convergence, that relies on discrete compactness techniques, is similar to the proof of [9, Theorem 3.3] and is therefore omitted. • Π Cm c m → c weakly- * in L ∞ (0, T ; L 2 (Ω)) and strongly in L r (0, T ; L 2 (Ω)) for all r < ∞, • ∇ Cm c m → ∇c weakly in L 2 (Q T ) d , where c is a weak solution of (1).