\documentclass[reqno]{amsart}
\usepackage{graphicx}


\AtBeginDocument{{\noindent\small 
{\em Electronic Journal of Differential Equations}, 
Vol. 2003(2003), No. 25, pp. 1--10.\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{9mm}}

\begin{document}

\title[\hfilneg EJDE--2003/25\hfil Stability of the interface between two
immiscible fluids]
{Stability of the interface between two immiscible fluids in porous media}

\author[C. I. Calugaru, D.-G Calugaru, J.-M. Crolet, \& M. Panfilov
\hfil EJDE--2003/25\hfilneg] 
{Cerasela I. Calugaru, Dan-Gabriel Calugaru,\\
 Jean-Marie Crolet, \& Michel Panfilov}

\address{Cerasela I. Calugaru \hfill\break
Equipe de Calcul Scientifique, Universit\' e de Franche-Comt\' e\\
16 Route de Gray, 25030 Besan\c con Cedex France}
\email{calugaru@math.univ-fcomte.fr}

\address{Dan-Gabriel Calugaru \hfill\break
Universit\'{e} Claude Bernard Lyon 1, ISTIL, MCS/CDCSP \\
15, Boulevard Latarjet, 69622 Villeurbanne Cedex, France}
\email{calugaru@cdcsp.univ-lyon1.fr}

\address{Jean-Marie Crolet \hfill\break
Equipe de Calcul Scientifique, Universit\' e de Franche-Comt\' e\\
16 Route de Gray, 25030 Besan\c con Cedex France}
\email{jmcrolet@univ-fcomte.fr}

\address{Michel Panfilov \hfill\break
LAEGO - ENS de Geologie  - INP de Lorraine \\
Rue M. Roubault, BP 40, F - 54501 Vandoeuvre-l\`{e}s-Nancy France}
\email{michel.panfilov@ensg.inpl-nancy.fr}

\thanks{ {\em Mathematics Subject Classifications:} 76B15, 35R35, 76T05.
\hfil\break\indent
{\em Keywords :} porous media, two-phase flow, interface, instability.
\hfil\break\indent
\copyright 2003 Southwest Texas State University. \hfil\break\indent
Submitted January 20, 2003. Published March 10, 2003.}
\date{}


\begin{abstract}
 A generalized model has been recently proposed in \cite{Pan}
 to describe deformations of the mobile interface separating two
 immiscible and compressible fluids in a deformable porous medium.
 This paper deals with a few applications of this model in
 realistic situations where it can be supposed that gravity
 perturbations are propagating much slower than elastic
 perturbations. Among these applications, one can include the
 classical well-known case of groundwater flow with free surface,
 but also more complex phenomena, as gravitational instability with
 finger growth.
\end{abstract}

\maketitle

\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}

\section{Introduction}

Evolution of the interface between two immiscible fluids is one of
the more complicated problem in mathematical analysis. It is
usually described by a system of partial differential equations
which are true on either side of the interface. The closure of the
system is assured by some dynamic and kinetic conditions on the
interface. The principal difficulty of this kind of problem is to
transform such a system into an explicit closed differential
equation which governs the interface movement.

For instance, if we deal with two superposed fluids and if
$h(x_1,x_2;t)$ denotes the thickness of the lower fluid, then the
problem is to deduce the closed differential equation for the
function $h(x_1,x_2;t)$ (here $x_1, x_2$ are the coordinates in
horizontal plane, $t$ is the time variable and $z=h(x_1,x_2;t)$ is
the equation of the interface considered as a mobile surface).

In natural gas-oil reservoirs, such a model is basic to describe
behavior of the stratified pair "gas-oil" or "oil-water". Indeed,
before any recovery, the natural non perturbed system is generally
a porous medium composed of three superposed layers (gas, oil and
water, in downward direction). The pumping of oil leads to an
ascent of water from below or/and to an overlapping of the oil
layer by broken gas which is coming from above. Then, for an
optimized recovery, it is absolutely necessary to have a good
modelling of the evolution of both interfaces (gas-oil and
oil-water).

Explicit differential equations describing the evolution of such
an interface have been already obtained for the case when the
viscosity of one of the two fluids can be neglected, as for
instance in the case of water-air flow \cite{Whit} in the
framework of shallow water theory, or in the case of the
groundwater flow in unconfined aquifers \cite{BarEnRy},
\cite{Bear}, \cite{Dagan}, \cite{LW}, \cite{Parl}. Whitham
\cite{Whit} has considered the case where the viscous forces can
be neglected for both fluids, but this situation is rarely
encountered in porous media. When both fluids are viscous, some
explicit equations were deduced in \cite{PolKoch}.

In all previous studies, some hypotheses (as vanishing of vertical
flow velocity (gravity equilibrium condition), constant horizontal
flow velocity along $z$, or steady-state flow) are imposed. Such
assumptions seem to be reasonable if the interface deformation is
rather small, if the boundary is free and if the transition
phenomena are not taken into account, but become excessively
strong in the case of two-liquids flow.

A new model has been recently proposed in \cite{Pan} for two
viscous fluids in porous media. It generalizes previous models by
removing the classical condition of gravity equilibrium and by
considering a non-stationary flow (the compressibility of the
fluids and of the porous medium is not neglected). Instead of
previous hypotheses, a more general assumption about linear
behavior of the vertical velocity along vertical coordinate is
used. This generalized model has been already applied to describe
the flow of rather compressible fluids, when an elastic
perturbation is propagating in the domain much slower than a
gravity perturbation.

However, the previous model is able to describe many other
situations, corresponding, for instance, to the case when an
elastic perturbation is propagating much faster than a gravity
perturbation. This paper deals with a few applications of the
model in this last case. The outline is as follows: in Section 1,
we give a brief description of the physical problem and recall the
general model deduced in \cite{Pan}. For the particular case
investigated here, i.e. fast elastic perturbations and slow
gravity perturbations, a reduced model is obtained in Section 2.
It can be converted into one nonlinear equation of third order. In
particular, a generalization of Boussinesq equation is obtained
for groundwater flow with free surface. Stable flow is taken into
account when the upper fluid is more light. If this condition is
not true, the evolution of instability is described.

\section{General model}

We consider an orthogonal coordinate system ($x_1, x_2, z$), where
$z$ is the vertical coordinate and ($x_1, x_2$) are coordinates of
the horizontal plane. In this system, a horizontal porous medium
with constant height $H$ is considered. Its inferior and superior
limits are some impermeable media. The voids of the porous medium
are filled with two immiscible fluids which are separated one from
other by an interface. Let the value $h(x_1, x_2; t)$ denotes the
height of the interface with respect to the bottom of the porous
stratum (Figure \ref{fi0}).

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig1.eps}
\caption{Schematization of the physical situation} \label{fi0}
\end{center}
\end{figure}

The initial position of the interface is supposed to be known (as
for instance, a horizontal surface) and denoted by $h_0$. At this
time, we suppose there is no flow in the domain. The initial state
is perturbed by an external factor (as pumping by a well or
seismic waves). The objective is to find the evolution of the
interface, i.e. of the function $h(x_1,x_2;t)$, generated by such
a perturbation. It is supposed that the interface does not cross
the top and bottom boundaries of the domain and it remains always
non perturbed at the lateral extremities.

\subsection{Physical parameters}

The porous medium is assumed to be homogeneous (its constant
porosity is denoted by $m$), but anisotropic, the permeability
tensor $K{\equiv}\{ K_{ij}\}_{i,j=1}^3$ being diagonal with
$K_{x}{\equiv}K_{11}=K_{22}$, $K_z{\equiv}K_{33}$.

Let the indexes $i = I,II$ correspond to the lower and the upper
fluid. For each fluid, $\rho^i$ denotes the density, $\mu^i$
denotes the viscosity and $\beta^i_*$ is a measure of fluid/medium
compressibility.

We need also to introduce the following parameters: $L$ is the
horizontal scale of the domain, \ ${\mathcal P}_0$ and \
${\mathcal P(x_1,x_2;t)}$ are respectively the common pressure of
the two fluids on the interface at initial state and for an
arbitrary instant time, $g$ is the gravity acceleration.

\subsection{Closed system in dimensionless formulation}

The time variable and the space variables are replaced by some
dimensionless variables:
\[
\tau{\equiv}{t}/{t_*}, \quad y_1{\equiv}x_1/L, \quad
y_2{\equiv}x_2/L
\]
%
and the unknowns are the thickness of the two fluids and the
common pressure on the interface, all of them in dimensionless
form:
\[
\varphi{\equiv}h/{h_0},\quad \psi{\equiv}(H{-}h)/(H{-}h_0), \quad
\xi{\equiv}{\mathcal P}/{\mathcal P}_0
\]
Then the dimensionless system deduced in \cite{Pan} is:
\begin{subequations}\label{eq2.8} 
\begin{gather} 
\left(\tau_*{\partial \over \partial \tau}{-}\Delta_{yy} \right)
\left( \varphi^2{+}\frac{\lambda_1^2\omega\tau_*}{3\varepsilon_z}
      \varphi^2{\partial \varphi\over \partial \tau} {+} \lambda_2\xi \right) =
  -\omega\tau_*{\partial \varphi\over \partial \tau},
                                 \label{eq2.8a}     \\ 
 \left( \frac{\tau_*\overline{\beta}}{\overline{\mu}}
         {\partial \over \partial \tau}{-}\Delta_{yy} \right) \left( {-}\psi^2 {+}
      \frac{\lambda_1^2\overline{\rho}\omega\tau_*}
           {3\varepsilon_z\overline{\mu}\lambda_0} \psi^2{\partial \psi\over \partial \tau} {+}
  \lambda_2\lambda_0\overline{\rho} \xi \right)  =
       -\frac{\overline{\rho}\lambda_0\omega\tau_*}{\overline{\mu}}
        {\partial \psi \over \partial \tau}      \label{eq2.8b},     \\ 
\psi = -\lambda_0\varphi{+}\lambda_0{+}1 \label{eq2.8c}
\end{gather}
\end{subequations} 
with the following notation
\begin{gather*}
\lambda_0 \equiv \frac{h_0}{H-h_0} \,,\quad
\lambda_1 \equiv \frac{h_0}{L} \,,\quad
\lambda_2 \equiv \frac{2{\mathcal P}_0}{\rho^I gh_0} \\
 \overline{\beta}\equiv\frac{\beta^I_*}{\beta^{II}_*}\,,\quad 
 \overline{\mu} \equiv \frac{\mu^I}{\mu^{II}}\,,\quad 
 \overline{\rho} \equiv \frac{\rho^I}{\rho^{II}}\,,\quad 
 \varepsilon_z \equiv \frac{K_z}{K_x} \\
t_{el}= \frac{\mu^I L^2}{K_x \beta^I_*} \,,\quad
t_{gr}= \frac{2\mu^I m L^2}{K_x \rho^Igh_0} \,,\quad 
\omega\equiv \frac{2m \beta^I_*}{\rho^I gh_0} =
\frac{t_{gr}}{t_{el}}\,,\quad \tau_*\equiv\frac{t_{el}}{t_*}\,,
\end{gather*}
where $\Delta_{yy}$ denotes the Laplace's operator with respect to the
variable $y=(y_1,y_2)$.

Two characteristic times have been introduced: $t_{el}$ signifies
the time of propagation of an elastic perturbation within the
scale $L$ and $t_{gr}$ is the time required for the complete
emptying of the medium occupied by a fluid (say the fluid $I$) if
the gravity drop is the only phenomenon taken into account.

In the general case, the scale of the time ($t_*$) used to get the
dimensionless time variable ($\tau$) can be arbitrary chosen.
However, for two particular but important cases, it can be defined
according to the relationship between previous characteristic
times, as follows:
\begin{equation}\label{eq2.7}
t_*=   \begin{cases}
   t_{el}, & \mbox{when }  t_{gr} \ll t_{el},\mbox{ or } \omega \ll 1 \\
   t_{gr}, & \mbox{when }  t_{el} \ll t_{gr},\mbox{ or } \omega \gg 1
   \end{cases}
\end{equation}
The first case, i.e. $t_{gr} \ll t_{el}$ has been already
investigated in \cite{Pan}.

\section{Gravity dominated flow. Fast elastic perturbations}
\label{sec3}

In this section, we deal with the case where the time of elastic
wave propagation is very small with respect to the gravity time,
i.e. $t_{el} \ll t_{gr}$. Therefore, we have:
$\omega \gg 1$
In this case, compressibility of fluid/medium does not play any
role and then, the non-stationarity is involved only by gravity
phenomena.

According to (\ref{eq2.7}), the scale of the time $t_*$ is chosen
as equal to $t_{gr}$. Then:
$$
\omega\tau_* = 1,\quad \tau_* = \omega^{-1} \ll 1$$


\subsection{General nonlinear equation for gravitational waves}

Taking into account previous relations, the general model
(\ref{eq2.8a})-(\ref{eq2.8c}) becomes:
%
\begin{subequations} \label{eq3.1} 
\begin{gather} 
 \Delta_{yy}\left(\varphi^2{+}\frac{\lambda_1^2}{3\varepsilon_z}
      \varphi^2{\partial \varphi\over \partial \tau} {+} \lambda_2\xi \right) =
  {\partial \varphi\over \partial \tau},
                                 \label{eq3.1a}     \\ 
 \Delta_{yy}\left( {-}\psi^2 {+}
      \frac{\lambda_1^2\overline{\rho}} {3\varepsilon_z\overline{\mu}\lambda_0}
      \psi^2{\partial \psi\over \partial \tau} {+}
  \lambda_2\lambda_0\overline{\rho}\xi \right)  =
      \frac{\overline{\rho}\lambda_0}{\overline{\mu}} {\partial \psi\over \partial \tau}
                                \label{eq3.1b},     \\ 
 \psi = -\lambda_0\varphi{+}\lambda_0{+}1 \label{eq3.1c}
\end{gather}
\end{subequations}
This system can be converted into one equation with respect to the
function $\varphi$, after excluding the functions $\xi$ and
$\psi$:
%
\begin{equation}
\Delta_{yy}\left( \varphi^2 {-} 2\alpha \varphi {+}
      \gamma_1 \left[ \varphi^2{+}\gamma_2
             \big( 1{+}\lambda_0{-}\lambda_0\varphi\big)^2\right] {\partial \varphi\over \partial \tau}
      \right)
       = \beta{\partial \varphi \over \partial \tau}        \label{eq3.2}
\end{equation}
%
where
%
\begin{equation}
\alpha{\equiv}\frac{\big( 1{+}\lambda_0 \big) }
                   {\lambda_0{+}\overline{\rho}},\quad
\beta{\equiv} \frac{\overline{\rho}\big( \lambda_0 {+}\overline{\mu}\big) }
    {\overline{\mu}\big( \overline{\rho}{+}\lambda_0\big) },\quad
\gamma_1{\equiv}\frac{ \lambda_1^2\overline{\rho} }
                     { 3\varepsilon_z(\lambda_0{+}\overline{\rho}) },\quad
\gamma_2{\equiv}\frac {1} { \lambda_0\overline{\mu} }
                  \label{eq3.2'}
\end{equation}
Due to the presence of a mixed derivative of third order, it is
difficult to determine the type of this equation, but in some
particular cases it can be analyzed.
%-----------------------------------------

\subsection{One-fluid flow with free boundary}

The following generalization of the classical nonlinear Boussinesq
equation can be deduced from equation (\ref{eq3.2}) for
single-phase flow of the lower liquid with free boundary. In fact,
assuming that the upper fluid has no density and viscosity, i.e.
$\overline{\rho}{\to}\infty$ and $\overline{\mu}{\to}\infty$, we
get: \ $ \alpha = 0$, \ $\beta = 1$, \
$\gamma_1 = \lambda_1^2{/}\big( 3\varepsilon_z \big)$,
$\gamma_2 = 0 $,\ and then
\begin{equation}
\Delta_{yy}\left( \varphi^2 {+}
      \gamma_1 \varphi^2 {\partial \varphi\over \partial \tau} \right)  = {\partial \varphi\over \partial \tau}   \label{eq3.3}
\end{equation}

When the porous layer is very thin, i.e. its horizontal dimension
is much larger than the vertical dimension, we have
\begin{equation}
\lambda_1 \ll 1\label{eq3.4}
\end{equation}
and then, the Boussinesq equation is obtained:
\begin{equation}
\Delta_{yy} \varphi^2  = {\partial \varphi\over \partial \tau}
\label{eq3.5}
\end{equation}
which is the classical object of the theory of shallow water flow
in porous layers \cite{BarEnRy,Bear,LW,PolKoch}.
%-------------------------------------------------

\subsection{Two fluids, thin layer. Gravitational instability }

Let us examine two fluids in a thin porous layer, such as the
assumption (\ref{eq3.4}) holds. Then (\ref{eq3.2}) becomes
\begin{equation}
\Delta_{yy}\left( \varphi^2 {-} 2\alpha \varphi \right)
 = \beta{\partial \varphi\over \partial \tau}, \quad
\mbox{or}\quad {\partial \over
\partial y_i}\left( \big(
\varphi{-}\alpha\big){\partial \varphi\over \partial y_i}\right)
                      = \frac{\beta}{2}{\partial \varphi\over \partial \tau}
                     \label{eq3.6}
\end{equation}
where parameters $\alpha$ and $\beta$ are defined according to
(\ref{eq3.2'}), and the summation is made on the index $i$.

This equation is parabolic when $\varphi>\alpha$. Otherwise,
when $\varphi{<}\alpha$, this is a nonlinear anti-diffusion
equation, which can describe the evolution of instability.

To check the nature of these phenomena, let us examine stability
of the unperturbed state of the interface, basing on the 1D
self-similar solutions of (\ref{eq3.6}). In fact, in 1D case,
(\ref{eq3.6}) has solutions in the form of $\varphi = \varphi
(y/\sqrt{\tau})$, which correspond to a boundary-value problem of
the interface evolution in an infinite horizontal domain. It is
easy to show that self-similar solutions satisfy the following
ordinary differential equation
%
\begin{equation}
\frac{d^2}{d \xi^2}\left( \varphi^2{-}2\alpha\varphi \right)
 = -\frac{1}{2}\beta \xi \frac{d\varphi}{d \xi},\quad
\xi{\equiv}\frac{y}{\sqrt{\tau}}
\end{equation}

Let us examine small perturbation of the unperturbed state, i.e.,
$\varphi = 1{+}\varepsilon\zeta$, where $\varepsilon$ is assumed
to be small. Then, using the perturbation method, we get:
$$
\zeta '' = -\frac{\beta\xi}{4-4\alpha }\zeta
^{\prime}, \ \ \mbox{or} \ \
\zeta(\xi) = C_1{+}C_2\int\limits_{\infty}^\xi
      e^{ \frac{\beta u^2}{8(\alpha{-}1)} }du $$
where the prime denotes derivation with respect to $\xi$, \ and
$C_1$, $C_2$ are some integration constants.

The evolution of perturbations is unlimited in time, when
%
\begin{equation}
\label{eqinstab}
 \alpha{\equiv}\frac{\big( 1{+}\lambda_0 \big) }
                   {\lambda_0{+}\overline{\rho}}>1 
\end{equation}
%
or $\overline{\rho}{<}1$.
Then, when the upper liquid is heavier than the lower liquid, 
(\ref{eq3.6}) can describe the classical phenomenon of
gravitational instability.
%-----------------------------------------------------------

\subsection{Stable flow}

In the case when the condition (\ref{eqinstab}) is not satisfied,
i.e. the upper fluid is more light, the flow is stable. In fact,
since $\alpha{<}1$, then the "diffusion coefficient"
$\varphi{-}\alpha$ is positive in the vicinity of the initial
state ($\varphi{\equiv}1$), and then (\ref{eq3.6}) is
parabolic.

Figure \ref{fi1} illustrates the evolution of an initial
perturbation defined as a single loop of sinus. The initial wave
tends to be dissipated, with finite rate of the wave propagation,
as usual for nonlinear parabolic equations.
%
\begin{figure}[htb]
\begin{center}
\includegraphics[height=5cm,width=10cm]{fig2.eps}
\caption{Evolution of the interface perturbation in stable case;
$\alpha = 0.888888...$ \label{fi1}}
\end{center}
\end{figure}

More precisely speaking, the following initial problem was
considered: 
\begin{gather*} 
\beta{\partial \varphi\over \partial \tau}
={\partial \over \partial y}\left( 2\big(\varphi-\alpha\big){\partial \varphi \over
\partial y}\right), \quad
    y{\in}(0,1),\; \tau>0 \\
 \varphi\vert_{\tau=0}=\begin{cases}
  1, & y{\notin}J \\
  1{+}a\sin{(by+c)}, & y{\in}J
  \end{cases} \\
\varphi\vert_{y=0; 1} = 1          
\end{gather*}
where $J{\subset}(0,1)$ is a small support located far from the
boundaries; $a,b,c$ are some constant parameters.

The time evolution of $\varphi$ has been obtained by numerical
simulation. For the time discretization, a semi-implicit scheme
has been employed. The discrete time instants are denoted by
$\tau_n=n \Delta \tau$, where $\tau_0=0$ and $\Delta \tau$ is the
time step. If the function $\varphi^n(y)\equiv\varphi(y;\tau_n)$
is already computed at time step $\tau_n$, then one computes the
corresponding values at the time instant $\tau_{n+1}$ by solving
the following semi-discretized problem
%
\begin{equation}
\label{schema} \beta \frac{\varphi^{n+1}-\varphi^{n}}{\Delta
\tau}={\partial \over \partial y}\left(
2\big(\varphi^n-\alpha\big){\partial {\varphi^{n+1}}\over \partial
y}\right), \quad
    y{\in}(0,1),\; \forall n\geq0
\end{equation}
For the space discretization, a conventional scheme of order 2,
obtained by a finite difference method, has been used.

The numerical values used to obtain the evolution shown in Figure
\ref{fi1} are
\[
\alpha = \frac{8}{9} \,,\quad \beta = \frac{55}{9} \,,\quad
a=\frac{1}{10}\,,\quad b=10\pi \,,\quad c=-4\pi
\]
These values of $\alpha$ and $\beta$ are obtained for the case of
the pair water-oil, i.e. for $\overline{\rho}=5/4$ ($\rho^I=1000 \
kg/m^3,\ \rho^{II}=800 \ kg/m^3$), $\overline{\mu}=1/10$
($\mu^I=0.001 \ Pa \cdot s, \ \mu^{II}=0.01 \ Pa \cdot s$) and by
considering that the initial thickness is the same for the two
fluids ($\lambda_0=1$). The consecutive time instants described in
Figure \ref{fi1} correspond to the value of the characteristic
time $t_* = t_{gr}= 1$ hour.

\subsection{Evolution of the instability}

If condition (\ref{eqinstab}) is true, then the function
$2(\varphi{-}\alpha )$ is negative in the vicinity of the initial
state $\varphi = 1$, and Eq. (\ref{eq3.6}) is anti-parabolic.
However, if the function $\varphi$ becomes larger than $\alpha$,
then this equation becomes parabolic, as shown in Fig.~\ref{fi2}.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig3.eps}
\caption{Variation of the parabolicity direction for
(\ref{eq3.6}) when $\alpha>1$ \label{fi2}}
\end{center}
\end{figure}

Thus, the evolution of a perturbation follows three basic stages:
\begin{itemize}
\item[i)] the linear unstable stage, when $\varphi$ is in the vicinity of
the initial state ($\varphi \sim 1$); anti-diffusion equation
governs the evolution;

\item[ii)] the nonlinear unstable stage, when $\varphi{\to}\alpha$; the
spatial derivative $\partial \varphi/\partial y$ becomes
unlimited;

\item[iii)] the final stable stage, when $\varphi$ becomes greater than
$\alpha$; diffusion equation describes this stage.
\end{itemize}
Equations with variation of the parabolicity direction were
considered in \cite{Plot}, as the model of hysteresis phenomena in
theory of phase transitions. However, this theory cannot be
applied to (\ref{eq3.6}), since the structure of nonlinearity
is very different. On the other part, some regularization methods
for equations similar to (\ref{eq3.6}) were examined there.
According to \cite{Plot}, the best regularization is obtained by
adding the term $\Delta_{yy}{{\partial_{\tau}}\varphi}$. Then, the
full equation (\ref{eq3.2}) is expected to regularize unstable
solutions.
%------------

\subsubsection*{The linear unstable stage}
Let us examine, in 1D framework, a small perturbation of the
initial plate interface, $\varphi = 1{+}\varepsilon f(t,y)$, where
$\varepsilon$ is the small amplitude of the perturbation, $f(y,t)$
is a function equal to zero everywhere, excluding a small support
$J = \{ y:\ y{\in}[y_0,y_1]\}$. Let the function $f$ behaves as
$f\sim\varepsilon (y{-}y_0)^2$, \ when $y{\to}y_0$, in such a way
that the perturbation is smooth in the vicinity of the point
$y_0$. Since the equation (\ref{eq3.6}) can be expressed as:
$$
2\big( \varphi{-}\alpha\big) {\partial^2 \varphi\over \partial
y^2} + \left({\partial {\varphi}\over \partial y}\right)^2
                      = \beta{\partial \varphi\over \partial \tau} $$
then, the evolution of such a perturbation may be easily analyzed.
Consequent time phases are shown in Fig. \ref{fi3}.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig4.eps}
\caption{Evolution of a smooth monotonous perturbation of the
interface in the unstable case
                 \label{fi3}}
\end{center}
\end{figure}

At initial time, for the point $A = (\varphi,y) = (1,y_0)$ it is
true:
$$
\partial\varphi /\partial y = 0,\quad \partial^2\varphi/\partial y^2>0, \quad
\mbox{then}\ \ (\varphi{-}\alpha)\partial^2\varphi/\partial
y^2{<}0, \ \rightarrow\ \ {\partial_{\tau}}\varphi{<}0 $$ Then, in
Figure \ref{fi3}, the plot I becomes the plot II. The point $A$ is
lowered downwards and then, the extreme point with zeroth spatial
derivative is displacing to the left (to the point $B$). The next
step of the evolution leads to the plot III, etc.

Thus, collecting all the results obtained about the unstable evolution, it
may be concluded that a small localized smooth and monotonous perturbation of
the interface leads to arising of three types of phenomena:
\begin{itemize}
\item[a)] lateral propagation of the perturbation along the interface;

\item[b)] fast development of the spatial and time oscillations of the interface;

\item[c)] fast grow of the amplitude of the oscillations.
\end{itemize}

Figure \ref{fi4} illustrates the evolution of an initial
perturbation defined as a single loop of sinus.
%
\begin{figure}[htb]
\begin{center}
\includegraphics[height=7.7cm]{fig5.eps}
\caption{Unstable evolution of the interface perturbation
                 \label{fi4}}
\end{center}
\end{figure}

Growth of the oscillations is qualitatively equivalent to the
finger growth observed in the experimental research, as for
instance in \cite{Max}.
We have used the same numerical scheme (\ref{schema}) as for the
stable flow, but with other numerical values, corresponding to the
case when the oil is the lower fluid and the water is the upper
fluid. In particular, the value of $\alpha$ is now $1.1111... >1$.
The support of the initial perturbation (the loop of the sinus)
has been enlarged to get a better representation of the finger
growth. On the other part, it must be mentioned that, using the
same value for the gravitational time $t_{gr}$, the evolution of
the interface is much faster. Indeed, the consequent instant times
shown in Figure \ref{fi4} are respectively $0.2 \ s$, $0.4 \ s$,
$0.6 \ s$, $0.8 \ s$.
%------------

\subsubsection*{The nonlinear unstable stage}
Growth of the fingers leads to the instant when the finger
approaches to the critical value $\varphi = \alpha$. This is
possible only if $|{\partial {\varphi}\over \partial y} |
{\to}\infty$ in this point, that follows from (\ref{eq3.6}).
Hence, each finger crossing the critical line $\varphi = \alpha$
has there an infinite spatial derivative. From this property, it
follows that a finger can cross the critical line, only if it
becomes a delta-function.

These results correspond well to the experimental results obtained
by Maxworthy \cite{Max}, where a similar behavior of fingers
growth has been observed.

\subsection*{Conclusion}

Applications of the general model developed in \cite{Pan} in the
case where gravity perturbations are propagating much slower than
elastic perturbations have been investigated. In this case, the
model can be reduced to a nonlinear evolution equation with
varying sense of parabolicity. This equation becomes the
well-known Boussinesq equation if the upper fluid has not
viscosity and density and can describe the gravity unstable system
if the lower fluid is more light.


\begin{thebibliography}{99} 

\bibitem{BarEnRy}
G.I. Barenblatt, V.M. Entov and V.M. Ryzhik, {\em Theory of Fluid Flows Through
Natural Rocks}, Kluwer Academic Publishers, Dordrecht, 1990.

\bibitem{Bear}
J. Bear, {\em Dynamics of Fluid in
Porous Media}, American Elsevier, New York, 1972.

\bibitem{Pan}
C. I. Calugaru, D.-G. Calugaru, J.-M. Crolet and M. Panfilov,
Evolution of fluid-fluid interface in porous media as the model
of gas-oil fields, submitted in {\it Electron. J. Diff. Eqns.}

\bibitem{Dagan}
G. Dagan, Second-order theory of shallow free surface flow in porous media, in:
{\it Q. J. Mech. Appl. Maths} 20 (1967), 517-526.

\bibitem{LW}
P.L.-F. Liu and J. Wen, Nonlinear diffusive surface waves in porous media, in:
{\it J. Fluid Mech.} 347 (1997) 119-139.

\bibitem{Max}
T. Maxworthy, The nonlinear growth of a gravitationally unstable
interface in a Hele-Shaw cell, in: {\it J. Fluid Mech.} 177 (1987)
207-232.

\bibitem{Parl}
J.-Y. Parlange, F. Stagniti, J.L. Starr and R.D. Braddock, Free-surface flow
in porous media and periodic solution of the shallow-flow approximation, in:
{\it J. Hydrology} 70 (1984) 251-263.

\bibitem{Plot}
P.I. Plotnikov, Equations with variable direction of parabolicity
and the hysteresis effect, in: {\it Doklady Rossiiskoi Academii
Nauk} 330 (1993) no. 6, 691-693.

\bibitem{PolKoch}
P.Ya. Polubarinova-Kochina, {\em Theory of Groundwater Flow.} Nauka, 
Moscow (in Russian), 1984.

\bibitem{Whit}
G.B. Whitham, {\em Linear and Nonlinear Waves}, John Wiley \& Sons, 
New York, 1974.

\end{thebibliography}

\end{document}