A new efficient Eulerian–Lagrangian localized adjoint method for solving the advection–dispersion equation on unstructured meshes

Advances in Water Resources 29 (2006) 1056–1074 www.elsevier.com/locate/advwatres

A new efficient Eulerian–Lagrangian localized adjoint method for solving the advection–dispersion equation on unstructured meshes A. Younes *, P. Ackerer, F. Lehmann Institut de Me´canique des Fluides et des Solides, Universite´ Louis Pasteur—CNRS—UMR 7507, 2 rue Boussingault, 67000 Strasbourg, France Received 20 October 2004; received in revised form 7 September 2005; accepted 7 September 2005 Available online 2 November 2005

Abstract A new and high efficient scheme is developed for the Eulerian–Lagrangian Localized Adjoint Method (ELLAM) to solve the Advection–Dispersion transport Equation (ADE) on unstructured triangular meshes. To obtain accurate results, the new method requires a very limited number of integration points (usually 1 per element). The scheme uses only strategic points as numerical integration points. With this scheme, locations of integration points and weights are assigned at the new time level and then backtracked to the old time level without redistributing the weights. Interpolation problems are minimized since we use continuous characteristics and only changes due to dispersion are interpolated to obtain the concentration at the foot of each characteristic. Different numerical experiments with a large range of grid Peclet numbers are presented to compare the new ELLAM to the standard one and to the Discontinuous Finite Element Method. The new ELLAM gives more accurate results and is much less CPU time consuming than all other methods especially with large time steps and highly unstructured meshes.  2005 Elsevier Ltd. All rights reserved. Keywords: ELLAM; ADE; Tracking; Quadrature; Strategic integration points; Unstructured meshes

vi vj jVj

1. Introduction

Dij ¼ ðaT jVj þ Dm Þdij þ ðaL  aT Þ

The most common mathematical model to simulate tracer transport in porous media is an advection–dispersion partial differential equation which has the following formulation:

L(C) is the differential operator defined above, V [L/T] is a given divergence free fluid velocity of components vi, aL and aT [L] are the longitudinal and transverse dispersivities, dij is the Kronecker delta function, and Dm [L2/ T] the molecular diffusion coefficient. T is the end of the simulation time period starting at time zero. Eq. (1) is subject to the initial and boundary conditions

oC þ r  ðVCÞ  r  ðD  rCÞ ¼ 0 ot x 2 X; t 2 ½0; T 

LðCÞ ¼

where C(x, t) [M/L3] is the unknown concentration at location x and time t and D [L2/T] the dispersion tensor defined by

Cðx; 0Þ ¼ C 0 ðxÞ Cðx; tÞ ¼ g1 ðx; tÞ

Corresponding author. E-mail address: [email protected] (A. Younes).

0309-1708/$ - see front matter  2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2005.09.003

ð2Þ

x2X ðx 2 oX1 ; t > 0Þ

ðD  rCÞ  noX ¼ g2 ðx; tÞ *

i; j ¼ 1; 2

ðx 2 oX2 ; t > 0Þ

ðVC  D  rCÞ  goX ¼ g3 ðx; tÞ ðx 2 oX3 ; t > 0Þ

ð3Þ

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

where X is a bounded, polygonal open set of R2, oX1, oX2 and oX3 are partitions of the boundary oX of X corresponding to Dirichlet, Neumann and total flux boundary conditions and goX the unit outward normal to the boundary oX. In many field situations, especially for small scale simulations, the hyperbolic advective part of the partial differential equation becomes dominant. In this case, standard numerical methods, such as finite elements or finite volumes, generate solution with numerical diffusion and/or non-physical oscillations, which can be avoided with a very fine discretization in both space and time. This can make standard numerical methods very expensive and therefore not suited for advection dominated transport problems. The Eulerian–Lagrangian Localized Adjoint Method (ELLAM), introduced by Celia et al. [2], is an interesting alternative to standard methods to solve (1). The ELLAM uses space–time test functions. It preserves the performance of characteristic methods and treats general boundary conditions naturally in their formulation. The ELLAM has been extended to multi-dimensions problems with varying test and trial functions [1,5,12,13,15]. A review of the researches done on the ELLAM is presented by Russell and Celia [13]. Almost all developed ELLAMs suffer from non-physical oscillations. The first reason of oscillations is the use of the consistent mass matrix which is non-diagonal. Mass lumping can be used to control oscillations but is known to add numerical diffusion for Eulerian– Lagrangian methods [10]. This phenomenon is avoided for the one-dimensional problem when ELLAM is combined with a moving grid procedure [16]. This approach cannot be extended to 2 or 3 dimensions. Russell and Binning [10] proposed a selective lumping scheme to avoid excessive numerical diffusion. The basic idea is to add the minimum numerical diffusion required to avoid oscillations. This method cannot be easily generalized for multidimensional unstructured meshes. The second reason of oscillations is the use of the forward tracking approach. With this approach, quadrature points and weights are defined on the old time level. Quadrature weights will not necessarily sum up to the correct volume for each mesh element, especially for unstructured meshes, on the new time level. To minimize this phenomenon, Healy and Russell [4,5] use Strategic Space Integration Points (SSIP) which are backtracked from the new time level to be included in the quadrature on the old time level. It is also known that Eulerian–Lagrangian methods perform well for problems in which they can successfully use a large time step, but, when using several time steps, they suffer from numerical dispersion introduced by interpolation at each time step [11]. On the other hand, in the literature, numerical results of ELLAM are usually presented on uniform rectangu-

1057

lar or block meshes. However, for practical problems, triangular meshes are preferred since they can handle domains with complex geometry and local mesh refinement. The ELLAM was used for the first time for triangular meshes in [17]. In this paper, the authors studied different spatial and temporal numerical integration techniques to evaluate boundary and spatial integrals with ELLAM. They suggest a new algorithm in order to avoid excessive numerical diffusion when using many time steps. The basic idea of this approach is to keep the same characteristics for all time steps and to interpolate only the concentration variations due to the dispersion process at the end of each time step. The effect of the mesh was also studied in [17]. It was shown that, when working with unstructured triangular meshes, the number of integration points has to be increased in order to avoid non-physical oscillations. Therefore, ELLAM becomes less efficient since it requires a lot of integration points for highly unstructured meshes. Moreover, the required number of integration points cannot be known a priori and is problem-dependent which can make ELLAM less attractive for practical cases. Therefore, the main objective of this work is to describe a new and efficient algorithm to solve the ADE with the ELLAM with a very limited number of integration points for unstructured triangular meshes. The new algorithm is shown to improve the methodÕs computational efficiency without non-physical oscillations. In this scheme, cumulative interpolation errors are very limited by continuously tracking integration points as in [17]. The article is outlined as follows. In the first part of this paper, we recall briefly the mathematical developments of the ELLAM for solving the two-dimensional ADE on triangular meshes. Numerical integration and improved implementation methods are described in the second part of the paper. Their accuracy and performance are compared in the third part of the paper. Finally, the robustness of the suggested new method is evaluated using different Peclet numbers and different time and space discretizations. The new ELLAM scheme is compared to a standard formulation (see [17]) and to another numerical model based on a combination of discontinuous Galerkin/mixte finite element in order to evaluate the efficiency of the developed ELLAM for unstructured meshes and grid Peclet numbers ranking from 1 to 1000.

