A Mixed Volume Element With Characteristic Mixed Volume Element for Contamination Transport Problem

Nonlinear systems of convection-dominated diffusion equations are used as the mathematical model of contamination transport problem which is an important topic in environ mental protection science. An elliptic equation defines the pressure, a convection-diffusion equation expresses the concentration of contamination, and an ordinary differential equation interprets the surface absorption concentration. The transport pressure appears in the equation of the concentration which determines the Darcy velocity and also controls the physical process. The method of conservative mixed volume element is used to solve the flow equation which improves the computational accuracy of Darcy velocity by one order. We use the mixed volume element with the characteristic to approximate the concentration. This method of characteristic not only preserves the strong computational stability at sharp front, but also eliminates numerical dispersion and nonphysical oscillation. In the present scheme, we could adopt a large step without losing accuracy. The diffusion is approximated by the mixed volume element. The concentration and its adjoint vector function are obtained simultaneously, and the locally conservative law is preserved. An optimal second order estimates in l2-norm is derived.


Introduction
Numerical simulation of the contamination transport problem is an efficient way to find how to protect the environment and treat the pollution. The physical interpretation is given by a mathematical system of partial differential equations. In this model, an elliptic equation defines the pressure, a convection-diffusion equation expresses the concentration of contamination, and an ODE interprets the surface absorption concentration. We provide a simplified physical background in the following. When the chemical substances are dissolved in the underground water, the substances are absorbed by the surface of the media, and then the chemical reaction takes place. The chemical reaction also affects the transportation of the solutes as the underground water moves. Consideration of the reaction helps us to understand how the contamination transport in terms of time and space. The reaction speed depends on the absorption. The equilibrium is donoted quick reaction, and non-equilibrium absorption is denoted by slow reaction speed. In this paper, we discuss the non-equilibrium absorption. The following mathematical model is discussed with the initial-boundary value conditions (Dawson, Van Duiji & Wheeler, 1994;Dawson, 1998;Hornung, 1988;Vogt, 1982): ∇ · u = −∇ · κ(c)∇p = q(X, t), X = (x, y, z) T ∈ Ω, t ∈ J = (0, T ], (1a) u = −κ(c)∇p, X ∈ Ω, t ∈ J, θ(X) ∂c ∂t + ρ(X) ∂s ∂t + ∇ · (uc) − ∇ · (θD∇c) = 0, X ∈ Ω, t ∈ J, where Ω is a bounded domain in R 3 . p(X, t) denotes the pressure, and u = (u x , u y , u z ) T is Darcy velocity. The functions c(X, t) and s(X, t) denote the contamination concentration and non-equilibrium absorption concentration. κ(c) is the permeability. θ(X) and ρ(X), two positive functions, are water storage rates of moving water and the density, respectively. D D mol (X) is a molecular diffusion matrix and I is a 3 × 3 identity matrix. α l and α t respectively represent longitudinal and transverse diffusion coefficients.û x ,û y andû z are three direction cosines of Darcy velocity u. Generally speaking, β is a constant not less than 2, and diffusion-dispersion matrix is positive definite. Here we only consider the molecular diffusion. Suppose that D * ≤ D mol (X) ≤ D * holds for positive constants D * and D * .α, a positive constant, is the exchange coefficient. φ(c) denotes the isothermal absorption, generally defined by two different methods (I) Langmuir isotherm, φ(c) = k 1 c 1+k 2 c , k 1 , k 2 > 0, (II) Freundlinch isotherm, φ(c) = k 3 c p , k 3 > 0, where p ∈ [0, 1].
Assume that the mixture is impermeable across the boundary, where ν denotes the outer normal vector at the boundary ∂Ω of Ω.
The initial conditions are defined by The following restriction is introduced to confirm the existence and uniqueness Ω q(X, t)dX = 0, The Lipschitz continuity does not hold for this physical model, and this nonlinearity decreases the regularity. It brings trouble on numerical analysis for the full discrete scheme. Dawson first presents and analyzes the method of characteristics-Galerkin for non-equilibrium absorption (Dawson, Van Duiji & Wheeler, 1994;Dawson, 1998). While numerical oscillation possibly appears and the local conservation of mass is lost. Furthermore, an upwind with mixed finite element method is put forward to discuss a nonlinear contamination transport equation with equilibrium absorption in one variable. Convergence analysis is shown only for a local low-order semidiscrete algorithm. It is well-known that standard finite element could not solve convection-diffusion problems well because of strong numerical oscillation. In order to treat this problem well, some efficient numerical methods are put forward such as characteristic difference, characteristic finite element (Ewing, 1983;Yuan, 1996Yuan, ,1999Yuan, ,2013, upstream-weighted finite difference (Todd, ODell & Hirasaki, 1972), high-order Godunov scheme (Bell, Dawson & Shubin, 1988), streamline diffusion method (Johnson, 1986), least squares mixed finite element (Yang, 2000), modified method of characteristics with Galerkin finite element (MMOC-Galerkin) (Dawson, Russell & Wheeler, 1989), Eulerian-Lagrangian localized adjoint method (ELLAM) ( Cella, Russell, Herrera & Ewing, 1990) and so on. This problem should consider an important physical conservation. Arbogast and Wheeler present a type of characteristic mixed finite element (Arbogast & Wheeler, 1995) to solve the convection-diffusion equation efficiently, where MMOC-Galerkin is conservative locally on every elements and a 3/2 order error estimates is derived. Lots of mapping integrals of test functions make the computations more complicated. Therefore, we develop and improve substantially the work of Arbogast and Wheeler (Arbogast & Wheeler, 1995), and put forward a mixed finite element-characteristic mixed element method. This method decreases computational work greatly, and the feasibility and effectiveness are tested experimentally (Sun & Yuan, 2009). This method could not be generalized for three-dimensional problems. Finite volume method (Cai, 1991;Li & Chen, 1994) is carried out simply, and it has high accuracy and conservation. Thus it is an effective tool for solving partial differential equations. The method of mixed finite element (Douglas, Ewing & Wheeler, 1983, 19832;Raviart & Thomas, 1977) can solve the pressure and Darcy velocity simultaneously, and develops the accuracy. Combined the above two methods, a mixed finite volume method is discussed (Russell, 1995;Vogt, 1982;Weiser & Wheeler, 1988) and numerical experiments are argued (Cai, Jones, Mccormilk & Russell, 1997;Jones, 1995). Theoretical analysis is given for elliptic problem and a general discussion frame is given for mixed finite volume method (Chou, Kawk & Vassileviki, 1998, 2000Chou & Vassileviki, 1999). Rui and Pan use this method to discuss numerical computation for hypotonic oil-gas flow problems (Pan & Rui, 2012, 20122). In this paper, a type of mixed volume element with characteristic mixed volume element method is put forward for three-dimensional underground water pollution problem. The pressure and Darcy velocity are computed simultaneously by a conservative mixed volume element and the accuracy of Darcy velocity is improved. A conservative characteristic mixed finite volume element is used to compute the concentration, where the convection is approximated along the characteristics and the diffusion is discretized by the mixed volume element. The method of characteristics obtains the values at the sharp fronts very well.
In actual computations, a large time step may be adopted and the efficiency is improved. Piecewise-defined constants are taken as test functions, and the local conservation of mass is preserved. Finally, we obtain an optimal order estimates in L 2 norm superior to 3/2-order (Arbogast & Wheeler, 1995). Thus, the potential is shown to be an efficient tool in solving some actual applications (Arbogast & Wheeler, 1995;Ewing, 1983;Shen, Liu & Tang, 2002;Yuan, 2013).
Some symbols and assumptions are introduced. The regularity conditions of (1)-(7) are defined by Suppose that the coefficients satisfy the following positive definite conditions where k * , k * , θ * , θ * , ρ * , ρ * , D * and D * are positive constants.
In this paper, let K and ε denote a generic positive constant and a generic small positive number, respectively. They may have different definitions at different places.

The Preparations
Two different partitions are supposed to be regular for showing the numerical composite scheme. The first partition with a large step is nonuniform for solving the pressure and Darcy velocity, and the second nonuniform one is for the concentration. For simplicity, take Ω = {[0, 1]} 3 and let ∂Ω denote the boundary. The first partition is defined by δ x ×δ y ×δ z , δ y : 0 = y 1/2 < y 3/2 < · · · < y N y −1/2 < y N y +1/2 = 1, (Weiser & Wheeler, 1988;Yuan, 2013).
The inner products and norms are defined as follows The difference operators and other notation are introduced, The properties are prepared for the following numerical analysis.
where M is a constant independent of q and h.
Another partition is obtained by refining the coarse partition of Ω = {[0, 1]} 3 uniformly. Generally, we take the step by Since the absorption concentration changes slow and stably, eq. (2) is discussed on the first partition. A first-order finite element space M h is constructed on a hexahedron element (Cialet, 1978;Jiang & Pang, 1979).

The Procedures
In order to use the mixed finite volume element, we rewrite eq.(1) in a normalized formulation For eq.(2), considering that the flow transports along the characteristic direction, thus we apply the method of characteristics to approximate the first-order hyperbolic term. The computational algorithm could use a large time step and has the strong stability and high order accuracy. Let ψ(X, u) = [θ 2 (X) (2) is rewritten as follows for applying the method of mixed volume element, where f (X, c) = −qc.
Let P, U, C, G and S denote the discrete solutions of p, u, c, g and s, respectively. The pressure and Darcy velocity are computed by the method of characteristic mixed volume element, The derivative along the characteristics in (13a) is approximated by the backward difference quotient Eq. (13) is discretized by the characteristic mixed volume element method whereĈ n = C n (X n ),X n = X − θ −1 U n+1 ∆t.
The adsorption concentration changes slow and stably with respect to t, so we obtain an explicit scheme to approximate eq.(3), Initial approximations are defined by where (C 0 ,G 0 ) is the elliptic projection of (c 0 , g 0 ), andS 0 is an L 2 projection of s 0 (see the definitions in the next section).
The composite procedures run as follows. First of all, the elliptic projection and the initial values are used to obtain {C 0 ,G 0 }, as the initial approximations of c 0 and g 0 = −θD∇c 0 , i.e. C 0 =C 0 , G 0 =G 0 . Use s 0 and the L 2 projection to get S 0 =S 0 . Secondly, S 1 is computed by (16). Then, {U 1 , P 1 } is computed by the method of conjugate gradient and mixed volume element (14). Similarly, {C 1 , G 1 } is obtained by (15). Next, S 2 , {U 2 , P 2 } and {C 2 , G 2 } are computed by (15), (14) and (16). The computations proceed repeatedly and all the numerical solutions are obtained. By the positive definite condition (C), the solutions exist and are unique.

The Conservation of Mass
Assume that the problem (1)-(7) has no source or sink, i.e., q ≡ 0, and assume that the boundary condition is impermeable.
Here we ignore the effect of adsorption. Thus, for simplicity we take l = 1, or suppose that two partitions are the same.
Theorem 1. If q ≡ 0, and the factor of surface absorption is ignored, then on every element J c ⊂ Ω, (15a) satisfies the conservation law Using the previous notation, we have Substituting (21) into (20), we could complete this proof.
Furthermore, we get the conservation of mass on the whole domain.
Theorem 2. Suppose q ≡ 0 and the boundary is impermeable. The factor of surface absorption is ignored. Then, Proof. Summing (19) on all the elements, we have Then, the proof is completed by using − i, j,k ∂Ω i jk G n+1 · γ J c dS = − ∂Ω G n+1 · γdS = 0.

Convergence Analysis
The elliptic projections are introduced first.Ũ ∈ V h andP ∈ S h are defined by where c denotes the exact solution of (1)- (7).
The L 2 finite element projectionS ∈ M h is defined by Suppose that (1)-(7) satisfies the positive definite condition (C), and suppose that exact solutions satisfy the regularity (R). From the discussions of Weiser and Wheeler (Weiser & Wheeler, 1988), we see that auxiliary functions {Ũ,P,G,C} of (24) and (25) exist and are unique.
Lemma 5 If the coefficients and exact solutions of (1)-(7) satisfy (C) and (R), then there exist two constantsC 1 > 0 and C 2 > 0 independent of h such that Under the same assumptions, using the property of L 2 projection (Cialet, 1978;Jiang & Pang, 1979), we have the following estimates.
The following theorem is concluded from (64), Lemma 5 and Lemma 6.

Modified Mixed Volume Element With Characteristic Volume Element and Its Numerical Analysis
In §3 and §4, we discuss the method of mixed volume element with characteristic volume element. While in actual applications, Darcy velocity changes more slow than the concentration with respect to t. Therefore, a large time step may be used for solving eq.(1). The interval J is partitioned by 0 = t 0 < t 1 < · · · < t L = T . Let ∆t m p = t m − t m−1 , and suppose that the partition is measured uniformly by ∆t m p = ∆t p , m ≥ 2 except ∆t 1 p . The concentration is computed by a small uniform step ∆t c = t n − t n−1 . Suppose that there exists a positive integer n corresponding to m such that the pressure and the concentration have the same nodes, t m = t n . Let j = ∆t p /∆t c and j 1 = ∆t 1 p /∆t c . Let ψ m (X) = ψ(X, t m ). Darcy velocity u n+1 at t n+1 , t m−1 < t n+1 ≤ t m , in (15) is defined as follows. The approximation is defined by a linear extrapolation of U m−1 and U m−2 for m ≥ 2 If m = 1, define EU n+1 = U 0 .
The procedures of (1) and (2) are formulated. Numerical solutions (U m , P m ) : (t 0 , t 1 , · · · , t L ) → S h × V h and (C n , G n , C n ) : k −1 (C x m )U x m , w x x + k −1 (C y m )U y m , w y y + k −1 (C z m )U z m , w z z − P m , D x w x + D y w y + D z w z For clarity, the previous m-norm is replaced bym-norm.
(1) and (2) are discretized with two different steps. It can shorten the computational work greatly without losing accuracy.
In a similar analysis to Theorem 3, we conclude the following result. ≤ M * * h 2 p + h 2 c + ∆t c + (∆t p ) 2 + (∆t 1 p ) 3/2 , where the constant M * * depends on p, c, s and their derivatives.

Conclusions and Discussions
Numerical approximation and theoretical analysis of three-dimensional contamination transport are considered in this paper. A type of mixed volume element-characteristic mixed volume element method is discussed. Several interesting conclusions are stated as follows.
(I) The method has the physical nature of conservation (local and whole), an important physical nature in numerical simulation of underground environmental science problems.
(II) The composite method combines a mixed volume element and the idea of characteristics, so it has strong stability and high accuracy. It is quite useful especially in large-scaled engineering computation on a three-dimensional complicated region.