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

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2017 (2017), No. 191, pp. 1--5.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2017 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2017/191\hfil Flow in porous media]
{Dynamics of three dimensional flow in \\ porous media}

\author[Nikola Konatar \hfil EJDE-2017/191\hfilneg]
{Nikola Konatar}

\address{Nikola Konatar \newline
Faculty of Mathematics and Natural Sciences,
University of Montenegro,
Cetinjski put bb, 81000 Podgorica, Montenegro}
\email{nikola.k@ac.me}

\thanks{Submitted April 18, 2017. Published July 30, 2017.}
\subjclass[2010]{76S05, 76T10}
\keywords{Three dimensional flow; Darcy law; immiscible flow}

\begin{abstract}
 We describe the dynamics of the interface between two immiscible fluids
 of different densities in three dimensional porous media.
 This generalizes previous results given in two dimensions in which
 existence of the stream function played a substantial role. The most
 important contribution of the paper is the introduction of a method which
 avoids usage of the stream function.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\allowdisplaybreaks

\section{Introduction}

Problems of flow in porous media are unavoidable in most of the branches 
of industry or science. We mention environmental engineering, petroleum research, 
traffic flow, sedimentation processes, radar shape-from-shading problems, 
blood flow etc. In the current contribution, we investigate dynamics of the 
interface between two (almost) immiscible fluids governed by Darcy law.

Our method is based on the level set procedure \cite{DO, MMN, MNK, Seth} 
in the frame of which we track the interface points. The motion of the interface 
is governed by a transport equation and the points move along characteristics. 
Although the velocity which determines the characteristics are unknown, 
they can be expressed via appropriate Green type functions \cite{KMN,MMN}.

An important property of the method to be presented here is the fact that
 we succeeded in avoiding the  use of the stream function (which does not
 exist for dimensions greater than two) by, roughly speaking, replacing 
it by the pressure function. The stream function is actually a potential 
corresponding to a two dimensional divergence free vector field 
(expressing the immiscibility of the involved liquids; see \eqref{darcy3} below). 
It played a substantial role in e.g. \cite{DO, MMN, TNL} where the research 
similar to ours is conducted. From the same reason, in e.g. \cite{Tch},
onset of convection in a gravitationally unstable diffusive boundary layer 
in porous media were described only in two-dimensional situation. 
We believe that our method has a potential for applications in different 
research directions when it comes to investigations in dimensions more than two.

In the next section, we use the method to prove rigorously that a tip of the 
interface directed downwards moves down if the heavier liquid is up and vice versa.

\section{The Darcy law}

Dynamics of interface between two immiscible fluids of different densities
in porous media is governed by the Darcy law, conservation of mass, and 
immiscibility (in three-dimensional environment):
\begin{gather}
\label{darcy1}
\mathbf{u} =-\frac{K}{\mu}(\nabla P - \rho g),\quad \mathbf{u}=(u,v,w),\\
\label{darcy2}
\partial_t(\Phi \rho)+\operatorname{div}(\rho \mathbf{u})=0,\\
\label{darcy3}
\operatorname{div}(\mathbf{u})=0.
\end{gather}
Here, $K$ represents a positive definite permeability matrix, $\mu$ is 
the viscosity of the fluids, and $\Phi$ is the porosity constant of the porous 
media. We assume that these parameters are constant, and they can be 
normalized to unity. $P$ is the pressure, $\rho$ is the density of fluids, 
and $g$ is vertically directed gravity. The divergence is taken with respect 
to the variables $x$, $y$ and $z$.

Next, we will assume that the flow of the fluid is periodic with respect to 
the $x$ and $y$ with periods $2L_1$ and $2L_2$, respectively, and that the 
interface between the fluids in the moment $t$ can be described by the set
\begin{equation*}
\Gamma_t=\{(x,y,z)\in [-L_1,L_1]\times [-L_2,L_2] \times [0,1] |\, \phi^t(x,y,z)=0 \},
\end{equation*}
where $\phi^t: [-L_1,L_1]\times [-L_2,L_2]\times [0,1] \to \mathbb{R}$. 
The initial interface position is given by $\Gamma_0=\{(x,y,z): \phi^0(x,y,z)=0 \}$.
 Using this we can define the density of the fluids as (in the sequel, $\rho_h$ and 
$\rho_l$ are constants):
\begin{equation}
\label{density}
\rho=\rho_l+(\rho_h-\rho_l)H(\phi^t(x,y,z)), \quad \rho_h>\rho_l.
\end{equation}
with $\rho_l$ and $\rho_h$ representing the densities of the lighter and 
heavier fluid respectively, and $H$ is the Heaviside function. From this 
definition we can see that the heavier fluid is in the area where $\phi^t>0$, 
while the lighter is in the area where $\phi^t<0$.