2. The Eulerian–Lagrangian localized adjoint method for ADE 2.1. ELLAM formulation The weak formulation of Eq. (1), using space–time test function x = x(x, t) leads to

1058

Z

T

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

Z

Z

LðCÞx X 0  Z TZ  oC x þ r  ðVCÞx  r  ðD  rCÞx dxdt ¼ 0 ¼ X ot 0

Z

T

E Z X

¼

The ELLAM selects x(x, t) so that for t 2 [t , t Dx ox ¼ þ V  rx ¼ 0 Dt ot

n+1

], ð5Þ

x(x, t) is discontinuous at t = tn, and x(x, t) = 0 for t 62 [tn, tn+1]. With this choice of x, the ADE to solve becomes (for details, see [1]) Z Z T Z nþ1 nþ1 C x ðxÞdx þ ðD  rCÞ  rx dx dt X

þ

Z

T

0

¼

Z

Z

0

oX

C n xn ðxÞdx

ð6Þ

X

The test functions x is constant along the characteristic since it verifies (5). At t = tn+1, we choose x as the piecewise linear function (chapeau function). ~; ~tÞ be the characteristic passing Let xðtÞ ¼ X ðt; x through a given point ð~ x; ~tÞ, with ~t 2 ½tn ; tnþ1 , and determined by the following system of equations:  dxðtÞ=dt ¼ VðxðtÞ; tÞ ð7Þ ~ xð~tÞ ¼ x

2.2.1. Total flux inflow boundary integral Let E be the number of edges e, of length le, for which nþ1 the total flux TQe is prescribed and constant for t 2 [tn, tn+1]. The boundary integral becomes

xi ðx; tÞdl dt

ð8Þ

p¼1 q¼1

To evaluate xi(xp, tq), the integration point xp at time tq is tracked forward to xp at tn+1 and we use xi ðxp ; tq Þ ¼ xi ðxp ; tnþ1 Þ. The contribution of this term is placed in the right-hand side vector of the final system to solve. 2.2.2. Dirichlet inflow boundary integral Let E be the number of edges e, of length le, for which the concentration is prescribed for t 2 [tn, tn+1]. The boundary integral becomes Z

Z

T

½ðVC  D  rCÞxi noX dl dt

oX

0

¼

Z

Z

T

0



Z 0

We use the mass lumping procedure for the storage term and evaluate the dispersion integral using a onestep backward Euler approximation in time. The two first integrals of Eq. (6) are therefore solved with standard finite element method, and are not developed here (see [17] for details). The whole difficulty of ELLAM lies in the accurate approximation of the third term, the boundary integral, and the fourth term at the right hand side of (6) which corresponds to the mass at the old time level. Recall that, in this work, we consider only the case of divergence free velocity without sources and sinks. Some care can be required to approximate strong point sources or sinks in ELLAM discretizations as stated in [13].

Z le

e¼1

This notation can refer to tracking forward or backward in time; so that (x, tn+1) backtracks to (x*, tn) and (x, tn) tracks forward to ðx; tnþ1 Þ. 2.2. Evaluation of the integrals

nþ1

TQe

Numerical integration over the edge e is used to evaluate the previous integral with Ns one-dimensional spatial integration points, along le, with weights W pe and Nt integration points in time with weights W qe . Hence, Z T Z ½ðVC  D  rCÞxi noX dl dt oX 0 ( ) Ns Nt E X  nþ1 X X  p q p q TQe xi ðx ; t ÞW e W e  ð9Þ

X

½ðVC  D  rCÞx  noX dx dt

tnþ1

tn

e¼1

ð4Þ n

½ðVC  D  rCÞxi noX dl dt

oX

0

½ðVCÞxi noX dl dt

oX T

Z

½ðD  rCÞxi noX dl dt

ð10Þ

oX

nþ1

Let C e be the prescribed concentration for the edge e nþ1 and Qe the water flux through e. As previously, the numerical integration leads to Z

Z

T

½ðVCÞxi noX dldt 

oX

0



E X

( nþ1 nþ1 Qe C e

e¼1



Ns X Nt E X X

Z

T 0

Ns X Nt X

Z

½ðD  rCÞ  noX xi dldt

oX p

) q

ðxi ðx ; t ÞW

p q eW eÞ

p¼1 q¼1

n o tnþ1 ne  ½ðD  rCÞxi xp W pe W qe

ð11Þ

e¼1 p¼1 q¼1

The first term is treated like the total flux boundary term. xi(xp, tq) is evaluated using xi ðxp ; tq Þ ¼ xi ðxp ; tnþ1 Þ and the contribution of this term is placed at the right-hand side of the final system. The second term is unknown and is approximated as in Binning and Celia [1] by forward-tracking the boundary values to the new time level. This term introduces an asymmetry into the matrix.

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

1059

40

30

20

10

0 0

10

20

30

40

50

60

70

80

90

100

Fig. 1. The domain X discretized with non-uniform triangular elements.

2.2.3. Outflow boundary condition In the literature, the outflow boundary is often treated using additional degrees of freedom depending on the Courant number (see [13] for example). In practice, we have almost exclusively a zero diffusive flux at the outflow boundary. In the purely advective case, the mass flowing out of the domain comes from the points at the old time level that are tracked to the outflow boundary during a time interval. Therefore, no additional degrees of freedom are required if the mass flowing out of the domain is not considered. In practice, all spatial integration points which reach the outflow boundary at a time less than Dt when tracked forward, will have zero as integration weighting factor. 2.2.4. Integral at the old time level (fourth term of (6)) A numerical integration is used to evaluate this term " # Z NA X X p n n n p n p n C xi dx ¼ ðC ðx ; t Þxi ðx ; t ÞW A Þ ð12Þ X

A

p¼1

NA integration points, with their corresponding 2D space weights W pA , are used for each element A. These points are located at xp at time tn and are tracked forward to xp at tn+1 to evaluate xi ðxp ; tn Þ ¼ xi ðxp ; tnþ1 Þ. Cn(xp, tn) in (12) is evaluated by linear interpolation of the known solution at nodes of the finite element mesh.

3. Implementation of the ELLAM In this part, we recall the standard approach used in [17] for numerical integration on unstructured triangular meshes. To improve the computational efficiency of ELLAM, a new approach is proposed to reduce the number of required integration points and to avoid non-physical oscillations. The new numerical integration scheme is used to approximate the boundary integral (11) and the mass

integral at the old time level (12). This scheme is described here and compared to a more usual scheme given in [17] using a standard test problem. 3.1. Test problem definition ELLAM implementations will be verified by solving a transport problem in the spatial domain X of rectangular shape (0, 100) · (0, 40) which is partitioned into triangular elements (Fig. 1). Time is divided into discrete intervals Dt = tn+1  tn. Fluxes at each edge are obtained by solving the stationary flow problem, solution of 8 DP ¼ 0 > > > > >

