Let be a bounded domain in with Lipschitz boundary , where . Let T be a positive constant. In this paper, we will consider the following linear convection-diffusion-reaction equations: find such that
Here, is the convection vector field; is the reaction function; and are given scalar functions; is a constant and is the outward unit normal vector to ; denotes the diffusion matrix function, where is the space of symmetric matrices, such that
where are positive constants.
Typically, these Equations (1) express the general mathematical model incorporating different types of transport phenomena in engineering and applied sciences, such as the dispersal of a pollutant through a moving viscous medium (e.g., a river or the atmosphere,  ), currents in semiconductor devices  , and airflow past an aerofoil (see  , for example). When the diffusive term is smaller than the convective one, these equations are the so-called convection dominated problems (see  ).
In many practical convection-diffusion processes, convection essentially dominates diffusion (e.g., in some financial models  ), and although the governing differential equation is parabolic, it displays several characteristics of hyperbolic problems. When applied to these problems, standard finite element and finite difference methods usually exhibit some combination of nonphysical oscillation and excessive numerical dispersion   . It is therefore logical to design numerical procedures that incorporate the parabolic/hyperbolic nature of these problems. One such method is the modified method of characteristics (MMOC) which was first formulated for a scalar parabolic equation by J. Douglas and T. F. Russell in  and then extended by Russell  to nonlinear coupled systems in two and three spatial dimensions. Similar schemes had been defined by Pironneau  for the incompressible Navier-Stokes equations and by Süli  and Morton, Priestley, and Süli  for first-order hyperbolic equations, with the latter technique being referred to as the Euler characteristic Galerkin method. The intent of the method is to obtain accurate approximations to convection- dominated problems. Basically, in the modified method of characteristics, the time derivative and the convection term are combined as a directional derivative. In other words, the procedure involves time stepping along the characteristics, allowing us to use large, accurate time steps.
Mixed finite element method has been proven to be an effective numerical method for solving fluid problems. It has an advantage to approximate the unknown variable and its diffusive flux simultaneously. There are many research articles on this method (     ). An algorithm combining the mixed finite-element method and the modified method of characteristics was first applied to the miscible displacement problem in porous media by Ewing, Russell, and Wheeler  . Then, the scheme had been extended by Wheeler and Dawson  to advection-diffusion-reaction problems. Numerical results have verified that large, accurate time steps are possible, and sharp fronts have been resolved (without oscillations or numerical diffusion) by coarser grids that standard procedures can use.
Arbogast and Wheeler  defined a characteristics-mixed method to approximate the solution of an advection-dominated transport problem. It used a characteristic approximation that is similar to that of MMOC-Galerkin method to handle advection in time and a lowest order mixed finite element spatial approximation for diffusion term. Piecewise constants were in the space of test function, so mass is conserved element-by-element. It was proved finally that the method was optimally convergent to order 1 in time and at least suboptimally convergent to order 3/2 in space. In  , we have considered a combined numerical approximation for incompressible miscible displacement in porous media. Standard mixed finite element was used for Darcy velocity equation and a characteristics-mixed finite element method was presented for approximating the concentration equation. Characteristic approximation was applied to handle the convection term, and a lowest order mixed finite element spatial approximation was adopted to deal with the diffusion term. Thus, the scalar unknown concentration and the diffusive flux can be approximated simultaneously. This approximation conserves mass globally. The optimal L2-norm error estimates were derived. Then, we extended this method to the slightly compressible miscible displacement problem in  .
It should be pointed out that the works mentioned above which involved the characteristic method all gave one order accuracy in time increment . That is to say, the first order characteristic method in time was analyzed. As for higher order characteristic method in time, Rui and Tabata  used the second order Runge-Kutta method to approximate the material derivative term for convection-diffusion problems. The scheme presented was of second order accuracy in time increment , symmetric and unconditionally stable. Optimal error estimates were proved in the framework of L2-theory. Numerical analyses of convection-diffusion-reaction problems with higher order characteristic/finite elements were analyzed in   , which extended the work  . The error estimates of second order in time increment were obtained.
The goal of this paper is to present a second order characteristic mixed finite element method in time increment to handle the material derivative term of (1). It is organized as follows. In Section 2, we formulate an approximate scheme that combines the second order characteristic finite element method for the material derivative term and mixed finite element method for the diffusion term. In Section 3, we prove the stability of the combined approximate scheme. In Section 4, we derive the L2-norm error estimates for both the unknown scalar variable and its flux. The scheme is of second order accuracy in time increment, symmetric, and unconditionally stable.
2. Formulation of the Method
In this paper, we adopt notations and norms of usual Sobolev spaces. Moreover, we adopt some notations for the functional spaces involved, which were used in    . For a Banach space X and a positive integer m, spaces and are abbreviated as and , respectively, and endowed with norms
where denotes the j-th derivative of with respect to time. The Banach space is defined by
equipped with the norm Similar spaces
and norms are considered for the boundary sets ΓR and ΓD. Denote by the closed subspace of defined by
2.1. The Characteristic Lines
Now, we define the characteristic lines associated with vector field and recall some classical properties satisfied by them. Thus, for given , the characteristic line through is the vector function solving the initial value problem
Next, assuming they exist, we denote by (respectively, by L) the gradient of (respectively, of ) with respect to the space variable x, i.e.,
We adopted some propositions and lemmas from  .
Proposition 1. If for an integer, then and it is with respect to the x variable.
In order to compute second order approximations of matrices and , we need the following equations:
Proposition 2. If , then
Proposition 3. If , then satisfies the Taylor expansion
and its inverse, , satisfies the Liouville’s theorem
By using the Liouville’s theorem and the chain rule, we obtain
Proposition 4. If , then
Proposition 5. If , then satisfies
2.2. Variational Formulation
From the definition of the characteristic lines and by using the chain rule, it follows that
By introducing the flux and using (12), we rewrite equation (1.1) at point and time as follows
Before giving a week formulation of (13), we adopted a lemma from  , which can be considered as a Green’s formula.
Lemma 6. Let be an invertible vector valued function. Let and assume that Then
with a vector valued function and a scalar function.
Now, we can multiply (13) by test functions and , integrate in respecetively, and apply the usual Green’s formula and (14) with , obtaining
Lemma 7.  Let be an invertible vector valued function satifying . Let and assume that Then
with and , where is the outward unit normal vector to .
Now, replacing in (15) formula (16) with , and replacing the Robin condition, we have
2.3. The Combined Approximate Scheme
We now present our time-stepping procedure for (17). Let N be the number of time steps, be the time step, and for . We will use the notation for a function. Moreover, for , we define
By fixing in (17) and using a Crank-Nicolson method  with respect to . Thus, from (18) we have
By using (8) and (11), we see that
Taking (20) into (19), we can obtain
We propose two explicit numerical schemes to approximate :
A similar notation to the one in Section 2.2 is used for the Jacobian of , namely,
Now, we restate three lemmas concerning properties of the characteristic line approximations. For this, we require the time step to be bounded and the velocity to satisfy the following assumption.
Claim 1. The velocity field and satisfies on .
Lemma 8.  Under Claim 1, if , we can see that
Lemma 9.  Under Claim 1, if , we have
Corollary 1. Under the assumptions of Lemma 9, , we have
Lemma 10.  Under Claim 1, if and , then there exists a positive constant such that
Thus, in the case where the characteristic lines and their gradients are not explicitly known, we propose the following time approximation of (21)
The time difference approximation (27) will be combined with a standard Galerkin finite element and mixed finite element in the space for and , respectively   . We discrete in space on a quasi-uniform finite element mesh of with maximal element diameter . For , we denote as and similarly. Let be finite element spaces with index k and l, respectively.
We define a bilinear form on and a linear form on for by
Then, the fully discrete scheme reads: Given , find such that
Throughout the analysis, K will denote a generic positive constant, independent of . Similarly, will denote a generic small positive constant.
3. Stability of the Approximate Scheme
In this section, we derive the stability of the approximate scheme (29). In order to develop the stability, some assumptions on the different terms of (1) are required.
Claim 2. The velocity field and satisfies on .
Remark. Throughout this paper denotes the maximum between the positive constant appearing in Lemma 10 and the norm of the velocity in .
Claim 3. The diffusion matrix coefficients, , belong to . More- over, A is a positive definite symmetric matrix and there exists a strictly positive constant which is a uniform lower bound for the eigenvalues of .
As a consequence of Claim 3, there exists a unique positive definite symmetric matrix function , such that and . Let us in-
troduce the constant . Clearly, Claim 3, we have
Claim 4. The reaction function, , satisfies in , where is a constant.
Under the previous claims, let
Claim 5. The source function . In Robin boundary condition, and .
For a given series of functions , we define the following norms
Similar definitions are considered for functional spaces and associated to the Robin boundary condition. Moreover, we define
Lemma 11. Let the above Claims 2 - 5 and be assumed. If be the solution of (29). Then it holds that
the constant c is given by
Proof. Substituting into (28), we have
Lemma 4.1 in  implies that
Next, by using , we obtain
Then when and , we have
Similarly, for we obtain the estimate
Analogous computations to term give
For the boundary integral term , we first use some properties of the inner product in the space and the inequality to get the estimate
for . Thus, we obtain
Then, by summing up from (33) to (38), inequality (31) follows.
By using Lemma 11 and following the arguments to Lemmas 5.6, 5.7 and Theorem 5.8 in  , we can get the following stability theorem:
Theorem 12 (Stability Theorem) Let the above Claims 2 - 5 assumed, and be the solution of (29) subject to the initial value . Then there exist two positive constants c and , such that, for , we have
4. Error Estimate Theorem
Now, we turn to derive a priori error estimate in L2-norm for the solutions of (29). In order to state error estimates, we need two following Lagrange interpolation operators (   ) , and .
Lemma 13. There exist positive constants and , independent of and respectively, such that
Corresponding to and , we introduce a bilinear form on and a linear form on for as follows
If and are the solutions of (1), we have for
We decompose as
In order to estimate the terms on the right-hand side of (44), we adopt the following lemmas.
Lemma 14.  Assume the above Claims 2 - 5 hold, and that the coefficients of the problem (1) satisfy ,
and that . Let the solution of (43) satisfy , , . Then, for each there exist two functions and ,
. Moreover, , and the following estimates hold:
where denotes a constant independent of .
Lemma 15.  Assume the above Claims 2 - 5 hold, and that the coefficients of the problem (1) satisfy , and . Then, for each there
exist function and , such that
Moreover, and the following estimates hold:
where denotes a constant independent of .
Now, we turn to bound the third term on the right-hand side of (44).
Lemma 16. Assume the above Claims 2 - 5 hold, and that the coefficients of the problem (1) satisfy and . There exists
with , a positive constant.
Proof. From the definition of , we have
From the arguments of Lemma 4.1 in  , we get similarly the bounds of the above terms as the followings:
To estimate , we divide it into three parts
Using the transformation , we have
Now, we replace in equality , where
with , and the function is defined by Thus, we get
Moreover, we have
By considering together from (54) to (58), we can state
Similar arguments, i.e. Lemma 5.4 in  lead to
Using these inequalities, and can be bounded as follows
Gathering (50), (51), (52), (53), (59), and (61) together, we complete the proof.
We now turn to estimate and . From the definition of and , we see
From Lemma 11, we obtain the lower bound for (62)
By jointly considering the lower bound (63), the upper bounds given in Lemmas 14 - 16, we deduce
Analogous computations to those developed in Lemma 16 give
Putting (65) into (64), multiplying (64) by , summing it about time form and , taking the initial values and using discrete Gronwall’s inequality, we can derive the following error estimate:
Theorem 17 (Error Estimate) Let Claims 2 - 5 be assumed. Let be the solution of (1) with
, be the solution of (29) subject to the initial values
. Then there exists a positive constant independent of and such that
This research was supported by the NSF of China (Grant Nos. 11271231 and 11301300).
 Frolkovic, P. and De Schepper, H. (2000) Numerical Modelling of Convection Dominated Transport Coupled with Density Driven Flow in Porous Media. Advances in Water Resources, 24, 63-72.
 Markowich, P.A. and Szmolyan, P. (1989) A System of Convection-Diffusion Equations with Small Diffusion Coefficient Arising in Semiconductor Physics. Journal of Differential Equations, 81, 234-254.
 Ewing, R.E. and Wang, H. (2001) A Summary of Numerical Methods for Time-Dependent Advection-Fominated Partial Differential Equations. Journal of Computational and Applied Mathematics, 128, 423-445.
 Douglas, Jr.J. and Russell, T.F. (1982) Numerical Methods for Convection-Dominated Diffusion Problems Based on Combining the Method of Characteristics with Finite Element or Finite Difference Procedures. SIAM Journal on Numerical Analysis, 19, 871-885.
 Russell, T.F. (1985) Time-Stepping along Characteristics with Incomplete Iterationfora Galerkin Approximation of Miscible Displacement in Porous Media. SIAM Journal on Numerical Analysis, 22, 970-1013.
 Johnson, C. and Thomée, V. (1981) Error Estimates for Some Mixed Finite Element Methods for Parabolic Type Problems. RAIRO Analyse Numérique, 15, 41-78.
 Raviart, P.A. and Thomas, J.M. (1977) A Mixed Finite Element Method for Second Order Elliptic Problems, Mathematical Aspects of the Finite Element Method. Lecture Notes in Mathematics 606, Springer, Berlin.
 Douglas, Jr.J. and Roberts, J.E. (1985) Global Estimates for Mixed Methods for Second Order Elliptic Equations. Mathematics of Computation, 44, 39-52.
 Ewing, R.E., Russell, T.F. and Wheeler, M.F. (1984) Convergence Analysis of an Approximation of Miscible Displacement in Porous Media by Mixed Finite Elements and a Modified Method of Characteristics. Computer Methods in Applied Mechanics and Engineering, 47, 73-92.
 Arbogast, T. and Wheeler, M.F. (1995) A Charcteristics-Mixed Finite Element Methods for Advection-Dominated Transport Problems. SIAM Journal on Numerical Analysis, 32, 404-424.
 Sun, T.J. and Yuan, Y.R. (2009) An Approximation of Incompressible Miscible Displacement in Porous Media by Mixed Finite Element Method and Characteristics-mixed Finite Element Method. Journal of Computational and Applied Mathematics, 228, 391-411.
 Sun, T.J. and Yuan, Y.R. (2015) Mixed Finite Element Method and the Characteristics-mixed Finite Element Method for a Slightly Compressible Miscible Displacement Problem in Porous Media. Mathematics and Computers in Simulation, 107, 24-45. https://doi.org/10.1016/j.matcom.2014.07.005
 Bermúdez, A., Nogueiras, M.R. and Vázquez, C. (2006) Numerical Analysis of Convection-Diffusion-Reaction Problems with Higher Order Characteristic/Finite Elements, Part I: Time Discretization. SIAM Journal on Numerical Analysis, 44, 1829-1853.
 Bermúdez, A., Nogueiras, M.R. and Vázquez, C. (2006) Numerical Analysis of Convection-Diffusion-Reaction Problems with Higher Order Characteristic/Finite Elements, Part II: Fully Discretized Scheme and Quadrature Formulas. SIAM Journal on Numerical Analysis, 44, 1854-1876.
 Brezzi, F. (1974) On the Existence, Uniqueness and Approximation of Saddle-Point Problems Arising from Lagrangian Multipliers. RAIRO Analyse Numérique, 2, 129-151.
 Wheeler, M.F. (1973) A Priori L2 Error Estimates for Galerkin Approximations to Parabolic Partial Differential Equations, SIAM Journal on Numerical Analysis, 10, 723-759.
 Arnolds, D.N., Scott, L.R. and Vogelus, M. (1998) Regular Inversion of the Divergence Operator with Dirichlet Boundary Conditions on a Polygonal. Annali della Scuola Normale Superiore di Pisa-Classe di Scienze, 15, 169-192.
 Russell, T.F. (1986) Finite Elements with Characteristics for Two-component Incompressible Miscible Displacement, SPE 10500. Proceedings of 6th SPE Symposium on Reservoir Simulation, New Orleans, 123-135.