We will assume that the speed of fluids vanishes at boundaries
\begin{equation}
\label{Boundaryconds}
\mathbf{u}|_{x=-L_1}=\mathbf{u}|_{x=L_1}=\mathbf{u}|_{y=-L_2}=\mathbf{u}|_{y=L_2}
= \mathbf{u}|_{z=0}=\mathbf{u}|_{z=1}=0
\end{equation}
With such assumptions, it is natural to assume hydrostatic pressure at boundaries:
\begin{equation}
\label{pressurebound}
P|_{x=\pm L_1}=P|_{y=\pm L_2}=\rho_l g (1-z).
\end{equation}
Note that we actually assume that there is no heavier fluid at the boundary.

We rewrite equations \eqref{darcy1}, \eqref{darcy2} and \eqref{darcy3} to get:
\begin{gather}
\label{darcy1u}
 u=-\frac{\partial \tilde{P}}{\partial x},\\
\label{darcy1v}
 v=-\frac{\partial \tilde{P}}{\partial y},\\
\label{darcy1w}
 w=-\frac{\partial \tilde{P}}{\partial z}+\tilde{\rho}g,\\
\label{darcy2new}
 \frac{\partial \tilde{\rho}}{\partial t}
+u \frac{\partial \tilde{\rho}}{\partial x}
+v \frac{\partial \tilde{\rho}}{\partial y} 
+ w\frac{\partial \tilde{\rho}}{\partial z}=0, \\
\label{darcy3new}
 \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}
+\frac{\partial w}{\partial z}=0,
\end{gather}
with $u$, $v$ and $w$ being the speeds of the fluid in the direction of the $x$, 
$y$ and $z$ axis respectively, and $\tilde{P}=P-\rho_l g (1-z)$, 
$\tilde{\rho}=(\rho_h-\rho_l)H(\phi^t)$. We can see from \eqref{pressurebound} 
that $\tilde{P}=0$ at the side boundaries.

Inserting \eqref{density} into \eqref{darcy2new} we obtain
\begin{equation}
\label{darcyinterface}
\frac{\partial \phi^t}{\partial t}+u|_\Gamma \frac{\partial \phi^t}{\partial x}
+v|_\Gamma \frac{\partial \phi^t}{\partial y} 
+ w|_\Gamma \frac{\partial \phi^t}{\partial z}=0.
\end{equation}
Looking at the characteristics of this equation we obtain $\dot{x}=u$,
 $\dot{y}=v$ and $\dot{z}=w$, with $(x,y,z)$ being the position of particles 
on the interface $\Gamma_t$. The initial conditions are given by $x(0)=x_0$, 
$y(0)=y_0$, $z(0)=z_0$, $(x_0,y_0,z_0) \in \Gamma_0$.

Next, we differentiate equation \eqref{darcy1u} with respect to $x$, 
\eqref{darcy1v} with respect to $y$, \eqref{darcy1w} with respect to $z$ and 
add the resulting equation we obtain:
\begin{equation}
\label{darcyfinal1}
0=\operatorname{div} \mathbf{u}
=-\Delta \tilde{P} + \frac{\partial \tilde{\rho}g}{\partial z}.
\end{equation}

Now introduce a non-decreasing function $\omega\in C^\infty(R)$ such that
$$
\omega(\xi)=
\begin{cases}
0, & \xi\leq 0\\
1, & \xi\geq 1
\end{cases}.
$$

It is well known that $H_\varepsilon(x)=\omega(\frac{x}{\varepsilon})$ 
converges weakly to the Heaviside function $H$ as $\varepsilon \to 0$. 
Also, $\delta_\varepsilon(x)=\frac{1}{\varepsilon}
H_\varepsilon'(\frac{x}{\varepsilon})\rightharpoonup \delta (x)$ as 
$\varepsilon \to 0$.

