\documentclass[twoside]{article} \usepackage[dvips]{graphics} % To input PostScript figures \usepackage{amssymb} % used for R in Real numbers \pagestyle{myheadings} \markboth{\hfil Degenerate two-phase incompressible flow III \hfil}% {\hfil Zhangxin Chen \& Natalia L. Khlopina\hfil} \begin{document} \setcounter{page}{29} \title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent {\sc 15th Annual Conference of Applied Mathematics, Univ. of Central Oklahoma}, Electronic Journal of Differential Equations, Conference~02, 1999, pp. 29--46. \newline ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu \newline ftp ejde.math.swt.edu (login: ftp)} \vspace{\bigskipamount} \\ Degenerate two-phase incompressible flow problems III: Perturbation analysis and \\ numerical experiments \thanks{ {\em 1991 Mathematics Subject Classifications:} 35K60, 35K65, 76S05, 76T05. \hfil\break\indent {\em Key words and phrases:} porous medium, degenerate elliptic-parabolic system, \hfil\break\indent perturbation method, finite element, regularization, two-phase flow. \hfil\break\indent \copyright 1999 Southwest Texas State University and University of North Texas. \hfil\break\indent Partly supported by National Science Foundation grants DMS-9626179, DMS-9972147, and INT-9901498 and a gift grant from the Mobil Oil Company, Dallas, Texas. \hfil\break\indent Published November 24, 1999.} } \date{} \author{Zhangxin Chen \& Natalia L. Khlopina} \maketitle \begin{abstract} This is the third paper of a three-part series where we develop and analyze a finite element approximation for a degenerate elliptic-parabolic partial differential system which describes the flow of two incompressible, immiscible fluids in porous media. The approximation uses a mixed finite element method for the pressure equation and a Galerkin finite element method for the saturation equation. It is based on a regularization of the saturation equation. In the first paper \cite{RckA} we analyzed the regularized differential system and presented numerical results. In the second paper \cite{RckB} we obtained error estimates. In the present paper we describe a perturbation analysis for the saturation equation and numerical experiments for complementing this analysis. \end{abstract} \section{Introduction} The flow of two incompressible, immiscible fluids in a porous medium $\Omega\subset{\mathbb R}^d$, $d\leq3$ \cite{Rbear, Rpeaceman} is given by $$ \begin{array}{l@{}l} &\phi\partial_t s-\nabla\cdot(\kappa\lambda_{w}(s) (\nabla p_{w}+\gamma_{w}))=q_{w}, \\[3pt] &-\phi\partial_t s-\nabla\cdot(\kappa\lambda_{o}(s) (\nabla p_{o}+\gamma_{o}))=q_{o}, \\[3pt] &p_{c}(s)=p_{o}-p_{w}, \end{array} \eqno(1.1) $$ where $w$ indicates a wetting phase (e.g., water), $o$ denotes a nonwetting phase (e.g., oil), $\phi$ and $\kappa$ are the porosity and absolute permeability of the porous system, $s$ is the (reduced) saturation of the wetting phase, $p_{\alpha}$, $\lambda_{\alpha}$, $\gamma_{\alpha}$, and $q_{\alpha}$ are, respectively, the pressure, mobility (i.e., the relative permeability over the viscosity), gravity-density vector, and external volumetric flow rate of the $\alpha$-phase ($\alpha=w, o$), and $p_{c}$ is the capillary pressure function. To separate the saturation equation from the pressure equation, we define the total mobility $$ \lambda(x,s)=\lambda_{w}+\lambda_{o}. $$ Also, following \cite{Ranton, Rcj} we define the global pressure $$ p=p_{o}-\int^{s}_{0}\left(\frac{\lambda_{w}} {\lambda}\frac{\partial p_{c}}{\partial s}\right)(x,\xi) \mbox{d}\xi, \eqno(1.2) $$ and following \cite{RchenJDE} the complementary pressure $$ \theta=D(s)=-\int^{s}_{0}\left(\kappa\frac{\lambda_{w}\lambda_{o}}{\lambda} \frac{\partial p_{c}}{\partial s}\right)(x,\xi) \mbox{d}\xi. \eqno(1.3) $$ Then, by (1.2), (1.3), and some manipulations, it follows from (1.1) that \cite{RchenJDE, RceNM} $$ \begin{array}{l@{}l} &-\nabla \cdot\left\{\kappa(\lambda(s)\nabla p+\gamma^\prime_{1}(s))\right\} =q\equiv q_{w}+q_{o},\\[3pt] &\phi\partial_t s -\nabla\cdot\big\{\nabla \theta+\kappa\big(\lambda_{w}(s)\nabla p +\gamma^\prime_{2}(s)\big)\big\}=q_{w}, \end{array}\eqno(1.4) $$ where $$ \begin{array}{l@{}l} &\gamma^\prime_{1}=-\lambda_{w}\nabla_x p_{c} +\lambda\int^{s}_{0}\nabla_x\left(\frac{\lambda_{w}}{\lambda} \frac{\partial p_{c}}{\partial s}\right)(x,\xi) \mbox{d}\xi+\lambda_{w}\gamma_{w}+\lambda_{o}\gamma_{o},\\[3pt] &\gamma^\prime_{2}=-\lambda_{w}\nabla_x p_{c} +\lambda_{w}\int^{s}_{0}\nabla_x\left(\frac{\lambda_{w}}{\lambda} \frac{\partial p_{c}}{\partial s}\right)(x,\xi) \mbox{d}\xi+\lambda_{w}\gamma_{w}\\[3pt] &\qquad\qquad +\int^{s}_{0}\nabla_x\left(\frac{\lambda_{w}\lambda_{o}}{\lambda} \frac{\partial p_{c}}{\partial s}\right)(x,\xi) \mbox{d}\xi. \end{array} $$ In (1.4), $s$ is related to $\theta$ through (1.3): $$ s={{\mathcal S}}(\theta),\eqno(1.5) $$ where ${{\mathcal S}}(x,\theta)$ is the inverse of $D(s)$ for $0\leq \theta\leq \theta^*(x)$ with $$ \theta^*(x)=-\int^{1}_{0}\left(\kappa\frac{\lambda_{w}\lambda_{o}}{\lambda} \frac{\partial p_{c}}{\partial s}\right)(x,\xi) \mbox{d}\xi. $$ The pressure equation is given by the first equation of (1.4), while the saturation equation is described by the second equation. They determine the main unknowns $p$, $s$, and $\theta$. The model is completed by specifying boundary and initial conditions. For simplicity, in this paper we only consider the Neumann boundary conditions $$ \begin{array}{r@{}l@{}ll} &-\kappa(\lambda(s)\nabla p+\gamma^\prime_{1}(s)) \cdot\nu =\varphi_{1}(x,t), &\quad& (x,t)\in \Gamma\times J,\\[3pt] &-\big\{\nabla \theta+\kappa\big(\lambda_{w}(s)\nabla p +\gamma^\prime_{2}(s)\big)\big\} \cdot\nu =\varphi_{2}(x,t), &\quad& (x,t)\in \Gamma\times J,\\[3pt] \end{array} \eqno(1.6) $$ where $\varphi_1$ and $\varphi_2$ are given functions, $J=(0, T]$ ($T>0$), $\Gamma$ is the boundary of $\Omega$, and $\nu$ is the outer unit normal to $\Gamma$. The boundary conditions in (1.6) come from those imposed for the phase quantities via the transformations (1.2) and (1.3) \cite{RceSIAMNUMB}. For other types of boundary conditions, refer to \cite{RceNM}. The initial condition is given by $$ s(x,0)=s_{0}(x), \quad x\in\Omega. \eqno(1.7) $$ In recent years, the interest in the numerical simulation of two-phase fluid flow in porous media has been rising rapidly (see \cite{Rewing} and the bibliographies therein). In conjunction, there has been intensive research into the error analysis of numerical methods used in the simulation (see the extensive references in \cite{RceSIAMNUMB}). However, in most previous works the error analysis has been carried out under the unrealistic assumption that the capillary diffusion coefficient is uniformly positive. Namely, it has been assumed that $-({\lambda_{w}\lambda_{o}}/{\lambda})({\partial p_{c}}/{\partial s})$ is uniformly positive. The case where this diffusion coefficient can be zero has been treated in \cite{Rcee, RceSIAMNUMB, RfsA, RfsB, Rsmy}. But, in these papers only simplied saturation equations have been analyzed; i.e., the second equation in (1.4) with $p$ (or $u$ the total velocity; see (3.1a) later) given and $q_w=\gamma_2^\prime\equiv0$ has been considered \cite{Rcee, RceSIAMNUMB, RfsA, RfsB, Rsmy}. More recently, in \cite{RceNM} the fully coupled system (1.4) has been analyzed for the finite element approximation used here. It has been shown \cite{RceNM} that when this approximation directly solves the degenerate system, optimal error estimates behave like $O(h)$, where $h$ is the discretization mesh size. The main purposes of this series of three papers are to analyze regularized versions of this fully coupled system, improve the error estimates in \cite{RceNM}, describe a perturbation analysis, and present numerical experiments. The error analysis presented in \cite{RckB} was based on a regularization of the saturation equation. The diffusion coefficient of this equation was perturbed to obtain a nondegenerate problem with smooth solutions. The regularized solutions were shown to converge to the original solution as the perturbation parameter goes to zero with specific convergence rates given. Then a finite element approximation was used to solve the regularized solutions of the differential system. This approximation, which follows \cite{RceNM}, combined a mixed finite element method for the pressure equation and a Galerkin finite element method for the regularized saturation equation. Then, for this approximation we proved that the norm of error estimates depended on the severity of the degeneracy in diffusivity. In particular, for the degeneracy under consideration error estimates we obtained can behave like $O(h^{1+\epsilon})$, where $0<\epsilon<1$. This result thus improved that in \cite{RceNM} in some cases. In the first paper \cite{RckA} we analyzed the regularized differential system, described the finite element approximation, and presented numerical results. In the second paper \cite{RckB} we obtained error estimates. Both semidiscrete (continuous in time) and fully discrete approximations were analyzed. In the present paper we describe a perturbation analysis for the saturation equation and carry out numerical experiments for complementing this analysis. We remark that since the differential system for the single-phase, miscible displacement of one incompressible fluid by another in porous media resembles that for the two-phase incompressible flow studied here \cite{RceSIAMMA}, the analysis presented in this paper extends to the miscible displacement problem. Also, due to its convection-dominated feature, more efficient approximate procedures should be used to solve the saturation equation. Characteristic-based finite element methods will be considered in forthcoming papers. The rest of the paper is organized as follows. After preliminaries in the next section we consider a regularization of the differential system (1.4) in the third section and a perturbation method for the saturation equation in the forth section. Then we describe a finite element approximation for this regularized system in the fifth section. Numerical experiments are given in the final section. \section{Preliminaries} In this section we present preliminary results used in later sections. In particular, we collect some results for the Poisson solution operator $E$, which is defined below, and for the regularized diffusion coefficient. For $g\in H^{-1}(\Omega)$ (the dual to $H^1(\Omega)$), set $$ g_\Omega=(g,1) $$ in the sense of duality between $H^{-1}(\Omega)$ and $H^1(\Omega)$. When $g$ is Lebesgue integrable on $\Omega$, we have $$ g_\Omega=\int_\Omega g\mbox{d}x. $$ For $g\in H^{-1}(\Omega)$, consider the Neumann boundary value problem $$ \begin{array}{l@{}l@{}l} &-\Delta \omega =g-g_\Omega &\qquad\mbox{ in } \Omega,\\[3pt] &\nabla \omega \cdot\nu=0 &\qquad\mbox{ on } \Gamma,\\[3pt] &\omega_\Omega=g_\Omega. &\qquad{ } \end{array}\eqno(2.1) $$ The Poisson solution operator $E : H^{-1}(\Omega)\to H^1(\Omega)$ is defined by $$ E(g)=\omega. $$ It follows from (2.1) that $$ (\nabla \omega ,\nabla v)=(g, v)-g_\Omega v_\Omega, \qquad v\in H^1(\Omega); $$ i.e., $$ (g, v)=(\nabla(E(g)),\nabla v)+g_\Omega v_\Omega. $$ Especially, take $v=E(g)$ to see that $$ (g, E(g))=\|\nabla(E(g))\|^2_{L^2(\Omega)}+(E(g)_\Omega)^2. \eqno(2.2) $$ It can be easily checked \cite{RfsA} that the operator $E$ is linear, symmetric, and positive definite. Furthermore, for $v\in H^1(\Omega)$ note that the norm $$ \left(\|\nabla v\|_{L^2(\Omega)}^2+(v_\Omega)^2\right)^{1/2} $$ is equivalent to the usual norm on $H^1(\Omega)$. Thus, by (2.2), we can define the norm on $H^{-1}(\Omega)$ $$ \|g\|_{H^{-1}(\Omega)}=(g, E(g))^{1/2}, $$ which is equivalent to the usual norm on $H^{-1}(\Omega)$. As mentioned in the introduction we improve the error estimate in \cite{RceNM} by studying a regularized version of the saturation equation. For this, we define the capillary diffusion coefficient $$ d(s)=-\kappa\frac{\lambda_{w}\lambda_{o}}{\lambda} \frac{\partial p_{c}}{\partial s}, $$ and assume that it satisfies $$ d(s)\ge \left\{\begin{array}{l@{}ll@{}l} &c_1|s|^{\mu_1}, &\quad 0\leq s\leq\beta_1,\\[3pt] &c_2, &\quad \beta_1\leq s\leq \beta_2, \\[3pt] &c_3|1-s|^{\mu_2}, &\quad \beta_2\leq s\leq 1, \end{array}\right. \eqno(2.3) $$ where $c_i$ ($i=1,2,3$) are positive constants, and $\mu_i$ and $\beta_i$ ($i=1,2$) satisfy $$ 0<\beta_1<1/2<\beta_2<1, \qquad 0<\mu_1,\, \mu_2\leq 2. $$ Set $$ \mu=\max\{\mu_1, \mu_2\}, $$ and $$ \gamma=\frac{2+\mu}{1+\mu}. $$ Note that $\gamma$ is the conjugate to $2+\mu$. As in (1.3), we also set $$ \theta=D(s)=\int^s_0 d(\xi)\mbox{d}\xi. $$ Remark that (2.3) assumes the nature of degeneracy in the coefficient $d(s)$ near zero and one. The following lemma can be found in \cite{RfsA}. \medskip \noindent {\bf Lemma 2.1.} {\it Let $d$ satisfy} (2.3). {\it Then there exists a positive constant C, depending only on the parameters in} (2.3), {\it such that} $$ C(s_2-s_1)^{1+\mu}\leq D(s_2)-D(s_1), \qquad 0\leq s_1\leq s_2\leq 1.\eqno(2.4) $$ \section{Regularization} For the convenience of the later analysis, we rewrite (1.4) as follows: $$ \begin{array}{l@{}l@{}l} &\nabla \cdot u=q, \quad u=-\kappa(\lambda(s)\nabla p+\gamma^\prime_{1}(s)) &\quad\mbox{ in } \Omega_T, \\[3pt] &\phi\partial_t s -\nabla\cdot\big\{\nabla D(s)-f_{w}(s) u +\gamma_{2}(s)\big\}=q_{w} &\quad\mbox{ in } \Omega_T, \end{array}\eqno(3.1a) $$ where $\Omega_T=\Omega\times J$ and $$ f_w(s)=\lambda_w(s)/\lambda(s), \quad \gamma_2(s)=\kappa\{\gamma^\prime_2(s)-f_w(s)\gamma^\prime_1(s)\}. $$ The boundary and initial conditions become $$ \begin{array}{r@{}l@{}ll} &u\cdot\nu =\varphi_{1}(x,t), &\quad& (x,t)\in \Gamma\times J,\\[3pt] &-\big(\nabla D(s)-f_{w}(s) u +\gamma_{2}(s)\big) \cdot\nu =\varphi_{2}(x,t), &\quad& (x,t)\in \Gamma\times J,\\[3pt] &s(x,0)=s_0(x), &\quad& x\in \Omega. \end{array} \eqno(3.1b) $$ Existence and uniqueness of a solution to (3.1) in the weak sense has been shown in \cite{RchenJDE} with $$ \theta=D(s)\in L^2(J;H^1(\Omega)), \quad s\in L^\infty(\Omega_T), \quad p\in L^\infty(J; \hat V), $$ where $$ \hat V=\{v\in H^1(\Omega) : v_\Omega=0\}. $$ Also, it was shown under physically reasonable assumptions that $u$ is bounded: $$ \|u\|_{L^\infty(\Omega_T)}\leq C. \eqno(3.2) $$ Property (3.2) and the following assumptions are implicitly used in this paper: $\phi\in L^\infty(\Omega)$ satisfies that $\phi(x) \ge \phi_*>0$, $\kappa(x)$ is a bounded, symmetric, and uniformly positive definite matrix, i.e., $$ 0<\kappa_*\leq |\xi|^{-2}\sum_{i,j=1}^d\kappa_{ij}(x)\xi_i\xi_j \leq \kappa^*<\infty,\quad x\in\Omega, \, \xi\neq 0\in{\mathbb R}^d, $$ and $\lambda(s)$ satisfies that $$ 0<\lambda_*\leq \lambda(s)\leq \lambda^*<\infty, \qquad s\in [0,1]. $$ Without loss of generality, we further assume that $\phi\equiv1$ (otherwise, we consider the new variable $\hat s=\phi s$ instead of $s$ and the subsequent analysis is the same). Also, the functions $d$, $f_w$, $\gamma_1^\prime$, and $\gamma_2$ are assumed to be bounded functions of $s$. Finally, all the functions of $s$ need to be defined only on $[0,1]$. We replace $d$ by a positive $d_\beta >0$ with $d_\beta\to d$ in some sense as $\beta\to0$; a specific example of $d_\beta$ will be given at the end of this section. For given $d_\beta >0$, define $$ D_\beta(s)=\int^s_0 d_\beta(\xi)\mbox{d}\xi. $$ The corresponding non-degenerate differential system is given by $$ \begin{array}{l@{}l@{}l} &\nabla \cdot u_\beta=q, \quad u_\beta=-(\lambda(s_\beta)\nabla p_\beta+\gamma^\prime_{1}(s_\beta)) &\quad\mbox{ in } \Omega_T, \\[3pt] &\partial_t s_\beta -\nabla\cdot\big\{\nabla D_\beta(s_\beta)-f_{w}(s_\beta) u_\beta +\gamma_{2}(s_\beta)\big\}=q_{w} &\quad\mbox{ in } \Omega_T, \end{array}\eqno(3.3a) $$ with the boundary and initial conditions $$ \begin{array}{r@{}l@{}ll} &u_\beta\cdot\nu =\varphi_{1}(x,t), &\quad& (x,t)\in \Gamma\times J,\\[3pt] &-\big(\nabla D_\beta(s_\beta)-f_{w}(s_\beta) u_\beta +\gamma_{2}(s_\beta)\big) \cdot\nu =\varphi_{2}(x,t), &\quad& (x,t)\in \Gamma\times J,\\[3pt] &s_\beta(x,0)=s_0(x), &\quad& x\in \Omega. \end{array} \eqno(3.3b) $$ We now determine in what manner $(p_\beta, u_\beta, s_\beta)\to (p, u,s)$ as $\beta\to0$. Toward that end, we make the assumption $$ \begin{array}{l@{}l} &\|\lambda(s_1)-\lambda(s_2)\|_{L^2(\Omega)}^2+ \|\gamma^\prime_{1}(s_1)-\gamma^\prime_{1}(s_2)\|_{L^2(\Omega)}^2\\[3pt] &+\|f_w(s_1)-f_w(s_2)\|_{L^2(\Omega)}^2 +\|\gamma_2(s_1)-\gamma_2(s_2)\|_{L^2(\Omega)}^2\\[3pt] &\leq C(D(s_1)-D(s_2), s_1-s_2), \qquad\qquad 0\leq s_1, \, s_2\leq 1. \end{array} \eqno(3.4) $$ A sufficient condition for (3.4) to hold will be described later in this section. The proof of Lemma 3.1 and Theorem 3.2 can be found in [15]. \medskip \noindent {\bf Lemma 3.1.} {\it Let $(p,u,s)$ and $(p_\beta, u_\beta, s_\beta)$ solve} (3.1) {\it and} (3.3), {\it respectively. Then there is a constant $C$ independent of $\beta$ such that} $$ \|p-p_\beta\|_{L^2(\Omega)}+ \|u-u_\beta\|_{L^2(\Omega)} \leq C\bigl\{ \|\lambda(s)-\lambda(s_\beta)\|_{L^2(\Omega)} +\|\gamma_1(s)-\gamma_1(s_\beta)\|_{L^2(\Omega)}\bigr\}. $$ \medskip \noindent {\bf Theorem 3.2.} {\it Assume that $d_\beta\ge d$ and conditions} (2.3) {\it and} (3.4) {\it are satisfied. Let $(p,u,s)$ and $(p_\beta, u_\beta, s_\beta)$ solve} (3.1) {\it and} (3.3), {\it respectively. Then} $$ \begin{array}{l@{}l} &\|p-p_\beta\|_{L^2(\Omega_T)}^2+ \|u-u_\beta\|_{L^2(\Omega_T)}^2 +\|s-s_\beta\|_{L^\infty(J;H^{-1}(\Omega))}^2\\[7pt] &\qquad+\int_J(D_\beta(s)-D_\beta(s_\beta),s-s_\beta)\mbox{d}\tau \leq C(\beta), \end{array} \eqno(3.5) $$ {\it where} $C(\beta)=C\|D(s)-D_\beta(s)\|^\gamma_{L^\infty[0,1]}$. \medskip \noindent {\bf Corollary 3.3.} {\it Under the assumptions of Theorem} 3.2, {\it we have} $$ \begin{array}{l@{}l} &\|D_\beta(s)-D_\beta(s_\beta)\|_{L^2(\Omega_T)}\leq C(\beta)^{1/2},\\[3pt] &\|s-s_\beta\|_{L^{\mu+2}(\Omega_T)}\leq C(\beta)^{1/(\mu+2)}. \end{array} $$ The first result follows from (3.5) and the obvious inequality $$ (D(s_1)-D(s_2))^2\leq \|d\|_{L^\infty[0,1]} (D(s_1)-D(s_2))(s_1-s_2), \qquad 0\leq s_1, \, s_2\leq1, $$ while the second result follows from (2.4) and (3.5). As in \cite{RfsA}, we now consider a specific example of the regularization $d_\beta$ given by $$ d_\beta(s)=\max\{d(s),\beta^\mu\}. \eqno(3.6) $$ Note that $$ \|D(s)-D_\beta(s)\|_{L^\infty[0,1]}\leq C\beta^{\mu+1}, $$ for $\beta$ small enough, so $$ C(\beta)\leq C\beta^{\mu+2}.\eqno(3.7) $$ We end this section with a remark on condition (3.4). Let $\eta$ represent one of the quantities $\lambda$, $f_w$, $\gamma_1^\prime$, and $\gamma_2$. It is clear that if $\eta$ satisfies that $$ |\eta(s_1)-\eta(s_2)|^2\leq C(D(s_1)-D(s_2))(s_1-s_2),\quad \mbox{ a. e. } \, 0\leq s_1, \, s_2 \leq 1, \eqno(3.8) $$ then assumption (3.4) is true for $\eta$. A necessary and sufficient condition for (3.8) to hold is that \cite{RceNM} $$ |\eta_s|^2\leq Cd(s), \qquad \mbox{ a. e. } \, s\in [0,1]. \eqno(3.9) $$ Inequality (3.9) means that $\eta_s$ vanishes with $d$. Below we give the conditions on $\eta$ so that (3.8) or (3.9) holds. The proof of the next proposition can be found in \cite{RchenB} or \cite{RfsA}. \medskip \noindent {\bf Proposition 3.4.} {\it Assume that $\eta\in C^1[0,1]$, $\eta_s(0)=\eta_s(1)=0$, $\eta_s$ is Lipschitz continuous at $0$ and $1$, and assumption {\rm(2.3)} is satisfied. Then there is a constant $C>0$ such that {\rm(3.8)} holds}. \section{Perturbation Analysis} In this section we report a formal application of the perturbation method for the saturation equation $$ {\partial_t s}-\nabla \cdot\left\{d(s)\nabla s\right\}+\nabla \cdot\left\{ f_{w}(s)u\right\}-\nabla \cdot \gamma_2(s)=q_w\,. \eqno(4.1) $$ The perturbation method in \cite{Rmhol} is applied to analyze numerical solutions of this problem. We assume that the numerical method we will develop produces a solution in the asymptotic form $$ s_{\epsilon}\sim s_1+\epsilon s_{2}+\cdots, \eqno(4.2) $$ where $s_{1}$ is the exact solution of (4.1), $s_2$ is a smooth function, and $\epsilon$ is a small constant. Substituting (4.2) into (4.1), we see that $$ \begin{array}{l@{}l} &{\partial_t s_1}+\epsilon {\partial_t s_2} -\nabla \cdot\left\{d(s_1+\epsilon s_2)\nabla(s_1+\epsilon s_2)\right\}\\[3pt] &\qquad+\nabla \cdot\left\{f_{w}(s_1+\epsilon s_2)u\right\}-\nabla \cdot \gamma_2(s_1+\epsilon s_2)+\cdots\sim q_w\,. \end{array} $$ Since $s_1$ satisfies (4.1), we have $$ \begin{array}{l@{}l} &\epsilon {\partial_t s_2}-\nabla \cdot\left\{ d(s_1+\epsilon s_2)\nabla (s_1+\epsilon s_2)-d(s_1)\nabla s_1\right\}\\[3pt] &\quad+\nabla\cdot\left\{\left(f_{w}(s_1+\epsilon s_2) -f_{w}(s_1)\right)u\right\}\\[3pt] &\quad-\nabla\cdot\left\{\gamma_2(s_1+\epsilon s_2)-\gamma_2(s_1)\right\}+\cdots\sim 0. \end{array} $$ That is, $$ \begin{array}{l@{}l} &{\partial_t s_2}-\nabla \cdot\left\{s_2\frac{d(s_1+\epsilon s_2)-d(s_1)}{\epsilon s_2}\nabla s_1\right\}-\nabla \cdot\left\{ d(s_1+\epsilon s_2)\nabla s_2\right\}\\[7pt] &\qquad+\nabla \cdot\left\{s_2\frac{f_w(s_1+\epsilon s_2)-f_{w}(s_1)} {\epsilon s_2}u\right\}-\nabla \cdot\left\{s_2\frac{\gamma_2(s_1+\epsilon s_2)-\gamma_2(s_1)}{\epsilon s_2}\right\}\sim 0. \end{array}\eqno(4.3) $$ Assuming that $d$, $f_w$, $\gamma_2\in C^1[0,1]$ as in Proposition 3.4, it follows from (4.3) as $\epsilon \rightarrow 0$ that $$ \begin{array}{l@{}l} &{\partial_t s_2}-\nabla \cdot \left\{s_2 d_s(s_1)\nabla s_1\right\}-\nabla\cdot\left\{d(s_1)\nabla s_2\right\}\\[3pt] &\quad+\nabla \cdot\left\{s_2 f_{ws}(s_1)u\right\}-\nabla \cdot\left\{s_2\gamma_{2s}(s_1)\right\}\sim 0; \end{array} $$ i.e., $$ {\partial_t s_2}-\nabla \cdot\left\{d(s_1)\nabla s_2\right\} + \nabla \cdot\left\{s_2\left(f_{ws}(s_1)u -\gamma_{2s}(s_1)-d_s(s_1)\nabla s_1\right)\right\}\sim 0\,. \eqno(4.4) $$ By (3.9), note that $f_{ws}$ and $\gamma_{2s}$ vanish with $d$, so (4.4) reduces to the following equation near the degeneracy of the function $d(s)$ with some assumptions for the functions $\gamma_{2}$ and $f_{w}$: $$ {\partial_t s_2}+\nabla \cdot\left\{s_2\left(-d_{s}(s_1)\nabla s_1\right) \right\}\sim d_{s}(s_1)\nabla s_1 \cdot \nabla s_2\,. $$ Hence the above formal analysis indicates that the behavior of errors close to the degeneracy is exponential and $-d_{s}(s_1)\nabla s_1$ determines the gross rate of the errors in time. This is the case as shown in our numerical experiments later. \section{Finite Element Approximation} For notational convenience, we consider the case of $\varphi_1\equiv0$ in the analysis below; otherwise, $\varphi_1$ can be incorporated into the differential equation, or the later mixed finite element method can be handled by introducing the space of Lagrange multipliers \cite{RceSIAMNUMB}. For $d=2$ or $3$, let $$ H(\mbox{div},\Omega)=\{v\in(L^2(\Omega))^d : \nabla\cdot v\in L^2(\Omega)\}, $$ and $$ V=\{v\in H(\mbox{div},\Omega) : v\cdot\nu=0\mbox{ on }\Gamma\},\quad W=\{w\in L^2(\Omega) : w_\Omega=0\}. $$ For $00$ for any $\beta >0$, so $D_\beta$ has a $C^1$ inverse function ${\cal S}_\beta$. Finally, let $P_h$ indicate the $L^2$-projection operator onto $M_h$. As remarked before, the pressure equation is approximated by the mixed finite element method. For each $t\in\bar J$, the mixed finite element solution $(u_h(\cdot,t), p_h(\cdot,t))\in V_h\times W_h$ satisfies $$ \begin{array}{l@{}l@{}l} &(\nabla\cdot u_h, w)=(q,w), &\quad \forall w\in W_h,\\[3pt] &(a(s_h)u_h,v)-(p_h,\nabla\cdot v)=(\gamma_1(s_h), v), &\quad \forall v\in V_h, \end{array}\eqno(5.3) $$ where $s_h$ is determined below. First, for each $t\in J$ we define $\theta_h(\cdot,t)\in M_h$ by $$ \begin{array}{l@{}l} &(\partial_t {\cal S}_\beta(\theta_h),v) +(\nabla\theta_h-f_w\big({\cal S}_\beta(\theta_h)\big)u_h\\[3pt] &\qquad+\gamma_2\big({\cal S}_\beta(\theta_h)\big), \nabla v) +(\varphi_2,v)_{\Gamma}=(q_w,v), \quad\forall v\in M_h, \end{array} \eqno(5.4) $$ with the initial approximation $$ P_h{\cal S}_\beta(\theta_h(0))=P_hs_0. \eqno(5.5) $$ Now, we determine $s_h$ by $s_h={\cal S}_\beta(\theta_h)$, which approximates $s_\beta$. Notice that approximating $D_\beta(s_\beta)$ by $\theta_h$, then $s_\beta$ by ${\cal S}_\beta(\theta_h)$, yields a higher rate of convergence than approximating $s_\beta$ by an element in $M_h$ directly \cite{RceNM, RfsB}. Also, note that if bases are introduced in $V_h$, $W_h$, and $M_h$, equations (5.3)--(5.5) can be written as a nonlinear system of ordinary differential equations for $s_h$ (after substituting (5.3) into (5.4)). With our assumptions on the data, this nonlinear system can be shown to have a unique solution from the fundamental theorem of ordinary differential equations \cite{Rcee}. An error analysis for (5.3)--(5.5) is given in the second paper \cite{RckB}. \subsection{A fully discrete approximation} For each positive integer $N$, let $0=t^0< t^1<\cdots< t^N=T$ be a partition of $J$ into subintervals $J^n=(t^{n-1},t^n]$ with length $\Delta t^n=t^n-t^{n-1}$, $1\leq n\leq N$. Also, set $v^n=v(\cdot,t^n)$. Finally, we indicate the time difference operator by $$ \partial v^n=\frac{v^n-v^{n-1}}{\Delta t^n}, \quad 1\leq n\leq N. $$ Now, the fully discrete approximation is given as follows. For any $1\leq n\leq N$, find $(u^n_h, p^n_h)\in V_h\times W_h$ such that $$ \begin{array}{l@{}l@{}l} &(\nabla\cdot u^n_h, w)=(q^n,w) &\quad \forall w\in W_h,\\[3pt] &(a(s^n_h)u^n_h,v)-(p^n_h,\nabla\cdot v)=(\gamma_1(s^n_h), v) &\quad \forall v\in V_h, \end{array}\eqno(5.6) $$ where $s_h={\cal S}_\beta(\theta_h)$ and for each $n$, $\theta^n_h\in M_h$ satisfies $$ \begin{array}{l@{}l} &(\partial {\cal S}_\beta(\theta^n_h),v) +(\nabla\theta^n_h-f_w\big({\cal S}_\beta(\theta^n_h)\big)u^n_h\\[3pt] &\qquad+\gamma_2\big({\cal S}_\beta(\theta^n_h)\big), \nabla v) +(\varphi^n_2,v)_{\Gamma}=(q^n_w,v), \quad\forall v\in M_h, \end{array} \eqno(5.7) $$ with the initial approximation given as in (5.5). \begin{center} \begin{tabular}{|c|c|c|c|c|} \hline \vbox to 11pt{} $1/\Delta x$ & $p-p_h$ & rate for $p$ & $u-u_h$ & rate for $u$\\[3pt] \hline 10 & 0.0585 & - & 0.0280 & - \\ \hline 20 & 0.0277 & 1.07 & 0.0073 & 1.94 \\ \hline 40 & 0.0135 & 1.04 & 0.0019 & 1.97 \\ \hline 80 & 0.0066 & 1.02 & 0.0005 & 1.98 \\ \hline \end{tabular} \\[6pt] {Table~1a. The error estimates for $p$ and $u$.} \end{center} The remark on existence and uniqueness of (5.6) and (5.7) can be made as above \cite{RceSIAMNUMB}. Also, an error analysis for this approximation is presented in \cite{RckB}. \begin{center} \begin{tabular}{|c|c|c|c|c|} \hline \vbox to 11pt{} $1/\Delta x$ & $s-s_h$ for & rate for & $s-s_h$ for & rate for \\ & $\beta=0$ & $\beta=0$ & $\beta=\beta_0$ & $\beta=\beta_0$ \\[3pt] \hline 10 & 0.0858 & - & 0.5443 & - \\ \hline 20 & 0.0452 & 0.91 & 0.2931 & 0.90 \\ \hline 40 & 0.0272 & 0.74 & 0.1236 & 1.24 \\ \hline 80 & 0.0162 & 0.75 & 0.0530 & 1.22 \\ \hline \end{tabular} \\[6pt] {Table~1b. The error estimates for $s$ in the first example.} \end{center} \medskip \section{Numerical Results} The numerical experiments are presented to show convergence of our approximation scheme, to demonstrate the qualitative behavior of error estimates, to compare the present regularization technique with the un-regularized version used in \cite{RceNM}, and to complement the perturbation analysis in the fourth section. Toward that end, we consider the pressure-saturation system of the form $$ \begin{array}{l@{}l@{}l} &-\nabla \cdot\left(\lambda(s)\nabla p\right) =q &\quad\mbox{ in } \Omega_T,\\[3pt] &\partial_t s -\nabla\cdot\big\{d(s)\nabla s+\lambda_{w}(s)\nabla p\big\}=q_{w} &\quad\mbox{ in } \Omega_T. \end{array}\eqno(6.1) $$ For simplicity we focus on the Dirichlet boundary conditions $$ p=p_D, \quad s=s_D \quad \mbox{ on } \Gamma\times J. $$ The initial condition is $$ s(x,0)=s_0(x),\quad x\in\Omega. $$ %\begin{figure}[ht] %\centering %\scalebox{.5}{\includegraphics{fig1.ps}} %\caption{The error $s-s_h$ for $\beta=\beta_0$ in the first example.} % \end{figure} \setcounter{figure}{1} \medskip \noindent{\bf Example 1.} The domain $\Omega$ is taken to be the unit square and other data are chosen as follows: $$ \lambda(s)=1, \quad d(s)=s(1-s), \quad \lambda_{w}(s)=1. $$ In the first example the exact solution of system (6.1) is of the form $$ p(x,t)=\sin(\pi x)\sin(\pi y),\quad s(x,t)=t\sin(\pi x)\sin(\pi y). $$ The boundary and initial data coincide with the exact solution on the boundary and at the initial time. Also, the functions in the right-hand side of system (6.1) result from the exact solution. The numerical experiments reported here are mainly to show the qualitative behavior of the finite element approximation for the degenerate saturation equation. This is why we consider the case where $\lambda(s)=1$, so the pressure equation is decoupled from it. We also carried out experiments with $\lambda(s)$ depending on $s$, which are not reported here and have similar results to those illustrated here. Furthermore, note that the data satisfy assumptions (2.3) with $\mu_1= \mu_2=\mu=1$ and (3.4). Uniform partitions of $\Omega$ into triangles are used, with $\Delta x=\Delta y$ as the lengths in the $x$ and $y$ directions, and the lowest-order Raviart-Thomas mixed finite element on the triangles \cite{Rrt} are exploited. The mixed method is used for the pressure equation and the standard finite element method with the backward Euler scheme for the time differentiation term is utilized for the saturation equation, as in (5.6) and (5.7). The time step is taken to be proportional to the space step. The mixed finite element method is implemented as in \cite{RchenEW}. Error estimates and convergence rates in the $L^\infty$-norm for the approximations to the pressure and velocity $u=-\lambda(s)\nabla p$ are presented in Table~1a. From this table we see that the mixed method is first-order accurate for the pressure and second-order accurate for the velocity. That is, a superconvergence rate occurs for the velocity. The error estimates in the $L^\infty$-norm for the saturation at $t=1$ are described in Table~1b. We consider the cases where the regularization parameter $\beta$ is taken to be zero or $\beta_0\equiv Ch^{\mu_0}$ with $\mu_0$ given by $$ \mu_0=\frac{4}{3\mu+2}, $$ which is chosen according to the error analysis performed in the paper \cite{RckB}. The convergence rates in the case of $\beta=\beta_0$ are better than those in the case of $\beta=0$. Also, the convergence rates in the former case coincide with the theoretical results obtained in \cite{RckB}. Finally, the error with $\Delta x=1/80$ for $s_h$ with $\beta={\beta_0}$ is shown as Figure 1 in a separate file. %Figure~\ref{fig1}. It can be seen that maximum errors occur in the center where $d(s)$ is zero and near the boundary of the domain where $d(s)$ is close to zero. This complements the perturbation analysis in the forth section. \medskip \begin{figure}[ht] \centering \scalebox{.5}{\includegraphics{fig2.ps}} \caption{ The error $s-s_h$ for $\beta=\beta_0$ in the second example.} \end{figure} \noindent{\bf Example 2.} In this example, we simulate a relatively simple two-phase flow problem. The data are given as follows: $$ \begin{array}{l@{}l@{}l@{}l@{}l} &k_{rw}=s, &\quad k_{ro}=1-s, &\quad \mu _{w}=0.5\mbox{ cp}, &\quad \mu _{o}=2\mbox{ cp},\\[3pt] &\rho _{w}=1 \mbox{ g/cm}^{3}, &\quad \rho_{o}=0.7\mbox{ g/cm}^{3}, &\quad \phi =0.2, &\quad k=0.05 \mbox{ darcy}. \end{array} $$ Moreover, the function $p_{c}$ is given by $$ p_{c}(s)=1-s\,. $$ With these, the system of equations (3.1a) now reduces to $$ \begin{array}{l@{}l} &\nabla \cdot u=q\,, \\[3pt] &u=-a(s)(\nabla p-G(s))\,, \\[3pt] &\phi {\partial_t s}-\nabla \cdot \left\{d(s)\nabla s-f_{w}(s)u-b(s)\right\}=q_{w}\,, \end{array} $$ where $$ \begin{array}{l@{}l@{}l} &a(s) =0.0125(3s+1), &\quad G(s)=(0.7+3.3s)\tilde g, \\[3pt] &d(s) =0.1s(1-s)/(3s+1), &\quad f_{w}(s)=4s/(3s+1),\\[3pt] &b(s)=0.03s(1-s)/(3s+1)\tilde g\,, &\quad \end{array} $$ where $\tilde g$ is the gravity vector. \begin{center} \begin{tabular}{|c|c|c|c|c|} \hline \vbox to 11pt{} $1/\Delta x$ & $p-p_h$ for & rate for & $p-p_h$ for & rate for \\ & $\beta=0$ & $\beta=0$ & $\beta=\beta_0$ & $\beta=\beta_0$ \\[3pt] \hline 10 & 0.075626 & - & 0.069228 & - \\ \hline 20 & 0.036445 & 1.0532 & 0.031686 & 1.1275 \\ \hline 40 & 0.017600 & 1.0501 & 0.014920 & 1.0866 \\ \hline 80 & 0.008615 & 1.0307 & 0.007110 & 1.0693 \\ \hline \end{tabular} \\[6pt] {Table~2a. The error estimates for $p$ in the second example.} \end{center} \bigskip \begin{center} \begin{tabular}{|c|c|c|c|c|} \hline \vbox to 11pt{} $1/\Delta x$ & $s-s_h$ for & rate for & $s-s_h$ for & rate for \\ & $\beta=0$ & $\beta=0$ & $\beta=\beta_0$ & $\beta=\beta_0$ \\[3pt] \hline 10 & 0.2905 & - & 0.2423 & - \\ \hline 20 & 0.1153 & 1.3331 & 0.0969 & 1.3222 \\ \hline 40 & 0.0528 & 1.1268 & 0.0430 & 1.1722 \\ \hline 80 & 0.0253 & 1.0614 & 0.0203 & 1.0829 \\ \hline \end{tabular} \\[6pt] {Table~2b. The error estimates for $s$ in the second example.} \end{center} \medskip The numerical results corresponding to those in Example 1 are given in Table~2 and Figure~2. Similar observations can be made here. \medskip \noindent{\bf Example 3.} In the final example, we test a more physically adequate set of data. We simulate a two-phase flow problem \cite{Rbear}. The function $p_{c}(s)$ is given by $$ p_{c}(s) =(1-s)\left\{\gamma (s^{-1}-1)+\theta\right\}\,, $$ where $$ \gamma=20,000\mbox{ dynes/cm}^{2},\quad \theta =100\mbox{ dynes/cm}^{2}\,. $$ The relative permeabilities are defined by $$ k_{ro} =\left\{ \begin{array}{l@{}l@{}l} &0 &\qquad \mbox{ if } s>s_{o}\,,\\[3pt] &s_{o}^{-2}(s_{o}-s)^{2} &\qquad\mbox{ if }0\leq s\leq s_{o}\,. \end{array}\right. $$ and $$ k_{rw} =\left\{ \begin{array}{l@{}l@{}l} &(s-s_{rw})^{2}(1-s_{rw})^{-2} &\qquad\mbox{ if }s\geq s_{rw}\,,\\[3pt] &0 &\qquad \mbox{ if } \;0\leq s