SIAM J. NUMER. ANAL. Vol. 33, No. 4, pp. 1669-1687, August 1996 () 1996 Society for Industrial and Applied Mathematics 021 A NONLINEAR MIXED FINITE ELEMENT METHOD FOR A DEGENERATE PARABOLIC EQUATION ARISING IN FLOW IN POROUS MEDIA* TODD ARBOGASTr, MARY F. WHEELERt, AND NAI-YING ZHANG Abstract. We study a model nonlinear, degenerate, advection-diffusion equation having application in petroleum reservoir and groundwater aquifer simulation. The main difficulty is that the true solution is typically lacking in regularity; therefore, we consider the problem from the point of view of optimal approximation. Through time integration, we develop a mixed variational form that respects the known minimal regularity, and then we develop and analyze two versions of a mixed finite element approximation, a simpler semidiscrete (time-continuous) version and a fully discrete version. Our error bounds are optimal in the sense that all but one of the bounding terms reduce to standard approximation error. The exceptional term is a nonstandard approximation error term. We also consider our new formulation for the nondegenerate problem, showing the usual optimal L2-error bounds; moreover, superconvergence is obtained under special circumstances. Key words, mixed finite element, degenerate parabolic equation, nonlinear, error estimates, porous media AMS subject classifications. Primary 65M60, 65M12, 65M15; Secondary 35K65, 76S05 1. Introduction. Let I2 C Rd, d 1, 2, or 3, be a bounded domain with sufficiently smooth boundary 012, and let 0 < T < cxz and J (0, T ]. We develop and analyze a mixed finite element approximation to the following nonlinear advection-diffusion problem inu(x,t)" (1.1a) (1.1b) (1.1c) v. [VP(u) +/(P(u))] (/(u)), Ot u u9, (x, t) 012 J, u=uo, (x,t)S2{O}, (x,t) EX2 J, where P(u) P(x, t; u) is strictly monotone increasing in u for each (x, t) E 12 x J, ?,(u) ?,(x, t; u),uo uo(x, t),uo uo(x), fl(u) =/3(x, t; u) is avector, andot or(x, t) is a d x d symmetric matrix that is uniformly positive definite with respect to (x, t) 6 x ]. These functions are tacitly assumed to be smooth enough for our purposes. We concentrate on the case in which OP(u)/Ou Pu(u) may be zero for some values of u. Since VP(u) Pu(x, t; u)Vu + Vx P(x, t; u), (1.1) is degenerate parabolic. Let (., .) denote the L2(f2)-inner product (or sometimes the duality pairing) and its norm. Our main assumptions are that there is a constant Co > 0, independent of time, such that (A1) IIP(qg)- P(g02)ll 2 < Co(P(q91)- P(p2), q)l -992), for qgl, 992 G L2(,Q), and both fl and F are Lipschitz continuous" (A2) II(go) (q92)ll + II’(goa) 9/(g02)11 Collq01 q9211 for (/91, (/92 G L2(S’2 ). *Received by the editors April 27, 1994; accepted for publication (in revised form) October 31, 1994. This work was partially supported by the Center for Research on Parallel Computation through National Science Foundation Cooperative agreement CCR-8809615 and by the State of Texas Governor’s Energy Office through contract 1059 for the Geophysical Parallel Computation Project. Department of Computational and Applied Mathematics, Rice University, Houston, TX 77251-1892. Current address: Department of Mathematics, University of Texas, Austin, TX 78712 (arbogast@math.utexas.edu). Department of Computational and Applied Mathematics, Rice University, Houston, TX 77251-1892. Current address: Texas Institute for Computational and Applied Mathematics, Taylor Hall 2.400, University of Texas, Austin, TX 78712 (mfw@ticam.utexas.edu). Center for Research on Parallel Computation, P.O. Box 1892, Rice University, Houston, TX 77251-1892. Current address: CLever Sys. Inc., 1334 Stokley Way, Vienna, VA 22182. 1669