Now, take instead of $\tilde{\rho}$ in \eqref{darcy1u}-\eqref{darcy3new} the function
$$
\tilde{\rho}=(\rho_h-\rho_l)H_\varepsilon (\phi^t(x,y,z)),
$$ and denote the velocity and pressure accordingly. In particular, 
instead of \eqref{darcyfinal1} we have
\begin{equation}
\label{darcyfinalconv}
\Delta \tilde{P}_\varepsilon = \frac{\partial \tilde{\rho}_\varepsilon}{\partial z}.
\end{equation}

We introduce a sequence of smooth functions 
$G_\varepsilon:[-L_1,L_1] \times [-L_2,L_2] \times [-1,1] \to \mathbb{R}$ such that:
\begin{gather}
\label{eq1}
-\Delta G_\varepsilon=\delta_\varepsilon(x)\delta_\varepsilon(y)
\delta_\varepsilon'(z),\\
\label{eq2}
G_\varepsilon|_{x=-L_1}=G_\varepsilon|_{x=L_1}=G_\varepsilon|_{y=-L_2}
=G_\varepsilon|_{y=L_2}=G_\varepsilon|_{z=-1}=G_\varepsilon|_{z=1}=0.
\end{gather}
It is not difficult to see that $G_\varepsilon$ is an odd function with 
respect to $z$, and therefore we can extend it periodically on the entire 
$\mathbb{R}$ with respect to the $z$ direction so that the extension remains 
a solution to \eqref{eq1}.

This sequence of functions converges in the sense of distributions to a 
distribution $G$, which satisfies
\begin{equation}\label{14}
\begin{gathered}
-\Delta G=\delta(x)\delta(y)\delta'(z),\\
G|_{x=-L_1}=G|_{x=L_1}=G|_{y=-L_2}=G|_{y=L_2}=G|_{z=\pm k}=0,~ k \in \mathbb{Z}.
\end{gathered}
\end{equation}

\begin{lemma} \label{FunctionGsign}
The distributions $G$ and $G_\varepsilon$ are positive for $z<0$ and negative 
for $z>0$.
\end{lemma}