in X in oX1L

> P ¼0 in oX1R > > > > : ðrP Þ  goX ¼ 0 in oX2

ð13Þ

where oX1L and oX1R are, respectively, the left and the right hand perpendicular sides of the domain and oX2 = oX/oX1. Accurate fluxes at each element edge are required to track the particles. Therefore, the flow equation (13) is solved with the mixed finite element approximation detailed in [18]. With the mixed formulation, the velocity is defined everywhere in the field and is continuous across the inter-element boundaries. Integration over each element gives the corresponding streamlines using the technique developed by Pollok [8]. The boundary conditions for the transport problem are of Dirichlet type at the inflow (left boundary of X), with 8 > < C ¼ 0 for x ¼ 0 and 0 6 y < 12 > :

C ¼ 1 for x ¼ 0 and 12 6 y 6 28 C ¼ 0 for x ¼ 0 and 28 < y 6 40

ð14Þ

1060

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

A zero diffusive flux is imposed at the outflow boundary. The corresponding analytical solution for an infinite domain is [6] ( " # Z T x y  12 3=2 Cðx; y; tÞ ¼ s erf 1=2 1=2 ð16paL Þ ð4aT sÞ 0 " #) " # 2 28  y ðx  sÞ þ erf exp  ds ð15Þ 1=2 4aL s ð4aT sÞ Rx with erfðxÞ ¼ p2ffiffip 0 expðs2 Þds. Solution of (13) leads to fluid velocity components equal to vx = 1.0 and vy = 0.0 over the whole domain. 3.2. Numerical integration of the boundary integral (11) The test problem is simulated with a single time step of 40 days. Since this problem corresponds to a constant injection into a domain originally free of contaminant, the spatial integration at the old time level (12) is not needed. 3.2.1. The standard scheme In the standard approach, the boundary integral (11) over each edge e, is approximated with a uniform distribution of numerical integration points (Fig. 2). The number of integration points is Ns · Nt where Ns and Nt are numbers of spatial and time integration points. Each integration point has a weight of W pe  W qe and is tracked forward in order to evaluate xi ðxp ; tq Þ ¼ xi ðxp ; tnþ1 Þ. Non-physical oscillations can appear in the solution since quadrature weights will not necessarily sum up to the correct volume for each mesh element. Several numerical integration methods are tested in [17] for the evaluation of the boundary integral (11) on unstructured triangular meshes. Strategic Time Integration Points (STIP) obtained from Strategic Spatial Integration Points (SSIP) when tracked to the boundary as defined by Healy and Russell [5], are used in [17]. The location and time of their arrival at the boundary become intermediate integration points for space and time, respectively. The trapezoid rule is used to obtain each weight of the integration points.

It is shown in [17] that the most accurate results are obtained when the number of time integration points Nt is depending on the number of spatial integration points Ns and the Courant number. Since we can have non-uniform meshes, we define CFLmax the greatest CFL over all elements for which a particle located at the center crosses the edge e in a time less than T when tracked backward. Nt is defined by Nt = 2 · CFLmax · Ns. SimpsonÕs rule is used for the numerical integration in time. This scheme was found to be accurate for unstructured meshes in [17] and will represent the standard scheme in the rest of the paper. 3.2.2. The new scheme Our aim is to develop a new algorithm with a limited number of integration points that leads to accurate results even for unstructured meshes. In order to avoid a uniform dense distribution of integration points, we use strategic integration points of Healy and Russell [4,5]. Integration points are located strategically at the later time level and then backtracked to the boundary (Fig. 3). Since we use only these strategic points as numerical integration points, weights are obtained via the compound trapezoidal rule as in [5]. Hence, our objective is to obtain accurate weights for these strategic points. In the case of zero dispersion, Eq. (6) written for the first time step leads to Z Z T Z C 1 x1 ðxÞdx þ ½ðVCÞx  noX dx dt ¼ 0 ð16Þ 0

X

oX

From (11), the boundary integral can be approximated by Z T Z nbp X ½ðVCÞxi noX dl dt  C j ðxp ; tq Þxi ðxp ; tq ÞW bound j 0

oX

j¼1

ð17Þ ¼ Qe W pe W qe and nbp is the total number of where W bound j boundary integration points. Since the domain is initially free of contaminant, the total mass in the domain at the end of the first time step is given by (17).

Δy

Δy

Δt

Δt

Fig. 2. Numerical integration of the boundary integral using a uniform distribution of integration points (Ns · Nt = 5 · 3).

Fig. 3. Numerical integration of the boundary integral using strategic integration points defined at the new time level.

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

On the other hand, for the second time step, the initial mass in the domain is approximated using (12) by Z nsp X C 1 x1 ðxÞdx  C j ðxp ; tn Þxi ðxp ; tn ÞW space ð18Þ j X

j¼1 space j

where W ¼ W pA is the 2D spatial weight and nsp is the total number of spatial integration points. Since dispersion is zero, the mass arriving from the boundary (Eq. (17)) is equal to the mass distributed over the domain (Eq. (18)). Moreover, if we use the same integration points (continuous characteristics for the two time steps) and since the concentration and the test function values are constant along characteristics, mass conservation leads to W bound ¼ Qe W pe W qe ¼ W space j j

ð19Þ

Therefore, the spatial–temporal integral at the boundary is converted to a two-dimensional spatial integral using (19). Then, weights of strategic points are spatial weights at the new time level. Hence, with the new scheme, we use P particles of weights W pA at the strategic locations over each triangle A (if P = 1, W pA corresponds to the area of A). Each particle starts at t = T and is tracked backward to t = 0. All particles j that cross the edge e (of prescribed flux or concentration) in a time tj > 0 will be used as boundary integration points. Using (19) for the conversion of boundary weights to spatial weights at the new time level, Eqs. (9) and (11) are evaluated with only strategic integration points Z T Z ½ðVC  D  rCÞxi noX dl dt oX 0 ! nþ1 nbp X TQe p q xi ðx ; t Þ ð20Þ  W space j nþ1 Qe j¼1 and Z

T 0

Z

½ðVCÞxi noX dldt 

oX nbp X



Z

T 0

Z

½ðD  rCÞ  noX xi dldt

oX

C j ðxp ; tq Þxi ðxp ; tq ÞW space j

j¼1

þ

nbp X 1 nþ1 j¼1 Qe

tnþ1

ne  ½ðD  rCÞxi xp W space j

ð21Þ

When working with only one time step (i.e. Dt = T), the algorithm of the new scheme is as following: • For each element, P particles are set at the strategic point locations xj at tj = Dt with their corresponding 2D spatial weights ðW space Þ. Weights are the area j corresponding to each particle. Indeed, if P = 1, xj corresponds to the location of the center of the element and W space corresponds to its area. If P = 4, j

1061

the triangular element is subdivided into 4 equal sub-triangles (by joining the edge midpoints), xj corresponds to the location of the center of each sub-triangle and W space corresponds to the area of the subj triangle. • Each particle is backtracked during Dt, from its position xj (at tj = Dt) to the departure position xj (at t ¼ tj Þ. xj can be located inside the domain if tj ¼ 0 or at the inflow boundary if tj > 0. • The boundary integral given by (20) or (21), is evaluated using only particles with tj > 0 and non-zero prescribed concentration at xj . Weights W space of strategic integration points are j defined at the new time level and can be obtained from quadrature rules including the well known Gauss scheme. Morton et al. [7] have shown that in most cases, quadrature rules lead to unstable schemes. An alternative area-weighting technique was presented in [7] for rectangular elements. The centroid of each element is tracked and the whole element deemed to move without distortion and rotation. The method was extended to triangular elements in Priestley [9]. The vertices of the triangle are first transported and the interior of the triangle is assumed to move in a linear manner. The method uses an exact projection from one arbitrary triangular grid to another. With the developed method, no mapping or projection is used. W space can be obtained with any quadrature j formula, since we lump the mass matrix. Mass lumping is known to stabilize the scheme but to add significant numerical diffusion. This phenomenon is minimized in our case by using a continuous tracking of characteristics (see next paragraph). ELLAM provides solutions with much more oscillations for unstructured meshes than for uniform meshes. In [12], flow orientated with 45 to the grid leads to results with non-physical oscillations. This phenomenon can be reduced by using finer grid [12] or increasing the number of boundary integration points. The mesh orientation effect is studied here by solving the previous test problem with three different triangulations: • a uniform triangulation (a) with 1071 nodes and 2000 elements (Fig. 4); • a structured triangulation (b) with 1107 nodes and 2076 elements (Fig. 1); • an unstructured triangulation (c) with 2071 nodes and 4004 elements (Fig. 4). For this last triangulation, the mesh refinement is located close to transverse concentration fronts in order to increase the sensitivity of the spatial resolution. The dispersion parameters are aL = 103, aT = 104 and Dm = 105.

1062

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074 40

30

20

10

0 0

10

20

30

40

50

60

70

80

90

100

0

10

20

30

40

50

60

70

80

90

100

40

30

20

10

0

Fig. 4. The domain X discretized with uniform and unstructured meshes.

Simulations are performed with both standard and new schemes. The standard scheme is used with different number of boundary integration points Ns = 3, Ns = 7, Ns = 15 and Nt = 2 · CFLmax · Ns. The new scheme is used with P = 1 and P = 4. Results of simulations are given in Table 1. This table gives the maximum of the concentration, the total number of boundary integration points NIb, the CPU time and the RMS error defined by rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N h i2ffi 1 X analytical numerical RMSError ¼ Ci ð22Þ  Ci N i¼1 where N is the total number of grid nodes. Table 1 gives also the relative mass error in the domain given by



Massnumerical  Massanalytical



ð23Þ MassError ¼

Massanalytical These experiments show that • For both methods, for a given number of integration points, the uniform mesh provides the most accurate results with less CPU time. • The accuracy of the results is much more dependent on the number of integration points used to evaluate

the boundary integral with the standard scheme than with the new scheme (see variation of Cmax versus NIb for the mesh (c) with the standard and the new schemes). • To obtain accurate results with the standard scheme, the number of integration points used to evaluate integrals is depending on the mesh. For the uniform mesh (a), Ns = 3 is sufficient to obtain a good solution (Cmax = 1.01). For the mesh (b), more integration points are required. A good solution (Cmax = 1.01) is obtained with Ns = 7. The last unstructured mesh (c) needs however much more integration points (Ns = 15) to obtain a good solution (Cmax = 1.05). Increasing the number Ns and so the number of particles at the boundary edge allows to achieve sufficient mass for each element at the new time level which reduces non-physical oscillations in the solution. The difficulty with the standard scheme is that the number Ns cannot be known a priori. Therefore, if the results present non-physical oscillations, the number Ns is increased and a new simulation is run until satisfactory results are obtained. This procedure is not needed with the new scheme of the ELLAM. • The use of strategic integration points with the new scheme allows a strong reduction of the number of

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

1063

Table 1 Effect of the mesh for the new and the standard schemes Cmax

MassError

RMSError

NIb

CPU

Mesh (a) Standard (Ns = 3) Standard (Ns = 7) Standard (Ns = 15) New scheme (P = 1) New scheme (P = 4)

1.01 1.01 1.01 1 1

<103 <103 <103 <103 <103

1.02 · 103 1.03 · 103 1.05 · 103 4.7 · 104 4.7 · 104

5100 11 900 25 500 320 1280

0.5 0.52 0.58 1.5 · 102 5 · 102

Mesh (b) Standard (Ns = 3) Standard (Ns = 7) Standard (Ns = 15) New scheme (P = 1) New scheme (P = 4)

1.09 1.01 1.005 1 1

<103 <103 <103 4 · 103 <103

1.25 · 102 1.1 · 102 1.1 · 102 1.3 · 102 1.1 · 102

8571 19 999 42 855 331 1319

0.625 0.7 0.9 3 · 102 8 · 102

Mesh (c) Standard (Ns = 3) Standard (Ns = 7) Standard (Ns = 15) New scheme (P = 1) New scheme (P = 4)

2.44 1.3 1.05 1.001 1

<103 <103 <103 8 · 103 <103

8.5 · 102 1.8 · 102 1.1 · 102 1.2 · 102 8 · 103

19851 46 319 99 255 1669 6707

1.9 2.1 2.6 6 · 102 0.2

integration points required to obtain the good solution. Indeed, for the unstructured mesh (c), the new scheme needs only NIb = 1669 integration points to obtain a solution with Cmax = 1.001. The standard scheme requires NIb = 99 255 to obtain a solution with Cmax = 1.05. • Even with P = 1, the new scheme gives accurate solution for all meshes. • The total mass in the domain with both methods is close to the analytical one. The new scheme with P = 4 leads to better mass balance than with P = 1. • For all meshes, the new scheme is much more efficient than the standard one since it gives a more accurate solution with less CPU time than the standard scheme. The efficiency of the new scheme is more pronounced for unstructured meshes.

3.3. Numerical integration of spatial integral at the old time level (12) The spatial integral (12) is also evaluated using strategic integration points defined at the new time level. For each element, P particles are set on the strategic point locations xj at tj = Dt with their corresponding 2D spatial weights ðW space Þ. Each particle is then backj tracked during Dt, from its position xj to the departure position xj (at t ¼ tj ). If tj ¼ 0, then xj is located inside the domain and the value of the concentration at the departure point is determined using interpolation based on concentrations at neighbouring mesh points, where they are known. This procedure is known to add significant numerical diffusion when used with several time steps.

4. Numerical diffusion with the ELLAM In the standard forward tracking approach, the integration points and weights are defined at the old time level. If strategic integration points are included, weights are reassigned (via the compound trapezoidal rule) once the strategic points are backtracked [4,5]. In the developed new scheme, strategic integration points and weights are assigned at the new time level and then backtracked to the old time level without redistributing the weights. This procedure is close to a backtracking numerical integration approach. On the other hand, it is known that the modified method of characteristics (MMOC) (e.g., [3]), based on the backtracking approach, introduces numerical dispersion, especially for sharp front problems, when a lower-order interpolation scheme is used. Higher-order interpolation schemes reduce numerical dispersion but are less computationally efficient and can lead to artificial oscillations for sharp front problems. Because interpolations are done at each time step, numerical errors are cumulated over time when many time steps are used. The algorithm developed here avoids this problem by continuously tracking the integration points. 4.1. Simulations with several time steps For simplicity, let us consider that T = m · Dt, that the inflow boundary flux is only due to advection and that the domain is initially free of contaminant. The idea of the algorithm is to evaluate the distribution of the concentration over the domain at the end of each time step with only the strategic integration points

1064

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

corresponding to that time step. Interpolation problems are minimized since we use continuous characteristics and only concentration changes due to dispersion are interpolated to obtain the concentration at the foot of each characteristic. The procedure is as following: • At each element, P particles are located at strategic locations xj at tj = T with the corresponding 2D area weights ðW space Þ. j • All particles are backtracked with T from their initial position xj (at tj = T) to their arrival position xj (at t ¼ t0j ). Only particles with t0j > 0 (i.e. xj is located at the inflow boundary), are kept. • Particles are then sorted according to their arrival time. For each time step k, nk is the number of boundary particles P kj such that t0j

is obtained by interpolating C disp from the C disp j i three nodes of the arrival element containing x1j . In this way, only the dispersive part of the transport is interpolated at each time step, which reduces the numerical diffusion. • At the second time step, previous particles j > n1 start at t = Dt at x1j and are forward tracked to x2j at t = 2Dt into the domain. – The boundary integral is updated (if the prescribed boundary concentration changes in time) with the first n1 particles using (26). – The mass at the old time level is obtained only with n2 particles: Z n2   X C 1 x1 ðxÞdx  C x1j ; Dt xi x1j ; Dt W space j X

j

ð24Þ

ð28Þ

We define also the local departure time tj of each boundary particle P kj by

where Cðx1j ; DtÞ for the n2 particles is obtained from (27). – Again, the concentration at the foot of each characteristic is updated to reflect the change due to dispersion

ðm  kÞDt 6

< ðm  k þ 1ÞDt

tj ¼ t0j  ðm  kÞDt

ð25Þ xj

• At the first time step, all particles start at at t ¼ tj 1 and are forward tracked to xj into the domain at t = Dt. – The boundary integral is evaluated using only the first n1 boundary particles P 1j with Z T Z ½ðVCÞxi noX dl dt 0

oX

n1   X C xj ; tj xi xj ; tj W space  j

ð26Þ

j

–

–

–

–

where xi ðxj ; tj Þ is evaluated using xi ðxj ; tj Þ ¼ xi ðx1j ; DtÞ. This allows us to solve Eq. (6) to obtain the concentration C iadvþdisp at each node of the mesh. Theses values of the concentration are due to both advection and dispersion phenomena. At the end of the time step, the value C adv i , which is the concentration due to the advection only, is computed at each node i. This is done using (6) with zero dispersion. This procedure is not CPU time consuming. It is performed at the element level since mass lumping is used. For each node, C disp , which is the concentration i due to the dispersion only, is computed using C disp ¼ C advþdisp  C adv i i i . For all boundary particles, except the first n1 particles, the concentrations are updated to reflect the change due to dispersion   C x1j ; Dt ¼ C xj ; tj þ C disp for j > n1 ð27Þ j ; where Cðxj ; tj Þ is the concentration at the foot of the characteristic (for the first time step, it corresponds to the prescribed boundary value) and

Cðx2j ; 2DtÞ ¼ Cðx1j ; DtÞ þ C disp j ;

for j > n2

ð29Þ

is obtained by interpolating C disp after the secC disp j i ond time step computation. • The other time steps are treated as the second time step. It is important to notice that contrarily to Eulerian methods, when concentration at Dirichlet boundaries is time varying, the ELLAM allows to take into account changes in the prescribed concentration when using a single time step simulation. In this case, boundary particles with different departure times will have different initial concentrations. 4.2. Simulations with non-zero initial conditions In the previous procedure, the distribution of the concentration at t = 0 is not taken into account. Moreover, at the end of each time step, the concentration in some elements may change only by dispersion and no particle reaches these elements. However, for the following time steps, the mass in these elements has to be transported by advection even if no particle will start from these elements. To solve this problem, the integral at the old time level is treated as following: • For each time step k (k = 1, . . . , m), P particles per element located at the strategic point locations xj at tj = T (with the corresponding 2D area weights W space ) are used. These particles are backtracked with j

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

(m  k)Dt from their position xj (at tj = T) to their departure position xj (at t ¼ t0j ). Only N ks particles where xj is located inside the domain ðt0j ¼ ðm kÞDtÞ are kept. Let us Ns the total number of Pcall m these particles ðN s ¼ k¼1 N ks Þ. The heads xj of the Ns characteristics correspond to the strategic point locations. The feet xj of the characteristics are located somewhere in the domain. • The concentration at the foot of each characteristic is obtained by linear interpolation of the initial concentration at the nodes of the mesh. • For the first time step, all Ns particles are forward tracked with Dt from xj to x1j . – Only the first N 1s particles are used to evaluate the spatial integral at the old time level with Z N 1s   X 0 0 C x ðxÞdx  C x0j ; 0 xi x0j ; 0 W space j X

j

ð30Þ – The concentration is modified for the other particles, as previously   C x1j ; Dt ¼ C xj ; 0 þ C disp for j > N 1s ð31Þ j ; where Cðxj ; 0Þ is the concentration at the foot of the characteristic (for the first time step, it corresponds to the initial concentration) and C disp j is obtained as previously by interpolating C disp of i the three nodes of the element containing x1j . • At the second time step, previous particles j > N 1s start at t = Dt at x1j and are forward tracked to x2j at t = 2Dt into the domain. – The mass at the old time level corresponds to the mass which arrived from the boundary (n2 particles) at the first time step given by Eq. (28) and the initial mass in the domain (N 2s particles) Z

C 1 x1 ðxÞdx 

X

n2   X C x1j ; Dt xi x1j ; Dt W space j j

þ

N 2s   X C x1j ; Dt xi x1j ; Dt W space j j

ð32Þ Cðx1j ; DtÞ

is obtained from (31). where – Again, the concentration at the foot of each characteristic is modified using   C x2j ; 2Dt ¼ C x1j ; Dt þ C disp for j > N 2s j ; ð33Þ C disp is obtained by interpolating C disp after the j i second time step computation. • The further time steps are treated as the second time step.

1065

4.2.1. Remark For many transport problems, only a small part of the domain is reached by the solute. Therefore, if several time steps are necessary, we can derive benefit of the fact that ELLAM allows large time steps as following: • The transport problem with a single time step equal to the final simulation time T is computed first. • All elements with a relative concentration greater than a given tolerance value (for example e = 105) will constitute the transport domain. • The same procedure, as previously, is done. However particles are backtracked only from elements of the transport domain. With this procedure, we have more computational effort since we add a simulation with a large time step T in order to define the transport domain. However, the computational domain can be strongly reduced. For practical problems, this procedure can be very attractive, since the transport domain is usually much smaller than the whole domain. Different formulations of the new scheme and socalled standard formulation are compared in the next section.

5. Efficiency of the new implementations with a uniform flow field 5.1. Computation with many time steps To check the efficiency of the new ELLAM when working with several time steps, the previous problem is solved with a final simulation time of T = 40 [T] divided into 20 time steps of 2 [T]. The discretization is the unstructured mesh (c). Three methods are used • MS_P1_o1: the ELLAM with standard linear interpolation at the end of each time step. No continuous tracking is used and the concentration at the feet of characteristics are obtained by linear interpolation of concentration at nodes. Strategic integration points with corresponding weights are defined at the new time level. We use P = 1 particle per element for the whole domain. • MS_P1_o4: similar to MS_P1_o1 with P = 4 strategic integration points per element. • MS_P1_a: the new implementation of the ELLAM with continuous tracking of characteristics and with P = 1 particle per element for the whole domain. • MS_P1_b: in this case, a first run with a single time step Dt = 40 [T] is done to estimate the transport domain. Then, only 1 particle per element is set for triangles located inside the transport domain.

1066

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

Results are shown in Fig. 5, the following conclusions can be made: • Significant numerical diffusion is obtained with the two first methods where the continuous tracking procedure is not used.

• A small improvement of the solution accompanied with a marked increase in the CPU time is obtained when the number of strategic integration points per element is increased. • No significant differences between the two last methods can be detected.

RMS-Error

CPU

0.11

1.8

0.1

6.7

1.2×10 -2

3.5

1.2×10-2

2.4

-0.1

-0.01

0.1

0.5

0.99

1.01

Concentration over the domain

Analytic

40 30 20 10 0

0

20

40

60

80

100

MS_P1_o1

40 30 20 10 0

0

20

40

60

80

100

MS_P1_o4

40 30 20 10 0

0

20

40

60

80

100

MS_P1_a

40 30 20 10 0

0

20

40

60

80

100

MS_P1_b

40 30 20 10 0

0

20

40

60

80

100

Fig. 5. Results after a simulation time of 40 days using 20 equal time steps.

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

• RMS errors remain small for MS_P1_a and MS_P1_b and are very close to the value obtained when using a single time step (Table 1). Numerical

diffusion due to the use of many time steps is not significant and a good mass balance is observed for simulations with one or several time steps.

Concentration at y = 20 m

a

1067

RMSError

CPU

1.0

Concentration

Analytic

0.8

0.6

0.4

0.2

0.0 0

20

40 X (m)

60

80

100

1.0

Concentration

DFEM

0.8

0.6

-2

34.7

1.1×10 -2

4.26

1.1×10 -2

0.23

3.7×10

0.4

0.2

0.0 0

20

40

60

80

100

X ( m)

0.8 Concentration

Standard ELLAM

1.0

0.6

0.4

0.2

0.0

0

20

40 X (m)

60

80

100

1.0

Concentration

New ELLAM

0.8

0.6

0.4

0.2

0.0 0

20

40 X (m)

60

80

100

Fig. 6. (a) Efficiency of the new ELLAM for case 1 (aL = 0.002, aT = 0.0005, Dm = 0). (b) Efficiency of the new ELLAM for case 2 (aL = 0.2, aT = 0.05, Dm = 0). (c) Efficiency of the new ELLAM for case 3 (aL = 2, aT = 0.5, Dm = 0).

1068

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

• MS_P1_a requires about 50% more CPU time than MS_P1_b. This ratio is expected to be much greater for practical problems where usually the transport domain is very small compared to the whole domain. Therefore, an important reduction of the computational effort should be obtained.

5.2. Efficiency of the new ELLAM The new implementation of the ELLAM is compared to two numerical codes. The first one is based on explicit Galerkin Discontinuous Finite Element Method (DFEM). An explicit scheme is used for the

Concentration at y = 20 m

b

RMSError

CPU

1.1×10-2

36.24

1.1×10-2

4.31

1.1×10-2

0.26

1.0

Concentration

Analytic

0.8

0.6

0.4

0.2

0.0

0

20

40 X (m)

60

80

100

1.0

Concentration

DFEM

0.8

0.6

0.4

0.2

0.0 0

20

40 X (m)

60

80

100

0.8 Concentration

Standard ELLAM

1.0

0.6

0.4

0.2

0.0 0

20

40 X (m)

60

80

100

1.0

Concentration

New ELLAM

0.8

0.6

0.4

0.2

0.0 0

20

40 X (m)

60

80

100

Fig. 6 (continued)

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

convective part of the solute transport equation, an implicit scheme for its dispersive part [14]. The dispersive part of the transport equation is solved with mixed hybrid finite element method. Due to the explicit scheme, the Courant criterion has to be fulfilled which

can lead to very small time steps when used for unstructured meshes. The second code is based on the standard ELLAM using Gauss points to perform the numerical integration in space and Simpson rules in time.

Concentration at y = 20 m

c

1069

RMSError

CPU

6×10 -3

43.5

2.2×10-2

4.37

2.1×10 -2

0.27

1.0

Concentration

Analytic

0.8

0.6

0.4

0.2

0.0 0

20

40 X (m)

60

80

100

1.0

Concentration

DFEM

0.8

0.6

0.4

0.2

0.0 0

20

40 X (m)

60

80

100

0.8 Concentration

Standard ELLAM

1.0

0.6

0.4

0.2

0.0 0

20

40 X (m)

60

80

100

1.0

Concentration

New ELLAM

0.8

0.6

0.4

0.2

0.0 0

20

40 X (m)

60

80

100

Fig. 6 (continued)

1070

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

The results of the new ELLAM, Standard ELLAM and DFEM are given for different ranges of grid Peclet number (Pe) using the following dispersion parameters: • Case 1: aL = 0.002 [L], aT = 0.0005 [L], Dm = 0 [L2/ T], i.e. a grid Peclet number varying from 101 to 1467. • Case 2: aL = 0.2 [L], aT = 0.05 [L], Dm = 0 [L2/T], i.e. a grid Peclet number varying from 1.01 to 14.67. • Case 3: aL = 2 [L], aT = 0.5 [L], Dm = 0 [L2/T], i.e. a grid Peclet number varying from 0.101 to 1.467.

ðx  xC Þ2 þ ðy  y C Þ2 C 0 ðx; yÞ ¼ exp  2r2

! ð34Þ

where xC, yC, and r are the centered and standard deviations of the Gaussian pulse. The corresponding analytical solution for (6) with a constant diffusion coefficient D is given by ! 2 2 2r2 ðx  xC Þ þ ðy  y C Þ exp  Cðx; y; tÞ ¼ 2 2r þ 4Dt 2r2 þ 4Dt ð35Þ

Simulations are done with the unstructured mesh (c) and the results are given in Fig. 6. The simulation is run with a single time step Dt = T = 40 [T] and the Courant number is varying inside the domain from 27.9 to 480.58. The simulations show that • For case 1, the DFEM provides more numerical diffusion and leads to more errors than both ELLAM schemes. The standard ELLAM gives very small oscillations, no oscillation was detected for the new implementation. • For case 2, the three methods provide accurate solutions with small and similar errors. • For case 3 with dispersive dominant transport, the DFEM leads to better results than both ELLAM. • Results show that for advective dominant transport, the new scheme is much more efficient than the standard scheme and than the DFEM. When dispersion is very dominant, the error of ELLAM increases if a single time step is used, because of the one-step backward Euler approximation in time used for the dispersion term. Indeed, when using MS_P1_a with four time steps, the error is similar to the single step ELLAM error for the two first cases. However, for dispersive dominant transport (case 3), the error of MS_P1_a with four time steps is 50% less than the error obtained with the one time step computation. Therefore, a minimum number of time steps has to be used when working with very diffusive transport problems.

6. Efficiency of the new implementations with a non-uniform flow field We consider the transport of a two-dimensional rotating Gaussian pulse which is a very standard test case [15]. The spatial domain is X = (0.5, 0.5) · (0.5, 0.5) [L2], the rotation field is V1(x, y) = 4y [L/ T], and V2(x, y) = 4x [L/T], the final time simulation is T = p/2 [T] which corresponds to the time period required for one complete rotation. The initial concentration is given by

where x ¼ x cosð4tÞ þ y sinð4tÞ and y ¼ x sinð4tÞ þ y cosð4tÞ. This problem provides ADE with variable velocity and known analytical solution. The problem is advection dominant in most of the domain and becomes diffusion dominant near the origin. 6.1. Structured meshes In the numerical experiments, we choose the same data as Wang et al. [15] with D = 104 [L2/T], xC = 0.25 [L], yC = 0 [L], and r = 0.0447 [L2/T]. A first triangulation (Mesh (d)) is obtained by subdivision of square elements of Dx ¼ Dy ¼ 321 (by joining the center of each element to its four corners). A second triangulation (Mesh (e)) is obtained in the same way starting with the finer spatial grid used by Wang et al. [15] with Dx ¼ Dy ¼ 641 . The initial concentration defined by (34) has a minimum value 0 and a maximum value 1. Due to diffusion, the analytical solution (35) after one complete rotation has a minimum value equal to 0 and a maximum value equal to 0.865 (Fig. 7). Table 2 gives the maximum value of the concentration, the RMS error, the maximum error, the CPU time and the relative mass error for both meshes. Results are given for the DFEM, the standard ELLAM and the new ELLAM when using 20 time steps and when using a single time step. In the standard ELLAM, we use 4 integration points (NA = 4) per triangle to approximate the integral (12). From this table, we can see that • The new ELLAM gives similar results when using one or 20 time steps. The method is very efficient, a good solution can be obtained for the mesh (d) in 0.39 s and in 3.2 s for the mesh (e). • The standard ELLAM requires more CPU time than the new ELLAM when using a single time step and less CPU time when using 20 time steps. With the new ELLAM, the CPU time increases faster, when the number of time steps increases, than with the standard ELLAM. With the so-called standard approach, the same particles are used to compute

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

1071

Fig. 7. Initial concentration, analytical solution and ELLAM solution for the rotating Gaussian pulse problem.

Table 2 Results for the two-dimensional rotating Gaussian pulse with structured meshes

Mesh (d) Standard_ELLAM New ELLAM

nDt

Max.

RMSError

Max. error

1 20

0.817 0.836

4 · 104 5.7 · 104

4.7 · 102 7.6 · 102

0.9 1.3

<103 3 · 103

1 20

0.844 0.853

4.6 · 104 4 · 104

3.6 · 102 3.3 · 102

0.39 2.5

2 · 103 3 · 103

0.6

1.8 · 103

0.22

11.1

<103

1 20

0.865 0.87

1.3 · 104 3.5 · 104

2 · 102 5.4 · 102

7.2 10

<103 2 · 103

1 20

0.868 0.86

104 1.2 · 104

6 · 103 1.4 · 102

3.2 17

<103 <103

0.78

3 · 104

7.8 · 102

DFEM Mesh (e) Standard_ELLAM New ELLAM DFEM

the concentration at all time levels. With the new implementation, each time step k uses its specific spatial integration points N ks . For uniform meshes with several time steps, the new ELLAM requires more integration points and become less efficient than the standard one. • For the coarse mesh (d), the new ELLAM gives better results (maximum of concentration equal to 0.843) than the standard ELLAM (maximum of concentration equal to 0.817).

CPU

105

MassError

<103

• For the fine mesh (e), results of both methods are very close to the analytical solution and the mass error remains small. • The mass error decreases for the finer mesh (e). This error is in general more important when we use 20 time steps than with a single time step. • For both meshes, DFEM gives less satisfactory results (maximum of concentration equal to 0.6 for mesh (d) and 0.78 for mesh (e)) than ELLAM. The required CPU time is significantly higher than for

1072

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

xC = 0.25 [L], yC = 0 [L], and r = 0.0447 [L2/T] is solved. The final time simulation is T = p/6 [T], which corresponds to the time required for a rotation of 2p/3. The final concentration has a minimum value equal to 0 and a maximum value equal to 1 since dispersion is zero. Standard ELLAM is used with different numbers of integration points per element (NA = 4, . . . , 192 in integral (12)). The comparison between both ELLAM formulations is given in Table 3. Following conclusions can be drawn:

both ELLAM formulations. Results show that DFEM is also more sensitive to the mesh size than ELLAM.

6.2. Unstructured mesh The same problem is now simulated with the unstructured mesh (f) shown in Fig. 8. This mesh contains a very fine discretization in the second quarter of the domain. A purely convective problem with D = 0 [L2/T],

1

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fig. 8. The unstructured mesh for the rotating pulse problem.

Table 3 Results for the two-dimensional rotating Gaussian pulse with the unstructured mesh

Standard_ELLAM NA = 4 NA = 16 NA = 24 NA = 48 NA = 96 NA = 192

nDt

Max.

RMSError

Max. error

1 1 1 1 1 1

12.62 6.38 3.69 2.08 1.4 1.12

0.12 8.82 · 102 2.69 · 102 1.19 · 102 5.26 · 103 2.7 · 103

11.65 8.08 2.7 1.1 0.44 0.14

0.883

4.98 · 103

0.12

0.999 0.999

1.4 · 103 1.4 · 103

3 · 102 4.5 · 102

CPU 1.98 5.18 7.25 13.58 25.98 52.31

MassError 1.4 · 103 1.4 · 103 1.4 · 103 1.4 · 103 1.4 · 103 1.4 · 103

DFEM New_ELLAM P=1

1 10

127 1.01 3.23

<103 2 · 103 2.1 · 103

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

• The standard ELLAM gives results with oscillations, which decrease as the number of spatial integration points increases. Indeed, with this scheme, the integration points and weights are defined at the old time level. For unstructured meshes, quadrature weights will not necessarily sum up to the correct volume for each mesh element, on the new time level which leads to oscillations. Increasing the number of spatial integration points allows to achieve sufficient mass for each element at the new time level which reduces non-physical oscillations in the solution. • The new ELLAM (maximum: 0.999) is more accurate than the DFEM (maximum: 0.883) and than the standard ELLAM (non-physical oscillations). Indeed, with this scheme, the strategic integration points and weights are defined at the new time level which avoids non-physical oscillations. • The new ELLAM gives accurate solution when used with several time steps on unstructured meshes since continuous tracking with specific strategic integration points for each time step is used. This make it less efficient than with a single time step. • The new ELLAM is highly accurate for unstructured meshes. Indeed, the good solution is obtained in only 1.01 seconds. The DFEM require 127 s and the standard ELLAM requires 57.31 s and both give less satisfactory results.

7. Conclusion A new implementation of the Eulerian–Lagrangian Localized Adjoint Method (ELLAM) for solving the Advection–Dispersion transport Equation on structured or unstructured triangular meshes has been developed and compared to another more standard ELLAM formulation and to a combination of discontinuous/mixed finite element methods. The standard ELLAM scheme is often used with a forward tracking approach where integration points and weights are defined at the old time level. For unstructured meshes, quadrature weights will not necessarily sum up to the correct volume for each mesh element, on the new time level which leads to oscillations. Increasing the number of integration points allows to achieve sufficient mass for each element at the new time level which reduces non-physical oscillations in the solution. Therefore, ELLAM becomes less efficient since it requires a lot of integration points for highly unstructured meshes. Moreover, the required number of integration points cannot be known a priori and this number is problem-dependent which can make ELLAM less attractive for practical cases. Contrarily to standard ELLAM formulations, the new scheme requires a very limited number of integra-

1073

tion points (usually one per element) even for unstructured meshes. With this scheme, only strategic points are used as numerical integration points. Location of strategic integration points and weights are assigned at the new time level and then backtracked to the old time level without redistributing the weights. Numerical diffusion due to interpolation when using the ELLAM with many time steps is also addressed. A new algorithm is developed in order to avoid excessive numerical diffusion. The idea of this algorithm is to evaluate the distribution of concentration over the domain at the end of each time step with only the strategic integration points corresponding to that time step. Interpolation problems are minimized since continuous characteristics are used and only changes due to dispersion are interpolated to obtain the concentration at the foot of each characteristic. When several time steps are necessary, an algorithm can be implemented in order to derive benefit of the fact that ELLAM allows large time steps to avoid excessive calculations in parts of the grid that will not contain significant concentrations. Numerical experiments for grid Peclet numbers ranking from 1 to 1000 show the efficiency of the developed scheme. The new scheme gives better solution and requires less CPU time than the standard scheme. The high performance of the developed scheme is more pronounced for unstructured meshes and large time steps. The new ELLAM is also shown to be much more efficient compared to a high order Eulerian scheme based on discontinuous and mixed finite element methods. The suggested implementations make the ELLAM very attractive, especially because large time steps can be used without fulfilling the Courant criterion on unstructured meshes. This makes the application of the new ELLAM to field cases very attractive and competitive. However, the computational implementation of this method is more complicated, due to its specific particles handling.

Acknowledgement The authors wish to thank the anonymous reviewers for their helpful remarks, which greatly improved the quality of this paper.

References [1] Binning P, Celia MA. A forward particle tracking Eulerian– Lagrangian localized adjoint method for solution of the contaminant transport equation in three dimensions. Adv Water Resour 2002;25:147–57. [2] Celia MA, Russell TF, Herrera I, Ewing RE. An Eulerian– Lagrangian localized adjoint method for the advection–diffusion equation. Adv Water Resour 1990;13:187–206.

1074

A. Younes et al. / Advances in Water Resources 29 (2006) 1056–1074

[3] Douglas Jr J, Russell TF. Numerical methods for convectiondominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures. SIAM J Numer Anal 1982;19:871–85. [4] Healy RW, Russell TF. A finite-volume Eulerian–Lagrangian localized adjoint method for solution of the advection–dispersion equation. Water Resour Res 1993;29:2399–413. [5] Healy RW, Russell TF. Solution of the advection–dispersion equation in two dimensions by a finite-volume Eulerian–Lagrangian localized adjoint method. Adv Water Resour 1998;21:11–26. [6] Leij FJ, Dane JH. Analytical solution of the one-dimensional advection equation and two or three-dimensional dispersion equation. Water Resour Res 1990;26:1475–82. [7] Morton KW, Priestley A, Su¨li E. Stability of the Lagrange– Galerkin method with non-exact integration. Math Model Numer Anal 1988;22:625–53. [8] Pollock DW. Semianalytical computation of path lines for finitedifference models. Groundwater 1988;26:743–50. [9] Priestley A. Exact projections and the Lagrange–Galerkin method: a realistic alternative to quadrature. J Comput Phys 1994;112:316–33. [10] Russell TF, Binning P. Oh No, not the Wiggles Again! In: Miller CT et al., editors. A revisit of an old problem and a new approach, proc computational methods in water resources. Amsterdam: Elsevier; 2004. p. 483–94. [11] Russell TF. Numerical dispersion in Eulerian–Lagrangian methods. In: Hassanizadeh SM et al., editors. Proc computational

[12]

[13]

[14]

[15]

[16]

[17]

[18]

methods in water resources. Amsterdam: Elsevier; 2002. p. 963–70. Russell TF, Heberton CI, Konikow LF, Hornberger GZ. A finitevolume ELLAM for three-dimensional solute-transport modeling. Groundwater 2003;41:258–72. Russell TF, Celia MA. An overview of research on Eulerian– Lagrangian localized adjoint methods (ELLAM). Adv Water Resour 2002;25:1215–31. Siegel P, Mose´ R, Ackerer P, Jaffre´ J. Solution of the advection dispersion equation using a combination of discontinuous and mixed finite elements. Int J Numer Methods Fluids 1997;24:595–613. Wang H, Dahle HK, Ewing RE, Espedal MS, Sharpley RC, Man S. An ELLAM scheme for advection diffusion equations in two dimensions. SIAM J Sci Comput 1999;20:2160–94. Younes A. An accurate moving grid Eulerian Lagrangian localized adjoint method for solving the one-dimensional variable-coefficient ADE. Int J Numer Methods Fluids 2004;45:157–78. Younes A, Ackerer P. Solving the advection–diffusion equation with the Eulerian Lagrangian localized adjoint method on unstructured meshes and non-uniform time stepping. J Comput Phys 2005;208:384–402. Younes A, Ackerer P, Mose´ R, Chavent G. A new formulation of the mixed finite element method for solving elliptic and parabolic PDE with triangular elements. J Comput Phys 1999;149:148–67.