\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