\begin{proof}
It is easy to see that 
$\operatorname{sign}(\delta_\varepsilon(x)\delta_\varepsilon(y)\delta_\varepsilon'(z))
=-\operatorname{sign}(z)$, so $\Delta G_\varepsilon>0$ for $z>0$ and 
$\Delta G_\varepsilon<0$ for $z<0$. Now we use the maximum principle, and we 
obtain that $G_\varepsilon$ is positive for $z<0$ and negative for $z>0$. 
Since $G$ is the limit of $G_\varepsilon$ as $\varepsilon \to 0$, 
we conclude that the same result applies for $G$. 
\end{proof}

Now, we  extend the system of equations \eqref{darcy1}, \eqref{darcy2},
 \eqref{darcy3} on the set $[-L_1,L_1]\times[-L_2,L_2]\times \mathbb{R}$ 
in the following manner (see the two-dimensional projection on Figure 1):

(a) we map the area $\Pi$ symmetrically onto 
$[-L_1,L_1]\times[-L_2,L_2]\times [-1,0]$ by introducing 
$\tilde{\rho}(x,y,-z)=\tilde{\rho}(x,y,z)$, $\tilde{P}(x,y,-z)=\tilde{P}(x,y,z)$, 
$g(-z)=-g$, $(x,y,z) \in \Pi$;

(b) we periodically extend this to $[-L_1,L_1]\times[-L_2,L_2]\times \mathbb{R}$.

\begin{figure}[ht]
\label{Expansion}
\begin{center}
 \includegraphics[width=0.5\textwidth]{fig1}
\end{center}
\caption{Extending the system from the area 
$\Pi=[-L_1,L_1]\times[-L_2,L_2]\times [0,1]$ onto 
$[-L_1,L_1]\times[-L_2,L_2]\times \mathbb{R}$ }
\end{figure}


We prove the following theorem.


\begin{theorem} \label{thm2}
Assume that the function $\phi^0$ has a minimum at the point 
$(x_0,y_0)$ and that $\phi^0(0,0)=z_0\in (0,1)$. Then the vertical 
component of the velocity of the fluid $w(0,0,z_0)$ at $(0,0,z_0)$ 
is less than zero (i.e. the fluid moves downward at the tip).
\end{theorem}

\begin{proof}
According to the choice of the approximation of the Heaviside function,
 we have at the initial moment $\tilde{\rho}_\varepsilon(0,0,z_0)\big|_{t=0}=0$. 
Therefore, from \eqref{darcyinterface} and \eqref{darcy1w} we obtain 
the following (at $t=0$ which we imply below):
\begin{equation} \label{15}
\begin{aligned}
\dot{z}(0,0,z_0)
&=w_\varepsilon(0,0,z_0) \\
&=-\frac{\partial \tilde{P}_\varepsilon}{\partial z}(0,0,z_0)
 +\tilde{\rho}_\varepsilon(0,0,z_0)\\
&=-\langle \frac{\partial \tilde{P}_\varepsilon}{\partial z},
 \delta(x) \delta(y) \delta(z-z_0) \rangle \\
&=\langle \tilde{P}_\varepsilon,
 \delta(x) \delta(y) \delta'(z-z_0) \rangle \\
&=\langle \tilde{P}_\varepsilon, \Delta G \rangle
 =\langle \Delta \tilde{P}_\varepsilon, G \rangle
=\langle g\frac{\partial \tilde{\rho}_\varepsilon}{\partial z}, G \rangle
\\
&\to ({\rho}_h-{\rho}_l)g\int_{-L_1}^{L_1}\int_{-L_2}^{L_2}
 G(x,y,\phi^0(x,y))\,dx\,dy \quad \text{as }  \varepsilon\to 0,
\end{aligned}
\end{equation}
 where, when integrating by parts, we took into account the periodicity of
$\tilde{P}$ and $G$ with respect to $z$. Here, $\langle \cdot ,\cdot \rangle$
represents the pairing of a distribution and a test function.

Since $(0,0)$ is the minimum of the function $\phi$, we can conclude that 
$\phi(x,y)>0$ for all $(x,y) \in [-L_1,L_1]\times [-L_2,L_2]$. 
Since $G < 0$ for $z=\phi(x,y)>0$, the whole integral \eqref{15} must be negative. 
Thus, we conclude that $w(0,0,z_0)=\lim\limits_{\varepsilon\to 0}w_\varepsilon(0,0,z_0)$ 
must be negative i.e. the tip of the interface will move downwards.
\end{proof}
 

\begin{thebibliography}{99}

\bibitem{DO} V. G. Danilov, G. A. Omelyanov;
 {\em Dynamics of the interface between two immiscible liquids with nearly 
equal densities under gravity.} Eur. J. Appl. Math., \textbf{13} (2002), 497--516.

\bibitem{KMN} H. Kalisch, D. Mitrovic, J. M. Nordbotten;
{\em Rayleigh-Taylor instability of immiscible fluids in porous media}, 
Continuum Mech. Thermodyn., \textbf{28} (2016), 721--731.

\bibitem{MMN} M.~Marohnic, D.~Mitrovic, A.~Novak;
 {\em On a front evolution in porous media with a source -- analysis and numerics}, 
Bull. Braz. Math. Soc., \textbf{47} (2016), 521--532.

\bibitem{MNK} D.~Mitrovic, J.~M.~Nordboten, H.~Kalisch;
 {\em Dynamics of the interface between immiscible liquids of
different densities with low Froude number}, Nonlinear Anal. Real
World Appl., \textbf{15} (2014), 361--366.

\bibitem{OL}  M.~Orlic, M.~Lazar;
 {\em Cyclonic versus Anticyclonic Circulation in Lakes and Inland Seas}, 
Journal of Physical Oceanography, \textbf{39} (2009), 2247--2263

\bibitem{Tch} A.~Riaz, M.~Hesse, H.~A.~Tchelepi, F.~M.~Orr;
 {\em Onset of convection in a gravitationally unstable diffusive boundary 
layer in porous media}, J. Fluid Mech., \textbf{548} (2006), 87--111.

\bibitem{Seth} J.~A.~Sethian, P.~Smereka;
 {\em Level set methods for fluid interfaces}, Annu. Rev. Fluid Mech., 
\textbf{35} (2003), 341--372.

\bibitem{TNL} A.~Tartakovsky, S.~Neuman, R.~Lenhard;
 {\em Immiscible front evolution in randomly heterogeneous porous media}, Phys.
Fluids, \textbf{15} (2003), 3331--3341 .

\end{thebibliography}

\end{document}