\documentclass[reqno]{amsart} \usepackage{graphicx} \AtBeginDocument{{\noindent\small {\em Electronic Journal of Differential Equations}, Vol. 2003(2003), No. 73, pp. 1--13.\newline ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu \newline ftp ejde.math.swt.edu (login: ftp)} \thanks{\copyright 2003 Southwest Texas State University.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE--2003/73\hfil Evolution of fluid-fluid interface in porous media] {Evolution of fluid-fluid interface in porous media as the model of gas-oil fields} \author[C.-I. Calugaru, D.-G Calugaru, J.-M. Crolet, \& M. Panfilov\hfil EJDE--2003/73\hfilneg] {Cerasela-Iliana Calugaru, Dan-Gabriel Calugaru, Jean-Marie Crolet, \& Michel Panfilov} \address{Cerasela-Iliana 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, MCS/CDCSP, ISTIL, \\ 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} \date{} \thanks{Submitted January 7, 2003. Published June 30, 2003.} \subjclass[2000]{76B15, 76S05, 76T05, 65N06} \keywords{Porous media, two-phase flow, interface, oil reservoirs} \begin{abstract} This article proposes a generalized model for describing deformations of the mobile interface separating two immiscible weakly compressible fluids in a weakly deformable porous medium. It describes a gravity non-equilibrium processes, including evolution of the gravitational instability and can be reduced in two cases. This paper deals with the first case in which elastic perturbations are propagating much slower than gravity perturbations. The obtained model has analytical solutions and is applied to simulate the behavior of oil-gas or water-oil interface in oil-gas reservoirs. \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} This paper studies the behaviour of the mobile interface between two immiscible fluids which can be distinguished by viscosity and density. Two-phase flow with mobile interface is one of the more complicated objects in mathematical analysis and is usually described by a system of partial differential equations which are true on either side of the interface. They are linked by some dynamic and kinetic conditions on the interface. The major difficulty is the transform of such a system to some explicit closed differential equation for a selected coordinate of the mobile surface. More generally speaking, if the interface equation is $F(x,y,z;t)=0 $, or $z=h(x,y;t)$, then the problem is to deduce the closed differential equation for the function $F(x,y,z;t)$ or $h(x,y;t)$. In fluid mechanics, it can be done in three basic cases. First of all, this is the case when the viscous forces can be neglected for both fluids \cite{Whit}, which does not concern flow in porous media. Secondly, this is the case when one of two fluids has no viscosity, as for instance in case of water-air flow. The theory of shallow water \cite{Whit} yields a classical example of such an explicit model for the function $h$ describing the surface waves on water. The model for $h$ gets the form of the Korteweg-de Vries equation. Another example corresponds to the groundwater flow in unconfined aquifer \cite{BarEnRy,Bear} which leads to the nonlinear parabolic Boussinesq equation with respect to $h$. In \cite{BarEnRy} this model has been deduced by a simple integration of flow equations over the vertical coordinate $z$ and assuming the hydrostatic pressure distribution along $z$. In \cite{LW} the same model has been obtained by asymptotic expansion method, assuming the vertical size of the porous reservoir is much smaller than its horizontal length. More general models have been obtained in \cite{Dagan,Parl} where the gravity equilibrium hypothesis has been removed. When both fluids are viscous, the explicit equations were deduced in \cite{PolKoch} assuming the steady-state flow for two the fluids and that the horizontal flow velocity is constant along $z$. The first condition is often replaced by the gravity equilibrium condition. In case of two-liquid flow, the condition of gravity equilibrium, i.e., the hydrostatic pressure distribution, becomes excessively strong. It does not allow the description of development of fast gravity perturbations. In particular, it is impossible to analyze the evolution of instability within the framework of such a model. Condition of steady-state flow does not allow the study of flow of rather compressible liquids, when the elastic perturbation is propagating slower than the gravity perturbation. In this paper both fluids are assumed to be viscous and compressible. The hypo- thesis of gravity equilibrium is removed. The flow is non-stationary. \section{Formulation of the problem} \subsection*{Physical model} Let us introduce an orthogonal coordinate system ($x_1, x_2, z$), where $z$ is the ``vertical" coordinate, which has the same direction as the gravity vector, while ($x_1, x_2$) are coordinates of the horizontal plane. Let us examine a horizontal porous stratum of a height $H$, in a domain $\Omega{\subset}R^3$ where two immiscible fluids are separated from each other by an interface, which does not cross the top and bottom boundaries, as is shown in Fig.\ref{fi1}. \begin{figure}[htb] \begin{center} \includegraphics[width=0.7\textwidth]{fig1.eps} \caption{Scheme of the process} \label{fi1} \end{center} \end{figure} Let the value $h(x_1, x_2; t)$ denote the height of the interface with respect to the bottom of the stratum, and the indexes $i=I, II$ correspond to the lower and the upper fluid. %-------------------------------- \subsection*{Mathematical formulation} The porous stratum is assumed to be homogeneous, but anisotropic, in such a way that the tensor of permeability $K{\equiv}\{ K_{ij}\}_{i,j=1}^3$ is diagonal with $K_z{\equiv}K_{33}$, $K_{x_1}{\equiv}K_{11}$ and $K_{x_2}{\equiv}K_{22}$. %-------- \subsubsection*{Equations of flow} Flow of the weakly compressible $i{-}th$ fluid in a weakly elastic medium can be described by the usual system of equations with respect to the fluid pressure $P^i$ and the flow rate $\vec{V}^i$: \begin{subequations}\label{eq1.2} \begin{gather} \frac{\partial {P^i}}{\partial t}={\partial \over \partial x_{k}} \big( \mbox{\ae}_{kj} {\partial {\Phi^i} \over \partial x_{j}} \big), \label{eq1.2a} \\ \vec{V}^i=-\frac{K}{\mu^i}\mathop{\rm grad} \Phi^i \label{eq1.2b} \quad i=I, II \end{gather} \end{subequations} where $\rho$ is the fluid density, $g$ the gravity acceleration, $\mu$ is the fluid viscosity and $\Phi^i{\equiv}P^i{+}\rho^i g z$. The summation is done with respect to repeated indexes $k$ and $j$. The piezo-conductivity tensor $\mbox{\ae}_{kj}$ is defined as $\mbox{\ae}_{kj}=(K_{kj}\beta^i_*)/\mu^i$, where the parameter $\beta^i_*$ has a dimension of pressure and is a measure of fluid/medium compressibility. %-------- \subsubsection*{Cinematic equation for the interface} For further deductions, an explicit cinematic equation for the height of the mobile interface $h(x_1, x_2; t)$ can be written in the following way. Let the interface equation be: \begin{equation} \label{eq1.3} z=h(x_1, x_2; t) \end{equation} This function is assumed to exist and to be unique, therefore, formation of some loops is excluded. Assuming the interface remains always smooth, we get the following cinematic equation: \begin{equation} \label{eq1.4} \frac{\partial {h}}{\partial t}+\vec{U}(h) \mathop{\rm grad} h=U_z(h) \end{equation} % with $\vec{U}$ being the true velocity of the surface, while $grad$ being the 2-D operator: $$ \mathop{\rm grad} \equiv \frac{\partial }{\partial x_1}\vec{q}_1+ \frac{\partial }{\partial x_2}\vec{q}_2 $$ where $\vec{q}_k$ is the basis vector along the axis $x_k$. \subsubsection*{Conditions at the interface} The following three necessary conditions should be imposed on the interface: \begin{itemize} \item [i)] continuity of the pressure; \item [ii)] continuity of the normal flow velocity; \item [iii)] equivalence between velocity of the interface and physical velocity of the fluid in each point of the interface. \end{itemize} Then \begin{subequations}\label{eq1.5} \begin{gather} P^I\Big\vert_h=P^{II}\Big\vert_h{\equiv}{\mathcal P} \label{eq1.5a} \\ V^I_n\Big\vert_h=V^{II}_n\Big\vert_h \label{eq1.5b} \\ m\frac{\partial {h}}{\partial t}{+}\vec{V^I}(h)\mathop{\rm grad} h=V^{I}_{z}(h) \label{eq1.5c} \\ m\frac{\partial {h}}{\partial t}{+}\vec{V^{II}}(h) \mathop{\rm grad} h=V^{II}_{z}(h) \label{eq1.5d} \end{gather} \end{subequations} because the true flow velocity is equal to $\vec{V}/m$, with $m$ being the medium porosity. Only two from three last relations are independent. We will use the equation (\ref{eq1.5c}) and the following one, which results from (\ref{eq1.5c}) and (\ref{eq1.5d}): \begin{equation} \vec{V^I}(h)\mathop{\rm grad} h - V^{I}_{z}(h) = \vec{V^{II}}(h)\mathop{\rm grad} h - V^{II}_{z}(h) \label{eq1.6} \end{equation} \subsubsection*{Conditions at the top and bottom of the medium} The top and the bottom surfaces of the stratum, which are horizontal, are assumed to be impermeable: \begin{equation}\label{eq1.7} V^{II}_z\big\vert_{z=H}=V^{I}_z\big\vert_{z=0}=0 \end{equation} \section{Averaged equations} Since the study of such a system of equations with a mobile boundary is too complex, it is easier to deduce an explicit equation for the interface height $h(x,t)$. This can be done by integrating (\ref{eq1.2a}) over the vertical coordinate $z$ within the intervals $0{\le}z{\le}h$ and $h{\le}z{\le}H$. This technic allows a simplification of the equations in such a way that the unknowns don't depend on $z$. However, this integration leads to some loss of information. Such a loss must be restored by introducing some hypothesis about the vertical distribution of the velocity or the pressure field. An assumption of linear behavior of the vertical velocity along $z$ is used. %---------------------------------------- \subsection*{Deduction of equations averaged over the vertical coordinate} The simple integration of (\ref{eq1.2}a) over intervals $0{\le}z{\le}h$ and $h{\le}z{\le}H$ yields the following set of equations: \begin{subequations}\label{eq2.1} \begin{gather} \begin{aligned} \frac{\partial {(h\overline{P}^I)}}{\partial t} - {\mathcal P} \frac{\partial h}{\partial t} =& \mathop{\rm div}{\left( \mbox{\ae}^I_x\mathop{\rm grad}{(h\overline{\Phi}^I)} \right)}- \mbox{\ae}^I_x \frac{\partial\Phi^I}{\partial x_j}\Big\vert_h \frac{\partial h}{\partial x_j} \\ &- \mathop{\rm div}{\left( \mbox{\ae}^I_x \Phi^I(h)\mathop{\rm grad}{h}\right) }+ \mbox{\ae}^I_z \frac{\partial\Phi^I}{\partial z}\Big\vert_h \end{aligned} \label{eq2.1a} \\ \begin{aligned} \frac{\partial {(h^{II}\overline{P}^{II})}}{\partial t} - {\mathcal P}\frac{\partial h^{II}}{\partial t} =& \mathop{\rm div}{\left( \mbox{\ae}^{II}_x \mathop{\rm grad}{(h^{II}\overline{\Phi}^{II})}\right)}+ \mbox{\ae}^{II}_x \frac{\partial\Phi^{II}}{\partial x_j}\Big\vert_h \frac{\partial h}{\partial x_j} \\ &-\mathop{\rm div}{\left( \mbox{\ae}^{II}_x \Phi^{II}(h)\mathop{\rm grad}{h^{II}}\right) }- \mbox{\ae}^{II}_z \frac{\partial\Phi^{II}}{\partial z}\Big\vert_h \end{aligned} \label{eq2.1b} \\ m\frac{\partial h}{\partial t} - \frac{K_x}{\mu^I}\frac{\partial\Phi^I}{\partial x_j}\Big\vert_h \frac{\partial h}{\partial x_j} = -\frac{K_z}{\mu^I}\frac{\partial\Phi^I}{\partial z}\Big\vert_h \label{eq2.1c} \\ \frac{K_x}{\mu^I}\frac{\partial\Phi^I}{\partial x_j}\Big\vert_h \frac{\partial h}{\partial x_j} - \frac{K_z}{\mu^I}\frac{\partial\Phi^I}{\partial z}\Big\vert_h = \frac{K_x}{\mu^{II}}\frac{\partial\Phi^{II}}{\partial x_j}\Big\vert_h \frac{\partial h}{\partial x_j} - \frac{K_z}{\mu^{II}}\frac{\partial\Phi^{II}}{\partial z}\Big\vert_h \label{eq2.1d} \\ h^{II} {\equiv} H{-}h\,, \label{eq2.1e} \end{gather} \end{subequations} where \[ \overline{f}^I \equiv \frac{1}{h}\int_0^h f dz,\quad \overline{f}^{II} \equiv \frac{1}{h^{II}}\int_h^H fdz,\quad f(h) \equiv f(x_1,x_2,h) \] for any function $f(x_1,x_2,z)$, and ${\mathcal P}$ is the pressure at the interface. %---------------------------------------- \subsection*{Hypothesis about the vertical distribution of velocity} Let us assume that the vertical distribution of flow velocity is linear: \begin{gather*} V^I_z(x_1, x_2, z)=\frac{K_z\eta^I(x_1, x_2)}{\mu^I}z \\ V^{II}_z(x_1, x_2, z)=\frac{K_z\eta^{II}(x_1, x_2)}{\mu^{II}}(z{-}H) \end{gather*} where conditions (\ref{eq1.7}) have been taken into account. Usually an hypothesis on hydrostatic pressure distribution which is equivalent to zero vertical flow velocity is retained \cite{BarEnRy,Bear}. Such an assumption seems to be sufficient if the interface deformation is rather small, if the boundary is free (the upper fluid has no viscosity and density) and if the transition phenomena are not taken into account. For a rather general case studied in this paper a more general assumption is needed. Moreover, one can show that the hypothesis dealing with a zero vertical flow velocity leads to an overdetermined, contradictory system of equations. On the other hand, any other law of velocity distribution (quadratic, etc) leads to a non closed system. Using this hypothesis we can define the functions $\Phi^I$ and $\Phi^{II}$ in the following form: \begin{equation} \label{eq2.2000} \Phi^I=\Phi^I(h){-}\frac{\eta^I\big( z^2{-}h^2\big)}{2} , \quad \Phi^{II}=\Phi^{II}(h){-}\frac{\eta^{II}}{2}\left[ (z{-}H)^2{-}(h{-}H)^2 \right] \end{equation} % and, at last, the derivatives $ \frac{\partial\Phi^i}{\partial z}$ and $ \frac{\partial\Phi^i}{\partial x_j}$ at the interface in (\ref{eq2.1}): \begin{subequations}\label{eq2.2} \begin{gather} -\frac{\partial\Phi^I}{\partial z}\Big\vert_h=\eta^I h\,,\quad \frac{\partial\Phi^I}{\partial x_j}\Big\vert_h= \frac{\partial\Phi^I(h)}{\partial x_j}{+}\frac{1}{2}\eta^I \frac{\partial h^2}{\partial x_j} \label{eq2.2a} \\ -\frac{\partial\Phi^{II}}{\partial z}\Big\vert_h=\eta^{II}(h-H), \quad \frac{\partial\Phi^{II}}{\partial x_j}\Big\vert_h= \frac{\partial\Phi^{II}(h)}{\partial x_j}{+}\frac{1}{2}\eta^{II} \frac{\partial (H{-}h)^2}{\partial x_j} \label{eq2.2b} \end{gather} \end{subequations} New parameters $\eta^i$ can be defined via $\overline{\Phi}^i{-}\Phi^i(h)$ after integrating equations (\ref{eq2.2000}): \begin{equation} \label{eq2.3} h\eta^I=\frac{3}{h}\left( \overline{\Phi}^{I}{-}\Phi^{I}(h)\right)\, , \quad h^{II}\eta^{II}=\frac{3}{h^{II}}\left( \overline{\Phi}^{II}{-}\Phi^{II}(h) \right) \end{equation} %---------------------------------------- \subsection*{Closed form of the averaged equations} After substituting (\ref{eq2.2}) and (\ref{eq2.3}) in (\ref{eq2.1}), we get \begin{gather*} L^I \big( R^I{+}\frac{\rho^Igh^2}{2}\big) {+} hL^I{\mathcal P}~=~ -\frac{3\mbox{\ae}^I_z R^I}{h^2} - \mbox{\ae}^{I}_x\big( \rho^{I}g + \frac{3R^{I}}{h^2}\big) (\mathop{\rm grad} h)^2 \\ L^{II} \big( R^{II} -\frac{\rho^{II}g(h^{II})^2}{2}\big) {+} h^{II}L^{II}{\mathcal P} = \mbox{\ae}^{II}_x\big( \rho^{II}g - \frac{3R^{II}}{(h^{II})^2}\big) (\mathop{\rm grad} h)^2 - \frac{3\mbox{\ae}^{II}_z R^{II}}{(h^{II})^2} \\ \frac{\mu^I m}{K_x}\frac{\partial h}{\partial t} - \mathop{\rm grad}{\mathcal P}\mathop{\rm grad} h - \big( \rho^{I}g + \frac{3R^{I}}{h^2}\big)(\mathop{\rm grad}{h})^2 - \frac{3\varepsilon_z R^I}{h^2}=0 \\ \begin{aligned} (1{-}\overline{\mu})\mathop{\rm grad}{\mathcal P}\mathop{\rm grad}{h} + 3\varepsilon_z \big( \frac{R^I}{h^2} + \overline{\mu} \frac{R^{II}}{(h^{II})^2}\big) &\\ + \Big[ \rho^{I}g + \frac{3R^I}{h^2} - \overline{\mu} \big( \rho^{II}g- \frac{3R^{II}}{(h^{II})^2}\big)\Big] (\mathop{\rm grad} h)^2 &= 0\,, \end{aligned} \end{gather*} where \begin{gather*} L^i \equiv \frac{\partial }{\partial t} - \mbox{\ae}^i_x\Delta ,\quad i = I,II \\ R^I \equiv h\big( \overline{\Phi}^I - \Phi^I(h)\big)\,,\quad R^{II}{\equiv}h^{II}\big( \overline{\Phi}^{II} - \Phi^{II}(h)\big) \end{gather*} Then, four equations define four functions, $\overline{\Phi}^I$, $\overline{\Phi}^{II}$, $\mathcal P$ and $h$. %---------------------------------------- \subsection*{Small deformations of the interface} Assuming the deformations of the interface and the lateral pressure gradients are small, and neglecting the small values of second order, we obtain the simplified system: \begin{equation}\label{eq2.400} \begin{gathered} L^I \big( R^I{+}\frac{\rho^Ig}{2}h^2\big) + hL^I{\mathcal P} = -\frac{3\mbox{\ae}^I_z R^I}{h^2} \\ L^{II} \big(R^{II}{-}\frac{\rho^{II}g}{2}(h^{II})^2\big) + h^{II}L^{II}{\mathcal P} = - \frac{3\mbox{\ae}^{II}_z R^{II}}{(h^{II})^2} \\ \frac{\mu^I m}{K_x}\frac{\partial h}{\partial t} = \frac{3\varepsilon_z R^I}{h^2}\,, \quad \frac{R^I}{h^2} {+} \overline{\mu} \frac{R^{II}}{(h^{II})^2}=0 \end{gathered} \end{equation} This system can be written in the following dimensionless form: \begin{subequations}\label{eq2.5} \begin{gather} {\mathcal L}^I\big( \varphi^2{+}\frac{\lambda_1^2\omega\tau_*}{3\varepsilon_z} \varphi^2{\partial {\varphi}\over \partial \tau}\big) + \lambda_2\varphi{\mathcal L}^I\xi = -\omega\tau_* {\partial {\varphi}\over \partial \tau} \label{eq2.5a} \\ {\mathcal L}^{II}\big( {-}\psi^2 + \frac{\lambda_1^2\overline{\rho}\omega\tau_*} {3\varepsilon_z\overline{\mu}\lambda_0} \psi^2{\partial {\psi}\over \partial \tau}\big) + \lambda_2\lambda_0\overline{\rho}\psi{\mathcal L}^{II}\xi = -\frac{\overline{\rho}\lambda_0\omega\tau_*}{\overline{\mu}} {\partial {\psi}\over \partial \tau} \label{eq2.5b} \\ \psi=-\lambda_0\varphi{+}\lambda_0{+}1 \label{eq2.5c} \end{gather} \end{subequations} where new variables are denoted as: $$ \varphi{\equiv}\frac{h}{h_0},\quad \psi{\equiv}\frac{h^{II}}{h^{II}_0},\quad \xi{\equiv}\frac{\mathcal P}{\mathcal P}_0,\quad \tau{\equiv}\frac{t}{t_*},\quad y{\equiv}\frac{x}{L} $$ where $L$ is the horizontal scale of the domain; $h_0$, $h^{II}_0$ are the heights of the lower and the upper layers in an unperturbed state; ${\mathcal P}_0$ is the pressure at the interface in the unperturbed state. The new operators are: $$ {\mathcal L}^I{\equiv}\tau_*{\partial \over \partial \tau}{-}\Delta_{yy} \,,\quad {\mathcal L}^{II}{\equiv}\frac{\tau_*\overline{\beta}}{\overline{\mu}} {\partial \over \partial \tau}{-}\Delta_{yy} $$ where the symbol $\Delta_{yy}$ denotes Laplace's operator written via variable $y$. The following set of parameters defines the process: % \begin{equation} \label{eq2.6} \begin{gathered} \lambda_0{\equiv}\frac{h_0}{h^{II}_0},\quad \lambda_1{\equiv}\frac{h_0}{L},\quad \lambda_2{\equiv}\frac{2{\mathcal P}_0}{\rho^Igh_0},\quad \omega{\equiv}\frac{2m\beta^I_*}{\rho^Igh_0}=\frac{t_{gr}}{t_{el}},\quad \tau_*{\equiv}\frac{t_{el}}{t_*} \\ \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} \end{gathered} \end{equation} Two characteristic times have been introduced: $$ 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}, $$ where $t_{el}$ defines the time of propagation of an elastic perturbation within the scale $L$, while $t_{gr}$ is the time of complete extraction of the fluid from the medium due to the gravity drop. The time $t_*$ may be chosen in two various ways: \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, where the elasticity of fluid/medium can be neglected respectively to the gravity governed motion, will be examined in section \ref{sec4}. %-------------------- \subsection*{Relation for flow rates and averaged pressures} For further consideration the relation for flow rates averaged over the layer thickness will be necessary to set boundary conditions. Let us examine any cylindrical surface ${\mathcal F}$ orthogonal to the plane $(x,y)$ and intersecting the top and the bottom of the domain $\Omega$. Let ${\mathcal G}$ be a closed plate line which results as the intersection of the surface ${\mathcal F}$ with any orthogonal horizontal plane. Volume flow rate of the upper and the lower fluids across the interface ${\mathcal F}$ is defined as: $$ Q^{I} {\equiv}\int_{\mathcal G}\int_0^h V_ndzd{\mathcal G} \,,\quad Q^{II}{\equiv}\int_{\mathcal G}\int_h^H V_ndzd{\mathcal G} $$ where $V_n$ is the component of flow velocity normal to ${\mathcal F}$. Let us introduce the dimensionless densities of the flow rates across ${\mathcal F}$ as: $$ q^{I} {\equiv}\frac{t_*}{h_0L}\int_0^h V_ndz\,,\quad q^{II}{\equiv}\frac{t_*}{h^{II}_0 L}\int_h^H V_ndz $$ which are related to $Q^i$ as \begin{equation}\label{eqdebit} Q^{I} =\frac{h_0L}{t_*}\int_{\mathcal G}q^I d{\mathcal G}\,,\quad Q^{II} =\frac{h^{II}_0L}{t_*}\int_{\mathcal G}q^{II} d{\mathcal G} \end{equation} For dimensionless flow rates it is easy to get the following relations via variables $\varphi$, $\psi$, $\xi$: \begin{equation} \label{eq2.67} \begin{gathered} q^I=-\frac{1}{\omega\tau_*}\Big\{ \frac{\partial}{\partial n} \big( \varphi^2{+} \frac{\lambda_1^2\omega\tau_*}{3\varepsilon_z}\varphi^2{\partial {\varphi} \over \partial \tau} \big) {+} \lambda_2\varphi\frac{\partial\xi }{\partial n} \Big\}\\ q^{II}=-\frac{\overline{\mu}} {\omega\tau_*\overline{\rho}\lambda_0^2} \Big\{ \frac{\partial}{\partial n}\big( {-}\psi^2 {+} \frac{\lambda_1^2\overline{\rho}\omega\tau_*} {3\varepsilon_z\overline{\mu} \lambda_0}\psi^2 {\partial {\psi}\over \partial \tau}\big) {+} \lambda_2\lambda_0\overline{\rho}\psi\frac{\partial\xi}{\partial n}\Big\} \end{gathered} \end{equation} where $\partial{/}\partial n$ means the derivation along the normal direction to the surface ${\mathcal F}$. For the averaged pressures $\overline{P}^I$ and $\overline{P}^{II}$ the following relations are true by definition: \[ h\overline{P}^I =h{\mathcal P}{+}\frac{\rho^Igh^2}{2} + R^I \, ,\quad h^{II}\overline{P}^{II}=h^{II}{\mathcal P} - \frac{\rho^{II}g(h^{II})^2}{2} + R^{II} \] and then we get from (\ref{eq2.400}): \[ h\overline{P}^I=h{\mathcal P}{+}\frac{\rho^Igh^2}{2}{+} \frac{h^2\mu^Im}{3\varepsilon_zK_x}\frac{\partial h}{\partial t} \,, \quad h^{II}\overline{P}^{II}=h^{II}{\mathcal P}{-} \frac{ \rho^{II}g(h^{II})^2 }{2}{-} \frac{(h^{II})^2\mu^Im}{3\varepsilon_zK_x\overline{\mu}}\frac{\partial h}{\partial t}, \] or in the dimensionless form: \begin{equation} \begin{gathered} \varphi p^I= \varphi^2{+}\frac{\lambda_1^2\omega\tau_*}{3\varepsilon_z} \varphi^2{\partial {\varphi}\over \partial \tau}{+}\lambda_2\varphi\xi , \\ \psi p^{II} = {-}\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}\psi\xi\,, \end{gathered} \label{eq2.500} \end{equation} where $$ p^I{\equiv}\frac{2\overline{P}^I}{\rho^Igh_0},\quad p^{II}{\equiv}\frac{2\overline{P}^{II}}{\rho^{II}gh^{II}_0} $$ Note the function $\xi $ can be excluded from (\ref{eq2.500}), if we multiply the first equation by $\lambda_0\overline{\rho}\psi$ , the second equation by $\varphi$, and substract the second equation from the first one. \begin{equation} \lambda_0\overline{\rho}p^I{-}p^{II}= \lambda_0\big( \overline{\rho}{-}1\big)\varphi{+}1{+}\lambda_0{+} \frac{\lambda_1^2\lambda_0\overline{\rho}\omega\tau_*}{3\varepsilon_z} \Big( \varphi{+}\frac{1}{\overline{\mu}\lambda_0}\big[1{+}\lambda_0{-} \lambda_0\varphi\big] \Big) {\partial {\varphi}\over \partial \tau} \label{eqpress} \end{equation} %--------------------------------------- \subsection*{Partial linearization} The condition of small perturbation being accepted, the following simplification is justified. In the left-hand part of system (\ref{eq2.5}), the second term can be linearized assuming that $\varphi{\simeq}1$, $\psi{\simeq}1$. Then system (\ref{eq2.5}) gets the form \begin{subequations}\label{eq2.8} \begin{gather} {\mathcal L}^I\big( \varphi^2{+}\frac{\lambda_1^2\omega\tau_*}{3\varepsilon_z} \varphi^2{\partial {\varphi}\over \partial \tau} {+} \lambda_2\xi \big) = -\omega\tau_*{\partial {\varphi}\over \partial \tau} \label{eq2.8a} \\ {\mathcal L}^{II}\big( {-}\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 \big) = -\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} In the next section only this partially linearized system of equations will be studied. Relations (\ref{eq2.67}) for the dimensionless flow rates across any vertical surface ${\mathcal F}$ take the form \begin{gather} \label{eq2.68} q^I=-\frac{1}{\omega\tau_*}\frac{\partial}{\partial n} \big( \varphi^2{+} \frac{\lambda_1^2\omega\tau_*}{3\varepsilon_z}\varphi^2{\partial {\varphi}\over \partial \tau} {+} \lambda_2\xi \big) \\ q^{II}=-\frac{\overline{\mu}} {\omega\tau_*\overline{\rho}\lambda_0^2} \frac{\partial}{\partial n}\big( {-}\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 \big) \end{gather} Note that the function $\xi$ might be excluded from these relations: \begin{equation} \begin{aligned} &\frac{\overline{\rho}\omega\tau_*}{\overline{\rho}{+}\lambda_0} \big( \frac{\lambda_0}{\overline{\mu}}q^{II}{-}q^I\big)\\ &= \frac{\partial}{\partial n}\Big( \varphi^2{-} \frac{2(\lambda_0{+}1)}{\overline{\rho}{+}\lambda_0}\varphi{+} \frac{\lambda_1^2\overline{\rho}\omega\tau_*}{3\varepsilon_z ({\overline{\rho}{+}\lambda_0}) } \big( \varphi^2{+}\frac{1}{\lambda_0\overline{\mu}}\big[ 1{+}\lambda_0{-} \lambda_0\varphi \big]^2\big) {\partial {\varphi}\over \partial \tau}\Big) \label{eqdebit2} \end{aligned} \end{equation} %############################################################ \section{Slow elastic perturbations} \label{sec4} Examine the case of strongly deformable fluids and medium, where the time of elastic wave propagation is large with respect to the gravity time ($\omega{\ll}1$). It is necessary to note that such situation may happen when the layer is rather thin, i.e., the condition $\lambda_1{\ll}1$ should be added. The scale of the time $t_*$ should be chosen as equal to $t_{el}$, according to (\ref{eq2.7}). Then $\tau_*=1$. Therefore, we get from (\ref{eq2.8}): \begin{equation} \label{eq4.1} \begin{gathered} {\mathcal L}^I\left( \varphi^2{+} \lambda_2\xi \right)= 0,\quad {\mathcal L}^{II}\left( {-}\psi^2 {+} \lambda_2\lambda_0\overline{\rho} \xi \right) = 0,\quad \psi={-}\lambda_0\varphi{+}\lambda_0{+}1 \\ {\mathcal L}^I{\equiv}{\partial \over \partial \tau}{-}\Delta_{yy},\quad {\mathcal L}^{II}{\equiv} \frac{\overline{\beta}}{\overline{\mu}} {\partial \over \partial \tau}{-}\Delta_{yy} \end{gathered} \end{equation} Let us introduce the new functions \begin{equation} S^I{\equiv}\varphi^2{+} \lambda_2\xi \quad \mbox{and} \quad S^{II}{\equiv}{-}\psi^2 {+} \lambda_2\lambda_0\overline{\rho} \xi \label{eq4.2} \end{equation} % which represent the averaged pressure over the corresponding layer height, as it follows from (\ref{eq2.500}). Then, the system (\ref{eq4.1}) can be written in the following form \begin{equation} \label{eq4.3} \begin{gathered} {\partial {S^I}\over \partial \tau}{-}\Delta_{yy} S^I=0, \quad \frac{\overline{\beta}}{\overline{\mu}} {\partial {S^{II}}\over \partial \tau}{-}\Delta_{yy} S^{II} =0, \\ \varphi=\alpha\Big[ 1{+} \sqrt{ \frac{ \overline{\rho}\lambda_0S^I{-}S^{II} } { \alpha\lambda_0\big( 1{+}\lambda_0\big) } {-} \frac{ \overline{\rho} }{ \lambda_0 } } \Big]\,, \end{gathered} \end{equation} where ${\alpha{\equiv}(1{+}\lambda_0)/(\overline{\rho}{+}\lambda_0)}$. By solving two linear parabolic equations with respect to $S^I$ and $S^{II}$, the function $\varphi$ is obtained as a simple solution of a quadratic equation. The function $\xi$ may be found as: ${\xi=\frac{1}{\lambda_2}\big( S^I{-}\varphi^2\big)} $. %------------------------------------------------- \subsection*{Problem of oil-water or oil-gas interface} Let us examine the problem of oil extraction from a porous reservoir by a well in the framework of the model (\ref{eq4.3}). The lower layer is saturated by water (index $I$). Construction of the well is assumed to be multitube, in such a way that each fluid can be extracted separately one from other through its own tube. The engineering problem consists of controlling deformations of the oil-water interface in order to reduce extraction of water. Such technology has been analyzed, for instance, in \cite{Zak}. The similar problem arises in a gas-oil system. Then, the indexes $I$ and $II$ are associated to oil and gas correspondingly. Examine the following problem of radial flow towards a single well located in the center of a cylindrical porous domain with the radius $R_*$, the height $H$, the porosity $m$ and the permeability $K$. Let the well be a vertical cylinder of radius $R_w$. Let $Q^i$, ($i=I,II$) be the volumic flow rate of extraction of the $i$-th fluid by the well, which are specified. Let us assume the normal flow velocity at the well border does not depend on the polar angle. Then, using (\ref{eqdebit}), we get the following equation relating $Q^i$ with $q^i|_{r=r_w}$, where $r^2=y_1^2{+}y_2^2$, $r_w=R_w/L$: % \begin{equation} \label{eqdebit4} Q^I= 2\pi R_w\int\limits^h_0 V_n\vert_{r=r_w} dz = \frac{ 2\pi R^2_* h_0 }{t_{el}}r_wq^I\big\vert_{r=r_w}, \quad Q^{II} = \frac{2\pi R^2_* h^{II}_0}{t_{el}}r_wq^{II}\big\vert_{r=r_w} \end{equation} Then using (\ref{eq4.3}), we obtain \begin{equation}\label{eq5.1} \begin{gathered} {\partial {S^I}\over \partial \tau}{-}\frac{1}{r}{\partial \over \partial r} \big( r{\partial {S^I}\over \partial r}\big) = 0, \quad \frac{\overline{\beta}}{\overline{\mu}} {\partial {S^{II}}\over \partial \tau}{-}\frac{1}{r}{\partial \over \partial r} \big( r{\partial {S^II}\over \partial r}\big) = 0, \\ S^I\big\vert_{\tau =0}=1{+}\lambda_2, \quad S^{II}\big\vert_{\tau =0}={-}1{+}\lambda_2\lambda_0\overline{\rho}, \\ S^I\big\vert_{r{\to}r_*}=1{+}\lambda_2, \quad S^{II}\big\vert_{r{\to}r_*}={-}1{+}\lambda_2\lambda_0\overline{\rho}, \\ r\frac{\partial S^I}{\partial r}\big\vert_{r=r_w}=\frac{\omega\varepsilon m}{2}, \quad r\frac{\partial S^{II}}{\partial r}\big\vert_{r=r_w}= \frac{\omega\varepsilon\gamma m\lambda_0}{2\overline{Q}}\,, \end{gathered} \end{equation} where $$ r=\frac{R}{R_*}, \quad \varepsilon={t_{el}\over T}={Q^I\mu^I\over\pi h_0 m K_x\beta^I_*}, \quad \overline{Q}={Q^I\over Q^{II}}, \quad T={\pi R^2_* h_0 m\over Q^I}, \quad \gamma=\frac{\overline{\rho}\lambda_0^2}{\overline{\mu}}\,, $$ $T$ is the time of full extraction of the fluid $I$ by the well. This value defines some new ``technological" time scale of the system. Boundary conditions at the well border represent the fixed flow rates of each fluid and are written using relations (\ref{eq2.68}). After this problem is solved, the interface height $\varphi (r,\tau )$ can be determined using the last equation in (\ref{eq4.3}). %----------------- \subsubsection*{Self-similar solution} When the reservoir is infinite $(R_*{\to}\infty)$, and the well radius is zero $(R_w{\to} 0)$, then the problem (\ref{eq5.1}) written in dimensional variables has the exact analytical solution in term of $S^I=S^I(\xi )$, $S^{II}=S^{II}(\xi )$, $\xi{\equiv}R/\sqrt{t}$. The problem takes the form \begin{gather*} - {\frac{t_*}{R^2_*}}{\xi\over 2}{dS^I\over d\xi}={1\over\xi}{d\over d\xi} \big(\xi{dS^I\over d\xi}\big) \\ - {\frac{t_*}{R^2_*}}{\xi\over 2}{\overline{\beta}\over\overline{\mu}}{dS^{II}\over d\xi} = {1\over\xi}{d\over d\xi}\big(\xi{dS^{II}\over d\xi}\big) \\ S^I\big\vert_{\xi=\infty}=1+\lambda_2,\quad S^{II}\big\vert_{\xi=\infty}=-1+\lambda_2\lambda_0\overline{\rho} \\ \xi{dS^I\over d\xi}\big\vert_{\xi=0}={\omega\varepsilon m\over 2},\quad \xi{dS^{II}\over d\xi}\big\vert_{\xi=0} = {\omega\varepsilon m\gamma\lambda_0\over 2\overline{Q}} \end{gather*} The solutions have the form \begin{gather*} S^I={\omega\varepsilon m\over 4}Ei\big(-{\xi^2 t_*\over 4R^2_*}\big){+}1{+} \lambda_2, \\ S^{II}={\omega\varepsilon m\gamma\lambda_0\over 4\overline{Q}}Ei \big(-{\xi^2\overline{\beta}t_*\over 4\overline{\mu}R^2_*}\big) {+}\lambda_2\lambda_0\overline{\rho}{-}1\,, \end{gather*} where $Ei$ is the integral exponential function defined as $$ Ei(x){\equiv}\int^x_{-\infty}{e^u\over u}du $$ Using the property of this function, we get for $\xi{\to} 0$: \begin{gather*} S^I\sim {\omega\varepsilon m\over 4} \big[\hbox{ln}{\xi^2 t_*\over 4R^2_*}{+}Ce{+}\dots \big] {+}1{+}\lambda_2\,,\\ S^{II}\sim{\omega\varepsilon m\gamma\lambda_0\over 4\overline{Q}} \big[\hbox{ln}{\xi^2\overline{\beta}t_*\over 4\overline{\mu}R^2_*}{+}Ce{+} \dots\big]{+}\lambda_2\lambda_0\overline{\rho}{-}1 \end{gather*} where $Ce=0.5772\dots$ is the Euler constant. Then the function $\varphi$ is \[ \frac{\varphi}{\alpha}=1{+} \sqrt{ \frac{{\overline{\rho}\lambda_0{+}1{+} \frac{\varepsilon\omega m}{4} \big[ \overline{\rho}\lambda_0 Ei\big( -{\xi^2 t_*\over 4R^2_*}\big) {-} \frac{\gamma\lambda_0}{\overline{Q}} Ei\big( -{\xi^2\overline{\beta}t_*\over 4\overline{\mu}R^2_*}\big) \big] }}{ \alpha\lambda_0\big( 1{+}\lambda_0\big) } {-}\frac{ \overline{\rho} }{ \lambda_0 }} \] %----------------- \subsubsection*{Numerical solution} A numerical solution of problem (\ref{eq5.1}) has been obtained. For the space discretization, a conventional scheme of order 2 obtained by a finite difference method has been used. For the time discretization a $\theta$ - scheme has been considered. In the presented numerical tests the value $\theta = 0.5$ has been used, which corresponds to the Cranck-Nickolson scheme. For all numerical tests we supposed that the initial position of the interface is horizontal and it corresponds to the unperturbed interface, such that $\varphi(r,0){\equiv}1$ (i.e. $h_0(\cdot)=h(\cdot,0){\equiv}20 m$). The physical dimensions of the reservoir are given by its depth ($H = 40 m$) and its radius ($100 m$). Figures 2, 3 and 4 show the evolution of the oil-water interface for different values of flow rates. In each case, the interface is presented at four time instants: 1 day, 20 days, 40 days, 60 days. For the first numerical test, the flow rates of the upper liquid (oil) and of the lower liquid (water) are considered to be equal, i.e. $Q^I = Q^{II} = 1000 m^3/day$. The pumping of oil leads to an ascent of water from below (fig. 2) even if water is extracted with the same flow rate. This characteristic of the process can be explained by the difference of viscosity (the viscosity ratio $\mu^{II}=10\mu^{I}$ has been used for all tests). \begin{figure}[htb] \begin{center} \includegraphics[width=10cm, height=6cm]{fig2.eps} \caption{Evolution of the interface oil-water for $Q^I = Q^{II} = 1000$ m$^3$/day} \label{fi2} \end{center} \end{figure} In fact, using the previous analytical solution, it can be proved the existence of a critical value (denoted $\overline{Q}^*$) of the ratio between the flow rates. If the flow rates verify $Q^I/Q^{II}< \overline{Q}^*$, then one can predict that, around the well, the interface evolves above its initial position. In the contrary case, the interface evolves below the initial position. When the initial layers are equal (as in our tests : $h_0 = h_0^{II} = 20m$), this critical value becomes $\mu^{II}/\mu^{I}$. A complete analysis dealing with the stability of the interface will be described in a forthcoming paper. The evolution shown in fig. 3 corresponds to the critical value $Q^I/Q^{II} = 10 = \overline{Q}^*$. One can observe that the interface remains relatively close to the initial state and tends towards this state. Therefore, the evolution of the interface can be easily controlled, but the amount of pumped water has to be ten times larger than the amount of pumped oil. \begin{figure}[htb] \begin{center} \includegraphics[width=10cm, height=6cm]{fig3.eps} \caption{Evolution of the interface oil-water for} \centerline{$Q^I = 1000$ m$^3$/day, $Q^{II} = 100$ m$^3$/day} \label{fi3} \end{center} \end{figure} In the last test, the water flow rate is maintained at 1000 m$^3$/day, but the oil flow rate is increased: $Q^{II} = 10000$ m$^3$/day. In fig. 4, one can observe an evolution similar to the first test. However, the deformation is more significant : the interface approaches the higher limit of the reservoir. The quantity of oil which was already pumped is significant, but the well can become not exploitable because of the complete covering of the oil layer by water near the well. \begin{figure}[htb] \begin{center} \includegraphics[width=10cm, height=6cm]{fig4.eps} \caption{Evolution of the interface oil-water for} \centerline{$Q^I = 1000$ m$^3$/day, $Q^{II} = 10000$ m$^3$/day} \label{fi4} \end{center} \end{figure} \section*{Conclusion} A new model to simulate interface evolution between two fluids in porous media has been developed. It generalizes previous models by removing the classical condition of hydrostatic pressure distribution and considering a non-stationary flow. Instead of these hypothesis, an assumption of linear behavior of the vertical velocity along vertical coordinate is used. The application of this model to simulate the oil-water or gas-oil interface deformations in oil reservoir is shown. The physical parameters of this problem are supposed to verify that gravity perturbations are propagating much faster than elastic perturbations. Then, the model consists of two linear diffusion equations respectively to the two averaged fluid pressures, while the vertical coordinate of the interface is related with them by a nonlinear algebraic equation. However it can describe many other real situations, including the classical well-known cases of groundwater flow with free surface. It is able also to take into account more complex phenomena, as gravitational instability with finger growth. These different applications of the model are studied in \cite{Calu}. \begin{thebibliography}{00} \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{Calu} C.-I. Calugaru, D.-G. Calugaru, J.-M. Crolet, M. Panfilov, {\em Stability of the interface between two immiscible fluids in porous media}, Electronic Journal of Differential Equations, Vol. 2003(2003), N. 25, 1-10 \bibitem{Dagan} G. Dagan, {em Second-order theory of shallow free surface flow in porous media}, Q. J. Mech. Appl. Maths., 20 (1967), 517-526. \bibitem{Zak} Yu. P. Korotaev, Yu.P. and S.N. Zakirov, {\em Theory and designing of gas and gas-condensate reservoir exploitation}, Nedra, Moscow (in Russian), 1981. \bibitem{LW} P. L.-F. Liu and J. Wen, {\em Nonlinear diffusive surface waves in porous media}, J. Fluid Mech., 347 (1997) 119-139. \bibitem{Parl} J.-Y. Parlange, F. Stagniti, J. L. Starr, and R.D. Braddock, {\em Free-surface flow in porous media and periodic solution of the shallow-flow approximation}, J. Hydrology, 70 (1984) 251-263. \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}