Stability and Convergence of a Finite Element Method for Reactive Transport in Ground Water.
An explicit finite element method is used to solve the linear  convection--diffusion-reaction equations governing contaminant transport in ground water flowing through an adsorbing porous medium. The use of discontinuous finite elements for the convective part of the equations combined with mixed finite elements for the diffusive part renders the method for the concentration solution, which displays strong gradients, trivially conservative and fully parallelizable. We carry out a stability and convergence analysis. In particular, the method is proven to satisfy a maximum principle, to be total variation bounded, and to converge to the unique weak solution of the equations. Special attention is paid to the convective part of the equations. Numerical simulations are presented and discussed.
Introduction
. In this paper we propose and analyze a finite element method
for solving the linear convection-diffusion-reaction equation:
@t
which describes the transport of a solute in a fluid phase flowing through a porous
medium [1], [16]. In this case, is the concentration of the solute in
the fluid phase for which we solve (1.1), is the Darcy velocity of
fluid, \Phi is the volume fraction-dependent constant, D is the diffusion constant, and
is the first-order chemical reaction rate. This equation, while
formally parabolic, is more nearly hyperbolic in practice [4]. In recent years many
finite element methods have been proposed to solve this important partial differential
equation. The classes of optimal spatial methods and characteristic methods have
been extensively studied [2], [9], [15], [17], [18], for example. However, all these finite
element methods are defined by taking advantage of the parabolicity of the equation
for the concentration u. As a result, the solution of the differential equation is required
very smooth in the derivation of error estimates, and the constants for the
error estimates blow up as the coefficient of the diffusion term goes to zero.
In this paper we propose and analyze a finite element method for numerically
solving (1.1). It is similar to a finite element method introduced in [5], [3], [6], [10],
[11] in that we approximate the convective part of the equation using a upwinding
discontinuous finite element method or a upwinding finite volume method [20], [19].
We use, however, a mixed finite element method for the diffusive part of (1.1) [8].
The main advantages of this method are that it is trivially conservative and fully
parallelizable, and that it can capture discontinuities within a couple of elements
without producing spurious oscillations.
* Department of Mathematics and the Institute for Scientific Computation, Texas A&M Univer-
sity, College Station, Current address of Z. Chen: Department of Mathematics, Box 156,
Dedman College, Southern Methodist University, Dallas, Texas 75275-0156. Partly supported by the
Department of Energy under contract DE-ACOS-840R21400. email: zchen@golem.math.smu.edu,
ewing@ewing.tamu.edu.
A stability and convergence analysis is carried out here for the finite element
method for equation (1.1) in two space dimensions. While a stability analysis was
completed for the similar approach for the two-dimensional semiconductor device equations
in [6], we are here able to prove much stronger results than those obtained in [6].
Namely, besides a strong maximum principle, the boundedness of the total variation
and the modulus of continuity in time of the approximate solution is proven here; only
an estimate on the weak derivatives of the approximate solution is given in [6]. These
properties suffice to show that the numerical method converges to the weak solution
of the differential equation; in [6], however, convergence of the approximate solutions
to the weak solution is proven under the assumption that there is a convergent subse-
quence. It is also emphasized that this paper contains the first stability analysis for the
two-dimensional equation (1.1) with the diffusion term included and the first convergence
analysis for (1.1) with the boundary conditions. The properties derived in this
paper will be exploited in a forthcoming paper where error estimates will be obtained
with minimum requirements on the solution and with the property that the constant
for the error estimates does not involve the small diffusion coefficient. Especially, the
error estimates apply to the case of D equal to zero.
The equation (1.1) is completed by specifying the boundary and initial conditions:
(1.2a)
(1.2c)
the endpoints of its segments, and - denotes the normal unit-vector to @
\Omega\Gamma The
boundary conditions need to be modified properly in the case of Namely,
only the inflow boundary condition is imposed for the concentration (see (3.9) below).
Moreover, in this case, note that, while the equation (1.1) is analogous to a classical
conservation law, the value of the Darcy velocity V at a point (t; x; y) contains the
information of all the values of the solution u(t; \Delta; \Delta) on \Omega\Gamma Hence a perturbation of
the solution u at any given point of the domain has a global effect immediately. This
is in sharp contrast with the classical conservation laws where local perturbations of
the solution have a local effect in finite time.
The rest of the paper is organized as follows. The finite element method is defined
in the next section. Then, in x3 we state and discuss our main results on a maximum
principle (Theorem 3.1), a total variation boundedness of the scheme (Theorem 3.2),
continuity with respect to data (Theorem 3.3), and convergence to the weak solution
(Theorem 3.4). The proofs of these properties are carried out in x4, x5, x6, and
x7, respectively. Numerical results are displayed in x8. These numerical results are
devised to test the performance of the method and to indicate the order of convergence.
Finally, a concluding remark is given in x9.
2. The finite element method. In this section we define the finite element
method for approximating the solution of the differential system (1.1). Toward that
end, let fx i+1=2 g nx
i=0 \Theta fy j+1=2 g ny
j=0 be a partition
of\Omega with x
n=0 be a partition of [0; T
we introduce the following notation
1-i-nx
1-i-ny
\Deltay
\Deltat
We tacitly assume that each exterior edge has imposed on it either Dirichlet or Neumann
conditions, but not both. Associated with these partitions, we introduce the
spaces
i;j x; a 3
i;j y); a k
i \Theta I y
represents the constant value w(x; y); (x; y) 2 I x
i \Theta I y
the constant w(t); t 2 J n , if w 2 W \Deltat . For notational and expositional convenience,
let \Deltax \Deltay ny , \Deltax
\Deltay j+1 )=2,
Finally, define the notation
Let PQh , PWh , and PW \Deltat denote the L 2 -projections into Q h , W h , and W \Deltat ,
respectively. To discretize (1.1), we first discretize the data as follows:
(2.1c)
The subscript h is omitted below when no ambiguity occurs. Then the approximate
solution
\Deltat\Omega W h is required to satisfy the equation, for
\Deltat n
\Deltay j
\Deltay j
i;j
(2.2a)
where
(2.2c)
and the function q
\Deltat\Omega Q h is the solution of
@\Omega D
After the mass matrix has been mass-lumped [22], the expression for the degrees of
freedom of q h is taken as follows:
(2.2d)
Finally, the Neumann boundary condition (1.2a) is discretized by the usual reflection
principle, and on
defined by uD;\Deltat . This implies that, if lies on
the Neumann boundary
0;j in (2.2) and the subsequent analysis is calculated
by
if it is on the Dirichlet boundary
0;j is computed by
D;\Deltat
Similar extensions hold for u n
nx+1;j , and u n
i;ny +1 in (2.2) and the subsequent
analysis.
Note that the lowest-order Raviart-Thomas mixed method [21] over rectangles
has been used in (2.2a). Since the elements in Q h have continuous normal components
on interelement edges, the numerical fluxes f n
1;i\Gamma1=2;j and f n
2;i;j \Gamma1=2 in (2.2b) and
(2.2c) are well defined. Furthermore, if appropriate approximations of the coefficient
are introduced and the mass-lumping technique is used as in (2.2d) and (2.2e),
the conservative scheme (2.2a) can be deduced from the discontinuous finite element
method [7], [12] or from the finite volume method [20], [19] combined with the mixed
finite element method [22]. Finally, the scheme applies to the case of
The following approximation properties are used later [14], [21]:
(2.3a)
for each n. Moreover, since the operator PQh is locally defined, we have for each n:
(2.3c)
is independent of i and j.
3. Stability and convergence results. In this section we state and discuss
the stability and convergence results of the scheme (2.2). Let Q We assume
that the data satisfy the following conditions:
(3.1a)
(3.1c)
(3.1d)
For expositional convenience, let
Theorem 3.1 (Stability). Suppose that (3.1a), (3.1b), (3.1f), and for
the following Courant-Friedrichs-Lewy (CFL) condition are satisfied:
\Deltat
\Deltax i+1=2
\Deltax i\Gamma1=2
\Deltay j
\Deltay j+1=2
\Deltay j\Gamma1=2
In addition, if
then we have
Obviously, since K - 0, (3.4) is satisfied if V 1 is nondecreasing in x and V 2 is
nondecreasing in y, or div V is uniformly positive by the definition of PQh [21].
nx
ny
and set
Theorem 3.2 (TVB). Assume that (3.1) and for
CFL condition are satisfied:
\Deltat
=\Deltay ?
Then there is a constant C 1 depending solely on the data and T such that
ae
+\Deltay ?
oe
We remark that either in the case of or in the case of uD being constant
in space, (3.7) shows that the total variation of the solution u h is bounded. The
latter case means that the total variation of the solution u h in the one-dimensional
case is always bounded since uD is constant in this case. The numerical experiments
given in x8 show that the bounds in (3.7) and (3.8) below are sharp when D 6= 0 and
in the sense that the left-hand sides of the inequalities (3.7)
and (3.8) blow up as \Deltax   or \Deltay   converges to zero (see Example 4 in x8).
In the following, v h stands for the approximate solution of (1.1) and (1.2) with
the data v init and v D satisfying the conditions (3.1a), (3.1c), and (3.1d).
Theorem 3.3 (Continuity with respect to data). Assume that the hypotheses
of Theorem 3.1 are satisfied for both sets of data. Then there exists a constant
depending only on the data and T such that
ae'
\Deltay ?
\Theta jju
oe
As for the convergence result, we now consider a simple case where In this
case Theorem 3.2 implies the total variation boundedness of the numerical scheme as
remarked above, which together with Theorem 3.3 yields the following convergence
result (see x7). For nonzero D, concrete error estimates for the numerical scheme (2.2)
will be obtained in the work mentioned earlier.
In the simple case where the boundary conditions (1.2a) and (1.2b) are
replaced by the following inflow boundary condition:
where
0g. We now extend the numerical flux
introduced in (2.2) to the general setting:
Also, we define
2\Omega g:
Then a weak solution of the differential equation given by (1.1) with
(1.2c) is defined to be a function u
BV(\Omega\Gamma2 satisfying the weak formulation:
ft=0g\Theta\Omega
where (\Delta; \Delta) S denotes the inner product in L 2 (S) for some set S. Note that the role
of the flux f is to select the correct boundary value for u, and that the smoothness
hypothesis on V guarantees the uniqueness of weak solution to (3.10).
Theorem 3.4 (Convergence). Assume that the hypotheses of Theorem 3.2
are satisfied. Then the sequence fu h gh?0 produced by the scheme (2.2) converges in
(\Omega\Gamma6 to the unique solution of (3.10). Moreover, u
BV(\Omega\Gamma3 .
4. Proof of the maximum principle. In this section we prove Theorem 3.1.
Let
Lemma 4.1. Suppose that
\Deltat n
\Deltat n
\Deltay j
D\Deltat n
\Deltax i+1=2
\Deltax i\Gamma1=2
D\Deltat n
\Deltay j
\Deltay j+1=2
\Deltay j \Gamma1=2
i;j \Deltat n - 0:
Then, if
we have, for
ae
I x
oe
In addition, if
we have
Proof. For it follows from (2.2) that
i;j
i\Gamma1;j
where
A n
\Deltat n
A n
\Deltat n
\Deltay j
\Deltay j \Deltay j+1=2
\Deltat n
\Deltat n
\Deltay j
D\Deltat n
\Deltax i+1=2
D\Deltat n
\Deltay j
\Deltay j+1=2
+\Deltay j \Gamma1=2
\Deltat n
\Deltat n
\Deltay j
\Deltay j \Deltay j \Gamma1=2
Then, by (4.1), we see that
A n
so that, by (4.2),
Furthermore, by the definition of Q h and (4.2),
\Deltat n
\Deltat n
\Deltay j
i;j \Deltat n
i;j \Deltat n
which implies (4.3) immediately since K - 0, and together with (4.4) yields (4.5).
Lemma 4.2. If for
ij is defined as in Theorem 3.1, then (4.1) is satisfied.
The lemma follows obviously from the inequality (4.1) and the definition of D ?
ij .
We are now ready to prove Theorem 3.1 by means of induction on n.
Proof of Theorem 3.1. For the results (3.3) and (3.5) follow trivially from
the assumption (3.1a). Let the results be true up to n. By Lemma 4.2 and (2.3a), (3.2)
clearly implies (4.1). Then iterating (4.3) on n and using the induction hypothesis
and (2.3b) yield that
Consequently, by (3.1a), (3.3) follows.
If (3.4) is true, so is (4.4). Then, in this case, it follows from (4.5) and the
induction hypothesis that
which implies (3.5) by (3.1a).
5. Proof of total variation boundedness. In this section we prove Theorem
3.2. In order to fix ideas, let
other cases can be treated similarly.
Lemma 5.1. For
\Deltat n
\Deltax 3=2 \Deltax 1
\Deltat n
\Deltay j
\Deltay j+1=2 \Deltay j
ae
\Deltat n
\Deltat n
\Deltay j
\Deltat n
\Deltay j
D\Deltat n
D\Deltat n
\Deltay j \Deltay j+1=2
D\Deltat n
\Deltay j \Deltay j \Gamma1=2
oe
\Deltay j
\Deltay j \Deltay j \Gamma1=2
\Deltat n
\Deltay j
\Deltat n
\Deltay j
\Deltay j \Deltay j+1=2
\Deltay j
\Deltay j \Deltay j \Gamma1=2
\Deltat n
\Deltat n
\Deltay j
\Deltay j \Deltay j+1=2
ae
\Deltat n
\Deltat n
\Deltay j
\Deltat n
\Deltat n
\Deltat n
\Deltay j
D\Deltat n
\Deltax i+1=2
D\Deltat n
\Deltay j
\Deltay j+1=2
\Deltay j \Gamma1=2
oe
\Deltay j
\Deltay j \Deltay j \Gamma1=2
\Deltat n
\Deltay j
\Deltat n
\Deltay j
\Deltat n (K n
and for
nx ;j =u n+1
ae
\Deltax nx
\Deltay j
\Deltat n
\Deltay j
D\Deltat n
\Deltax nx \Deltax nx+1=2
D\Deltat n
\Deltay j \Deltay j+1=2
D\Deltat n
\Deltay j \Deltay j \Gamma1=2
oe
\Deltax nx
\Deltax nx \Deltax nx \Gamma1=2
\Deltat n
\Deltay j
\Deltay j \Deltay j+1=2
\Deltay j
\Deltay j \Deltay j \Gamma1=2
\Deltax nx
1;nx \Gamma1=2;j
\Deltat n
\Deltay j
1;nx ;j \Gamma1=2
\Deltat n
\Deltay j
\Deltay j \Deltay j+1=2
\Deltay j
\Deltay j \Deltay j \Gamma1=2
Similar expressions hold for u n+1
i;j .
Proof. From (2.2), we see that
\Deltat n
\Deltat n
\Deltay j
\Deltay j \Deltay j+1=2
\Deltay j
\Deltay j \Deltay j \Gamma1=2
\Deltat n
\Deltat n
\Deltay j
Then, the proof is completed by simple algebraic manipulations on u n+1
i;j .
Lemma 5.2. Assume that
and, for
(5.1b) \Deltat n -D ?
=\Deltay ?
Then there is a constant C
ny
\Deltay j+1=2
ny
\Deltay j+1=2
ny
ny
Proof. From (5.1b) we see that the coefficients of the terms between the brackets
fg in the expressions of Lemma 5.1 are nonnegative. Then the estimate of a typical
term is given as follows:
\Deltat n
\Deltat n
\Deltay j
\Deltay j \Deltay j+1=2
ae
\Deltat n
\Deltat n
\Deltay j
\Deltat n
\Deltat n
\Deltat n
\Deltay j
D\Deltat n
\Deltax i+1=2
D\Deltat n
\Deltay j
\Deltay j+1=2
+\Deltay j \Gamma1=2
oe
\Deltay j
\Deltay j \Deltay j \Gamma1=2
\Deltat n
\Deltay j
\Deltat n
\Deltay j
\Deltat n jK n
Thus simple algebraic manipulations and use of (2.3c) yield the desired result.
We are now in a position to prove Theorem 3.2.
Proof of Theorem 3.2. Note that the CFL condition (3.6) implies (5.1b) by (2.3a).
Then the result (3.7) follows by iterating on n the inequality in Lemma 5.2 and using
Theorem 3.1 and (2.3b).
6. Proof of continuity with respect to data. In this section we prove Theorem
3.3 and a result on equicontinuity in time of the approximate solution, Proposition
6.4 below. We recall that v h stands for the solution of (2.2) with the data v D and
v init .
Lemma 6.1. For
where A n
i\Gamma1;j , and E n
i;j are defined as in the proof of Lemma
4.1.
The result easily follows from (2.2).
Lemma 6.2. Supposed that (4.1) is satisfied. Then
+\Deltay ?
\Theta jju D \Gamma v D jj L 1
Proof. Since, by (4.1), the coefficients in the equality of Lemma 6.1 are nonneg-
ative, we see that
Then multiplying by \Deltax i \Deltay j , adding over i, j, and rearranging terms imply the
desired result.
Now Theorem 3.3 can be easily seen from Lemma 6.2.
Lemma 6.3. Assume that the CFL condition (4.1) is satisfied. Then
+\Deltay ?
Proof. By (2.2), we observe that
\Deltat 0
\Deltat 0
\Deltay j
\Deltay j \Deltay j+1=2
\Deltay j
\Deltay j \Deltay j \Gamma1=2
\Deltay j
Then the lemma follows by multiplying this inequality by \Deltax i \Deltay j and adding the
resulting one over i, j.
Proposition 6.4 (Equicontinuity in time). Under the assumptions of Theorem
3.1, there is a constant C 4 depending only on the data and T such that for
+\Deltay ?
\Theta
Proof. Take v n+1
h in Lemma 6.2 and use Lemma 6.3 to obtain the result

7. A convergence analysis. In this section we prove Theorem 3.4 by applying
the ideas used in [10] for analyzing the one-dimensional drift-diffusion semiconductor
device equations. We point out that the analysis here is much simpler than that
given in [10]. The reason is that we are here using the standard entropy j \Delta j, while
a smoother entropy has been used there, which requires muck work to estimate the
distance between the smooth entropy and the standard one. We also emphasize the
difference between the present analysis and that used in classical conservation laws; in
the present case the delicate part is how to handle the boundary terms in the 'entropy
form' \Theta (see (7.3) below), while an unbounded domain is treated in the classical
conservation laws.
It should be emphasized that this whole section concerns the case of
that, although the differential equation (1.1) is linear, techniques which have been
originally developed for nonlinear hyperbolic conservation laws will be used.
The proof of Theorem 3.4 proceeds as follows. First, we prove that there is a
subsequence fu h 0 g h 0 ?0 converging to a limit u. Then, we show that
lim
(7.1a)
\Theta\Omega\Gamma2 where R(\Delta; \Delta) defines the left-hand side of (3.10). Since the
weak solution of (3.10) is assumed unique, this completes the proof of Theorem 3.4.
As in classical conservation laws, (7.1) follows from the following result [13]:
lim
(7.2a)
where \Theta is defined in (7.3) below. Most part of this section is devoted to proving this
result.
7.1. The entropy form. The entropy form \Theta(u; c; V ; ') with boundary terms
included is defined as follows:
ft=0g\Theta\Omega
and the `entropy flux' G and the function H are defined
by
The motivation of the form \Theta can be given as in the one-dimensional case [10].
7.2. A convergent subsequence. In this subsection we prove the existence of
a convergent subsequence fu h 0 g h 0 ?0 .
Lemma 7.1. Assume that the hypotheses of Theorem 3.2 are satisfied. Then
there exists a subsequence fu h 0 g h 0 ?0 converging in L 1
1(\Omega\Gamma9 to a limit u in
1(\Omega\Gamma3 .
Proof. We note that the ideas in [13] can be used to prove the lemma. In [13], a
discrete version of Azcoli-Arzel'a Theorem was used. In the present case with
the equicontinuity in time is provided by Proposition 6.4, and the compactness of
the range is given by Theorem 3.2. Also, the regularity result on u follows from the
convergence and Theorem 3.2.
7.3. Proof of (7.2a). Here we prove (7.2a) under a condition.
Lemma 7.2. Suppose that for c 2 IR and nonnegative
Then
Proof. First, for every nonnegative ' 2 C 1
\Theta\Omega\Gamma3 (7.5) follows from Lemma
7.1 and the standard argument in the classical conservation laws [13]. Also, since
1(\Omega\Gamma7 by Lemma 7.1, the same result holds for
We now consider the case where ' 2 C 1
(J \Theta [0; 1) \Theta (0; 1)). Since we are mainly
concerned with the boundary term associated with the edge
suffices to consider '(t; x; y) of this form !(t; y)-(x). Then, set
Z TZ 1ae
oe
and rewrite \Theta(u h
@t
QT
@y
QT
Since the sequence fu h 0 (\Delta; 0+; \Delta)g h 0 ?0 is bounded in L 1 (J \Theta (0; 1)) by Theorem 3.1,
there is a subsequence fu h 00 (\Delta; 0+; \Delta)g h 00 ?0 converging in L 1 (J \Theta (0; 1))-weak ? to a
limit e u. Let fl t;y be the Young measure corresponding to e u. Then, by Lemma 7.1 and
(2.3a), we see that
lim
@t
QT
@y
QT
where
(7.7a)
Z TZ 1ae
oe
where
Z
(7.8a)
Thus, to prove (7.5), it suffices to prove that
Z TZ 1ae
oe
!dtdy:
Take - such that its support is contained in [0; ffl]. Then, by (7.4) and Theorem 3.1, it
follows from (7.6) and (7.7) that
Since ffl is arbitrary, this inequality yields that
Choose c 2 IR such that ju \Gamma c) for some ff 2 IR. Then, by (7.8) and the
definition of fl t;y ,
so that, by (7.7), (7.10) becomes
ff
aeZ TZ 1(u(t; 0+;
oe
0:
Since the sign of ff is arbitrary and this inequality is true for any nonnegative ! 2
1)), we have
a.e. in J \Theta (0; 1);
(7.11a)
Finally, by (7.7)-(7.9) and (7.11), we see that
Z TZ 1/fi fi fi fi fi
Z
which together with (7.10) implies that This completes the proof of the case
1)). The same argument applies to the remaining three
cases.
It is now clear that it suffices to prove (7.4). This is done in the next two
subsections.
7.4. A discrete entropy inequality. The following discrete entropy inequality
will be needed for obtaining an upper bound for \Theta(u h
Lemma 7.3. Under the CFL condition (3.2), we have, for c 2 IR,
i;j \Deltat n )ju n
\Deltat n
\Deltat n
\Deltay j
where
Proof. From (2.2) and the definition of the mixed finite element space Q h , we
have, for c 2 IR,
\Deltat n
\Deltat n
\Deltay j
ae
\Deltat n
\Deltat n
\Deltay j
i;j \Deltat n
oe
\Deltay j
Note that the term between the brackets is nonnegative by (3.2). Thus the lemma
follows by multiplying this expression by sign (u n+1
7.5. An upper bound of entropy form. In this section we obtain an upper
bound for \Theta(u h which implies the inequality (7.4). We first have the
following decomposition of \Theta(u h
Let, for
Z y j+1=2
Z y j+1=2
Z y j+1=2
Z y j+1=2
Lemma 7.4. (decomposition of \Theta). We have
ent
where (with arguments omitted)
\Theta ent =
nx
ny
ae
i;j \Deltat n )ju n
\Deltat n
\Deltat n
\Deltay j
oe
and
nx
ny
ae
oe
\Deltat n \Deltay j
nx
ny
ae
oe
\Deltat n \Deltax i
nx
ny
nx
ny
nx
ny
nx
ny
i;j \Deltat n \Deltax i \Deltay j
nx
ny
Proof. From the definition of \Theta and the fact that div V h is piecewise constant, we
have
where
nx
ny
nx
ny
nx
ny
nx
ny
\Deltat n \Deltay j
nx
ny
\Deltat n \Deltax i
ny
nx+1=2;j \Deltat n \Deltay
ny
1=2;j \Deltat n \Deltay j
nx
i;ny +1=2 ' n+1=2
i;ny +1=2 \Deltat n \Deltax
nx
nx
ny
i;j \Deltat n \Deltax i \Deltay j
nx
ny
i;j \Deltat n \Deltax i \Deltay j
nx
ny
Then simple algebraic manipulations yield the desired result.
Lemma 7.5. (Upper bound of \Theta). Suppose that the conditions of Theorem
3.2 are satisfied. Then there is a constant C 5 depending solely on the data and T such
that for any ' 2 C 1 (Q T
\Theta ent -0;
\Theta com -C 5 (1
ae
oe
Proof. The first inequality follows immediately from Lemmas 7.3 and 7.4. Also,
observe that
\Deltay
\Deltat
Then, if an integration by parts on n is applied to the last two terms in the expression
of \Theta com , the second inequality follows from Theorems 3.1 and 3.2, Proposition 6.4,
Lemma 7.4, and (2.3a).
We are now in a position to prove Theorem 3.4.
Proof of Theorem 3.4. From Lemma 7.1 there exists a subsequence fu h 0 g h 0 ?0
converging in L 1
1(\Omega\Gamma2 to a limit u. Now, by Lemma 7.5, we have
lim
for every c 2 IR and nonnegative ' 2 C 1 (Q T ). Thus, by Lemma 7.2, we see that
lim
This implies that u is the unique solution of (3.10). Consequently, the whole sequence
converges to u, and thus Theorem 3.4 is proven.
8. Numerical results. This section reports on numerical results with the finite
element method (2.2) for three problems. They are designed to show the performance
of the method and to indicate the convergence properties. In all examples the CFL
condition (3.6) is required to hold.
FIG. 1. The '-', and exact u
and approximate solution u h with
Example 1. In this example we consider a convecting Gaussian hill in one space
dimension. Specifically, we solve equation (1.1) with
on the interval [0; 6]. The initial datum u init is given by
As a pure initial-value problem, this leads to the analytical solution
u a (t;
We obtain an initial-boundary-value problem with the same solution by imposing the
Dirichlet boundary condition:
In

Figure

1 we display the analytical solution u a and the approximate solution
u h at time 0:25. In Table 1 we display the errors and their respective order of
convergence at the same time. From the table we see that the scheme is first-order
accurate both in L 1 and in L 1 for the concentration. This shows that the scheme
(2.2) is first-order accurate in both spaces when the solution of the differential equation
is smooth. Also, Figure 1 agrees with the stability property given in Theorem 3.1.
Finally, our numerical experiments report (not shown here) that, if the CFL condition
(3.2) is violated, then the stability result (3.5) and the TVB boundedness (3.7) are no
longer valid.
100 3.57 0.76 2.50 1.03
200 1.89 0.91 1.20 0.99

Table

1. Convergence of u h in (0; 6) at
Example 2. In this example we consider a problem whose solution displays a
discontinuity. The problem has the data: The boundary
and initial conditions are given by
The exact and approximate 'nonviscous' solution (i.e., in the case of and the
'viscous' solution with are displayed in Figure 2. Notice that the
biggest error in the approximation of u occurs around the location of the discontinuity
In

Table

2 we show the errors and their convergence orders with
Note that the orders of convergence in L 1 and L 1 are nearly 1=2. This implies that
the presence of discontinuity has an effect on the convergence. Finally, from Figure 2,
we see that the 'nonviscous' solution is quite close to the `viscous' solution.
100 6.15 0.48 8.09 0.26
200 5.02 0.29 6.32 0.37
400 3.76 0.43 4.36 0.54

Table

2. Convergence of u h in (0; 1) at
FIG. 2. The '-' ,
and the approximate solution u h with
Example 3. In the third example we consider a two-dimensional problem which
has a shock. The data are set as
. The Neumann and Dirichlet boundaries
@\Omega 2 are
defined by
and the boundary and initial data by
ae
ae
elsewhere:
The approximate solution of this problem obtained using the method (2.2) with
shown in Figure 3. The graph clearly shows that the
method can capture the shock around the location
FIG. 3. The approximate solution u h on (0; 1) 2 .
Example 4. In the final example we test the sharpness of the bounds appearing
in (3.7) and (3.8) when D 6= 0 and jju D jj L 1 (J;BV
The same set of data are
chosen as in Example 3 except that the initial and boundary data are determined by
uD
where
Uniform partitions
of\Omega into rectangles are exploited. The TVB bounds on different
meshes at are given in Table 3. From this table we see that the left-hand side
of the inequality (3.7) blows up as \Deltay converges to zero. Similar results
are observed for the bound in (3.8) (not shown here).

Table

3. TVB Bounds of u h
in\Omega at
9. A concluding remark. A new finite element method for numerically solving
the two-dimensional convection-dominated transport equation in ground water has
been formulated and analyzed in this paper. The primary computational advantage
of the method is that it is local and thus fully parallelizable, and is conservative. The
stability properties of this method and its convergence in a suitable topology have
been established. Moreover, the numerical results have shown that the method is
first-order accurate when the solution is smooth and is one-half order accurate when
the solution has discontinuities, and that the method is non-oscillatory and shock-
capturing. Future work will be devoted to obtaining error estimates for both cases of
zero and nonzero coefficient D.

Acknowledgments

. The authors wish to thank the referees for their comments
leading to an improved presentation of this paper.



--R


Approximate symmetrization and Petrov-Galerkin methods for diffusion-convection problems
Finite element analysis of the 1D full drift diffusion semiconductor model
Numerical simulation of multiphase flow in groundwater hydrology
estimates for a finite element method for the drift-diffusion semiconductor device equations
Convergence of a finite element method for the drift-diffusion semiconductor device equations: the multidimensional case


Finite element methods for second order differential equations with significant first derivatives
Convergence of a finite element method for the drift diffusion semiconductor device equations
estimates for a finite element method for the drift diffusion semiconductor device equations: the zero-diffusion case
TVB Runge-Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws II: general framework
Monotone difference approximations for scalar conservation laws


Research Frontiers in Applied Mathematics
Multistep Galerkin methods along characteristics for convection-diffusion problems
Convergence analysis of an approximation of miscible displacement in porous media by mixed finite elements and a modified method of characteristics
Finite Volume
Multilevel Adaptive Methods for Partial Differential Equations
A mixed finite element method for second order elliptic problems
Finite element and finite difference methods for continuous flows in porous media
--TR
conservation law;finite element and volume method;convergence;mixed method;convection-diffusion-reaction equation;stability