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

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

\begin{document}
\title[\hfilneg EJDE-2015/264\hfil A model of miscible liquids]
{A model of miscible liquids in porous media}

\author[K. Allali, V. Volpert, V. Vougalter \hfil EJDE-2015/264\hfilneg]
{Karam Allali, Vitaly Volpert, Vitali Vougalter}

\address{Karam Allali \newline
Department of Mathematics, Faculty of Sciences and
Technologies, University Hassan II,
PO Box 146, Mohammedia, Morocco}
\email{allali@fstm.ac.ma}

\address{Vitaly Volpert \newline
Institute Camille Jordan, University Lyon 1,
UMR 5208 CNRS,
69622 Villeurbanne,  France}
\email{volpert@math.univ-lyon1.fr}

\address{Vitali Vougalter \newline
University of Toronto, Department of Mathematics,
Toronto, ON, M5S 2E4, Canada}
\email{vitali@math.toronto.edu}

\thanks{Submitted June 3, 2015. Published October 12, 2015.}
\subjclass[2010]{35A01, 35A02, 76D03, 76S05}
\keywords{Darcy approximation; Korteweg stress; miscible liquids; porous media}

\begin{abstract}
 In this article we study the interaction of two  miscible liquids in
 porous media. The model consists of hydrodynamic equations with the
 Korteweg stress terms coupled with the reaction-diffusion equation for the
 concentration. We assume that the fluid is incompressible and its
 motion is described by the Darcy law. The global existence and
 uniqueness of solutions is established for some optimal conditions
 on the reaction source term and external force functions.
 Numerical simulations are performed to show the behavior of two miscible
 liquids subjected to the Korteweg stress.
\end{abstract}

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

\section{Introduction}

There exists transient interfacial phenomena between two miscible liquids
similar to interfacial tension \cite{Jos}.
However they are rather weak and they decay in time because of the mixing
of the two liquids because of the molecular diffusion \cite{Jos,PetMax}.
Investigation of such phenomena is motivated by enhanced oil
recovery, hydrology, frontal polymerization, groundwater pollution
and filtration \cite{Bes1,Hut,Off,Sch,Sch1}.

In $1901$ Korteweg \cite{Kor} introduced additional stress terms in the Navier-Stokes
equations to describe the influence of the composition gradients on
fluid motion.
In $1949$, Zeldovich \cite{Zel} studied the existence of a transitional interfacial
tension and described it with the  expression
$$
\sigma = k \frac{[C]^2}{\delta},
$$
where $[C]$ is the variation of mass fraction through the
transition zone, and $\delta$ is the width of this zone. This
relationship was generalized by Rousar and Nauman \cite{RN} to  systems
far from equilibrium, for linear concentration gradients.
In $1958$, Cahn and Hilliard \cite{CH} introduced the free energy
density for a non-homogeneous fluid,
$$
e = e_0 + k |\nabla \rho|^2,
$$
where $e_0$ is the energy density of a homogeneous fluid and
$\rho$ denotes the density of the fluid.


A miscible liquid model with fully incompressible
Navier-Stokes equations is studied in \cite{Kos}.
Modelling and experiments of miscible
liquids in relation with microgravity experiments were carried out
in \cite{Bes1,Bes2,Bes4,Poj}.
The existence and uniqueness of solutions for miscible liquids model in porous
media is studied in \cite{Kar}.

In this article, we continue the studies of miscible liquids in
porous media. We consider a three-dimensional formulation and
introduce the source terms in the equation of motion and in the
equation for the concentration. The paper is organized as follows.
The next section is devoted to the model presentation, while
Section 3 deals with the existence of  solutions. We establish the
uniqueness of solutions in Section 4 followed in Section 5 by
numerical simulations.


\section{Model presentation}

The model describing the interaction of two miscible liquids is
written as follows:
\begin{gather}\label{prob1}
\frac{\partial C}{\partial t} + u\cdot\nabla C = d \Delta C - C g, \\
\label{u}
\frac{\partial u}{\partial t} + \frac{\mu}{K}u = - \nabla p +
\nabla \cdot T(C)  + f,\\
\operatorname{div}(u) = 0 .
\end{gather}
We consider the boundary conditions:
\begin{equation}
\frac{\partial C}{\partial n} =0, \quad u.n=0, \quad\text{on } \Gamma,
\end{equation}
and the initial conditions:
\begin{equation}\label{prob2}
C(x,0) = C_0(x),\quad u(x,0) = u_0(x),\quad x \in \Omega.
\end{equation}
 Here $u$ is the velocity, $p$ is the pressure, $C$ is
the concentration, $d$ is the coefficient of mass diffusion, $\mu$
is the viscosity, $K$ is the permeability of the medium, $\Gamma$
is a Lipschitz continuous boundary of the open bounded domain
$\Omega$, $n$ is the unit outward normal vector to $\Gamma$, $f$
is the function describing the external forces such as gravity and
buoyancy while the term, $g$ stands for the reaction source term.
The stress tensor terms are given by the relations:
\begin{gather*}
 T_{11} = k \Big(\frac{\partial C}{\partial x_2}\Big)^2, \quad
 T_{12} = T_{21} = -k \frac{\partial C}{\partial x_1}\frac{\partial C}{\partial x_2},
 \quad
 T_{13} = T_{31} = -k \frac{\partial C}{\partial x_1}\frac{\partial C}{\partial x_3},
\\
 T_{23} = T_{32} = -k \frac{\partial C}{\partial
x_2}\frac{\partial C}{\partial x_3},\quad
 T_{22} = k \Big(\frac{\partial C}{\partial x_1}\Big)^2, \ T_{33} = k
\Big(\frac{\partial C}{\partial x_3}\Big)^2,
\end{gather*}
where $k$ is nonnegative constant. We set
\begin{equation}
 \nabla \cdot T(C)
= \begin{pmatrix}
 \frac{\partial T_{11}}{\partial x_1}+\frac{\partial
T_{12}}{\partial x_2}+\frac{\partial T_{13}}{\partial x_3}\\
 \frac{\partial T_{21}}{\partial x_1}+\frac{\partial
T_{22}}{\partial x_2}+\frac{\partial T_{23}}{\partial x_3}\\
 \frac{\partial T_{31}}{\partial x_1}+\frac{\partial
T_{32}}{\partial x_2}+\frac{\partial T_{33}}{\partial x_3}
\end{pmatrix}.
\end{equation}


To state the problem in a variational form we need to introduce function
spaces:
\begin{gather*}
S_u=\{u\in H({\rm div};\Omega); {\rm div} (u)=0, u.n=0\text{ on } \Gamma\} , \\
S_C=\{C\in H^2(\Omega); \frac{\partial C}{\partial n} =0\text{ on } \Gamma\}.
\end{gather*}
The variational form of the problem is to find  $C$, $u$ such that for
all $B$, $v$ the following equalities hold:
\begin{gather}\label{1}
(\frac{\partial C}{\partial t},B) + d (\nabla C,\nabla B) +
(u.\nabla C,B)+(g C,B)=0, \\
\label{2}
(\frac{\partial u}{\partial t},v)+\mu_p(u,v)-(div\
T(C),v)-(f,v)=0.
\end{gather}
Here $ \mu_p = {\mu}/{K}$. The functions $f(x,t)$ and $g(x)$
are assumed to be a sufficiently regular in $\Omega$ such that the
first is bounded in $L^{\infty}(0,t;L^2(\Omega))$, the second is
bounded in $L^{\infty}(\Omega)$ and both of them are positive.

\section{Existence of global solutions}

We begin the proof of existence of solutions with the following lemmas.

\begin{lemma} \label{lem3.1}
The concentration $C$ is bounded in the $L^{\infty}(0,t;L^2)$
space.
\end{lemma}

\begin{proof}
Choosing $C$ as test function in \eqref{1} and taking into account that $g$ is a
positive function, we get the inequality:
$$
\frac{1}{2} \frac{\partial}{\partial t} (C,C) + d (\nabla C,
\nabla C) + (u.\nabla C,C) \le 0.
$$
Since $u\in S_u$, the last term vanishes. The second term is
positive, so integrating by time we obtain:
$$
\|C(t=s)\|_{L^2}^2 \le \|C_0\|_{L^2}^2 .
$$
From this inequality, it follows that $C$ is bounded in
$L^{\infty}(0,t;L^2)$.
\end{proof}


\begin{lemma} \label{lem3.2}
The concentration $C$ is bounded in $L^{\infty}(0,t;H^1)$ and the
velocity $u$ is bounded in $L^{\infty}(0,t;L^2)$.
\end{lemma}

\begin{proof}
Choosing  $-k \Delta C$ as test function in equation \eqref{1}, we have
$$
(\frac{\partial C}{\partial t},-k \Delta C) + (u.\nabla C,-k
\Delta C) = d(\Delta C,-k \Delta C)+(gC,k \Delta C) .
$$
Next, since the reaction source term $g$ is bounded, we get from
the previous estimate:
$$
\frac{k}{2}\frac{\partial}{\partial t}(\nabla C,\nabla C) +
dk(\Delta C,\Delta C) -k(u.\nabla C,\Delta C) \le  kg_0 (\nabla
C,\nabla C).
$$
Then
\begin{equation}\label{3}
\frac{1}{2}\frac{\partial}{\partial t}(\nabla C,\nabla C) +
d(\Delta C,\Delta C) \le (u.\nabla C,\Delta C) + g_0(\nabla
C,\nabla C).
\end{equation}
Also, by choosing in \eqref{2}, $u$ as test function we obtain
\begin{equation}\label{4}
\frac{1}{2} \frac{\partial}{\partial t} (u,u) + \mu_p (u,u)  -
(\nabla.T(C), u)=(f,u).
\end{equation}
To have an explicit expression of $\nabla.T(C)$, we
calculate its first component:
\begin{equation}\label{star}
\begin{aligned}
 \frac{\partial T_{11}}{\partial x_1}+\frac{\partial
T_{12}}{\partial x_2}+\frac{\partial T_{13}}{\partial x_3}
&=2k\frac{\partial C}{\partial x_2}\frac{\partial^2 C}{\partial
x_1\partial x_2}-k\frac{\partial^2 C}{\partial x_1\partial
x_2}\frac{\partial C}{\partial x_2}-k\frac{\partial C}{\partial
x_1}\frac{\partial^2 C}{\partial x_2^2}\\
&\quad -k\frac{\partial^2 C}{\partial x_1\partial x_3}
 \frac{\partial C}{\partial x_2}
-k\frac{\partial C}{\partial x_1}\frac{\partial^2 C}{\partial x_3^2} .
\end{aligned}
\end{equation}
Hence
$$
\frac{\partial T_{11}}{\partial x_1}+\frac{\partial
T_{12}}{\partial x_2}+\frac{\partial T_{13}}{\partial
x_3}=k\frac{\partial C}{\partial x_1}\frac{\partial^2 C}{\partial
x_1\partial x_2}+k\frac{\partial C}{\partial x_1}\frac{\partial^2
C}{\partial x_1\partial x_3}+k\frac{\partial C}{\partial
x_1}\frac{\partial^2 C}{\partial x_1^2}-k\frac{\partial
C}{\partial x_1}\Delta C .
$$
Then
$$
\frac{\partial T_{11}}{\partial x_1}+\frac{\partial
T_{12}}{\partial x_2}+\frac{\partial T_{13}}{\partial
x_3}=\frac{k}{2}\frac{\partial}{\partial x_1}(\nabla C)^2
-k\frac{\partial C}{\partial x_1}\Delta C.
$$
Following the same steps for the second component, we conclude:
$$
\nabla \cdot T = \frac{k}{2}\nabla(\nabla C)^2 -k\Delta C\nabla C\,.
$$
Replacing this last equality in \eqref{4} and since $u\in S_u$, we
have
\begin{equation}\label{5}
\frac{1}{2} \frac{\partial}{\partial t} (u,u) + \mu_p (u,u) -
k(\Delta C\nabla C, u)=(f,u).
\end{equation}
Adding \eqref{5} to the inequality \eqref{3}, and with the fact $u
\in S_u$ and $C \in S_c$, we have
\begin{gather*}
\frac{1}{2}\frac{\partial}{\partial t} \left((u,u)+k(\nabla
C,\nabla C)\right)+\mu_p (u,u)+dk(\Delta C,\Delta C)\le (f,u)+
g_0(\nabla C,\nabla C).
\\
\begin{aligned}
&\frac{1}{2}\frac{\partial}{\partial t} \left((u,u)+k(\nabla
C,\nabla C)\right)+\mu_p (u,u)+dk(\Delta C,\Delta C)\\
&\leq  \frac{1}{2}(f,f)+ \frac{1}{2}(u,u) + g_0(\nabla C,\nabla C).
\end{aligned}
\end{gather*}
Since the third and the fourth terms in the left hand side of this
inequality are positive, we have
$$
\frac{\partial}{\partial t} \left((u,u)+k(\nabla C,\nabla
C)\right) \le (f,f) + (u,u) + 2g_0(\nabla C,\nabla C).
$$
Therefore,
$$
\frac{\partial}{\partial t} \left((u,u)+(\nabla C,\nabla C)\right)
\le \frac{(f,f)}{\min(1;k)} +
\frac{\max(1;2g_0)}{\min(1;k)}\left((u,u) + (\nabla C,\nabla
C)\right).
$$
By integrating over time, and since $f$ is bound in
$L^\infty(0,t;L^2)$, we have
\begin{align*}
&\|u(t=s)\|_{L^2}+\|\nabla C (t=s) \|_{L^2}\\
&\le \|u_0\|_{L^2}+\|\nabla C_0 \|_{L^2} + \frac{f_0}{\min(1;k)}
 + \frac{\max(1;2g_0)}{\min(1;k)}\int_0^t \left(
\|u(s)\|_{L^2}+\|\nabla C (s) \|_{L^2}  \right)ds.
\end{align*}
From the Gronwall's Lemma, it follows that
\begin{align*}
&\|u(t=s)\|_{L^2}+\|\nabla C (t=s) \|_{L^2} \\
&\le (\|u_0\|_{L^2}+\|\nabla C_0 \|_{L^2}
+ \frac{f_0}{\min(1;k)}) \exp\big(\frac{\max(1;2g_0)t}{\min(1;k)}\big).
\end{align*}
We conclude that $C$ is bounded in $L^{\infty}(0,t;H^1)$ and $u$
is bounded in $L^{\infty}(0,t;L^2)$ for $t \in [0;T]$.
\end{proof}


\begin{lemma} \label{lem3.3}
The time derivative of the concentration 
$\frac{\partial C}{\partial t}$ is bounded in $L^2(0,t;L^2)$.
\end{lemma}

\begin{proof}
From  \eqref{1}, since $g$ is a positive function and
by the triangle inequality, we have
$$
\|\frac{\partial C}{\partial t}\|_{L^2} \le d\|\Delta C\|_{L^2} +
\|u\cdot \nabla C\|_{L^2}.
$$
Using H\"{o}lder inequality, we obtain
$$
\|\frac{\partial C}{\partial t}\|_{L^2} \le d\|\Delta C\|_{L^2} +
\|u\|_{L^4}\|\nabla C\|_{L^4},
$$
and from the Gagliardo-Nirenberg inequality, it follows that $\exists
N>0$ such that
$$
\|\frac{\partial C}{\partial t}\|_{L^2} \le d\|\Delta C\|_{L^2} +
N\|u\|_{L^2}^{1/2}\|\nabla u\|_{L^2}^{1/2}\|\nabla
C\|_{L^2}^{1/2}\|\nabla C\|_{H^1}^{1/2}.
$$
We conclude that $ \frac{\partial C}{\partial t}$ is
bounded in $L^2(0,t;L^2)$.
\end{proof}

\begin{lemma} \label{lem3.4}
The time derivative of the velocity $ \frac{\partial
u}{\partial t}$ is bounded in $L^2(0,t;L^2)$.
\end{lemma}

\begin{proof}
To prove this lemma, it is sufficient to remark that $\nabla\cdot T(C)$
is a sum of the expressions of the form
$\lambda D_i(D_j C D_l C)$, where $ D_i=\frac{\partial}{\partial x_i}$,
$i=1,2,3$ and $\lambda$ depending on $i,j$ and $l$ (see for example
\eqref{star}). We have
\begin{align*}
\| D_i(D_j C D_l C)\|_{S'_C}
& \le  \| D_j C D_l C\|_{L^2(\Omega)}\\
& \le  \| D_j C \|_{L^4(\Omega)}\|  D_l C\|_{L^4(\Omega)}\\
& \le  \mathcal{M} \| D_j C \|_{L^2(\Omega)}^{1/2}\|  D_l
C\|_{L^2(\Omega)}^{1/2}\| D_j C \|_{H^1(\Omega)}^{1/2}\|  D_l
C\|_{H^1(\Omega)}^{1/2}.
\end{align*}

We notice that $f$ is a bounded function. Using the the same
reasoning as for the previous lemmas, we prove that $
\frac{\partial u}{\partial t}$ is bounded in $L^2(0,t;L^2)$.
\end{proof}

We can now formulate the main result of this section.

\begin{theorem} \label{thm3.5}
Problem \eqref{prob1}-\eqref{prob2} admits a global solution.
\end{theorem}

\begin{proof}
It is easy to see that the problem admits a finite-dimensional
solutions $C_m$ and $u_m$ defined on the interval of time
$[0;T_m[$. From the previous Lemmas applied to $C_m$ and $u_m$ we
deduce the global existence of those solutions.

Furthermore, the previous Lemmas provide the existence of
subsequences, still denoted by $C_m$ and $u_m$, such that
\begin{gather*}
C_m \to C  \quad \text{weakly in }L^2(0, T; S_C), \\
C_m \to C  \quad\text{weak-star in } L^{\infty}(0, T ;H^1), \\
C'_m \to C' \quad\text{weakly in } L^2(0, T ; S'_C), \\
u_m \to u  \quad \text{weakly in } L^2(0, T ; S_u), \\
u_m \to u  \quad \text{weak-star in } L^{\infty}(0, T ;H_{\rm div}), \\
u'_m \to u' \quad \text{weakly in } L^2(0, T ; S'_u)\,.
\end{gather*}
By classical compactness theorems (see for example
\cite{lio,tem}), we also obtain the strong convergence of  $(C_m;
u_m)$ and by passing to the limit we obtain the existence of
solutions.
\end{proof}

\section{Uniqueness of solution}

To prove uniqueness of solution, we  assume that problem
\eqref{prob1}-\eqref{prob2} has two  solutions
$(C_1,u_1)$ and $(C_2,u_2)$.
From \eqref{prob1}, we have
\begin{equation}\label{6}
\frac{\partial}{\partial t}(C_1-C_2) - d \Delta (C_1-C_2) +
u_1\nabla C_1 - u_2\nabla C_2 + g (C_1-C_2)=0,
\end{equation}
and from \eqref{u}, we also have
\begin{equation} \label{7}
\begin{aligned}
&\frac{\partial}{\partial t}(u_1-u_2) + \mu_p (u_1-u_2)
+ \nabla (p_1-p_2) \\
&= \frac{k}{2} \nabla\left((\nabla
C_1)^2-(\nabla C_2)^2\right)
-k(\Delta C_1 \nabla C_1 - \Delta C_2 \nabla C_2).
\end{aligned}
\end{equation}
Multiplying \eqref{6} by $-k\Delta(C_1-C_2)$ and integrating, we
obtain
\begin{align*}
& (\frac{\partial}{\partial
t}(C_1-C_2),-k\Delta(C_1-C_2)) + dk(\Delta (C_1-C_2),\Delta
(C_1-C_2)) \\
&+ (u_1\nabla C_1,-k\Delta(C_1-C_2))+(u_2\nabla
C_2,k\Delta(C_1-C_2)) \\
&+ ( g (C_1-C_2),-k\Delta(C_1-C_2)) =0.
\end{align*}
Similarly, multiplying \eqref{7} by $u_1-u_2$ and integrating, we
have
\begin{align*}
&(\frac{\partial}{\partial t}(u_1-u_2),u_1-u_2) + \mu_p
(u_1-u_2,u_1-u_2)\\
&= \frac{k}{2}(\nabla\left((\nabla C_1)^2-(\nabla
C_2)^2\right),u_1-u_2) -k(\Delta C_1 \nabla C_1 - \Delta C_2
\nabla C_2,u_1-u_2).
\end{align*}
Adding the two last equalities, using Green's formula and the fact
that $u_i \in S_u$, we conclude that
\begin{align*}
&\frac{1}{2}\frac{\partial}{\partial t}
(\|u_1-u_2\|_{L^2}^2+k\|\nabla C_1- \nabla C_2\|_{L^2}^2) + \mu_p
\|u_1-u_2\|_{L^2}^2+kd \|\Delta(C_1-C_2)\|_{L^2}^2\\
&= k(u_1\nabla(C_1-C_2),\Delta (C_1-C_2))
 +k((u_1-u_2)\nabla C_2,\Delta(C_1-C_2)) \\
&\quad - k(\Delta C_1\nabla(C_1-C_2),u_1-u_2)
 +k(-\Delta C_1 \nabla C_2+\Delta C_2\nabla C_2,u_1-u_2)\\
&\quad +  k( g (C_1-C_2),\Delta(C_1-C_2)) .
\end{align*}
Therefore,
\begin{align}
&\frac{1}{2}\frac{\partial}{\partial t}
(\|u_1-u_2\|_{L^2}^2+k\|\nabla C_1- \nabla C_2\|_{L^2}^2) + \mu_p
\|u_1-u_2\|_{L^2}^2+kd \|\Delta(C_1-C_2)\|_{L^2}^2 \nonumber \\
&=k(u_1\nabla(C_1-C_2),\Delta (C_1-C_2))-k(\Delta
C_1\nabla(C_1-C_2),u_1-u_2) \label{8}\\
&\quad + k( g (C_1-C_2),\Delta(C_1-C_2)). \nonumber
\end{align}
We  now estimate the right-hand side of this equality. We put
$C=C_1-C_2$ and $u=u_1-u_2$. From the H\"{o}lder inequality it follows
that
\[
|(\Delta C_1\nabla C,u)| \le \|\Delta C_1\|_{L^2}\|\nabla C.u\|_{L^2}
\le \|\Delta C_1\|_{L^2}\|\nabla C\|_{L^4}\|u\|_{(L^4)^2}.
\]
Also, from the Gagliardo-Nirenberg inequality we obtain
$$
|(\Delta C_1\nabla C,u)| \le N_1 \|\Delta C_1\|_{L^2}\|\nabla
C\|_{L^2}^{1/2}\|\Delta C\|_{L^2}^{1/2}\|u\|_{L^2}^{1/2}\|\nabla
u\|_{L^2}^{1/2}.
$$
Next, applying the Young's inequality, we obtain
$$
|(\Delta C_1\nabla C,u)| \le \frac{N_1}{4}\|\Delta
C\|_{L^2}^{2}+\frac{3N_1}{4}\|\Delta C_1\|_{L^2}^{4/3}\|\nabla
C\|_{L^2}^{2/3}\|u\|_{L^2}^{2/3}\|\nabla u\|_{L^2}^{2/3}.
$$
Using that same technics, we obtain the  inequality
\begin{align*}
|(u_1\nabla C,\Delta C)| \le \|\Delta C\|_{L^2}\|\nabla
C.u_1\|_{L^2}
 \le \|\Delta C\|_{L^2}\|\nabla
C\|_{L^4}\|u_1\|_{(L^4)^2} .
\end{align*}
Therefore,
$$
|(u_1\nabla C,\Delta C)| \le N_2 \|\Delta C\|_{L^2}^{3/2}\|\nabla
C\|_{L^2}^{1/2}\|u_1\|_{L^2}^{1/2}\|\nabla u_1\|_{L^2}^{1/2} .
$$
Finally,
$$
|(u_1\nabla C,\Delta C)| \le \frac{3N_2}{4} \|\Delta
C\|_{L^2}^{2}+\frac{N_2}{4}\|\nabla
C\|_{L^2}^{2}\|u_1\|_{L^2}^{2}\|\nabla u_1\|_{L^2}^{2}.
$$
From \eqref{8} and assuming that $N_1+3N_2\le 4d$, we have
\begin{align*}
&\frac{1}{2}\frac{\partial}{\partial t} (\|u\|_{L^2}^2+k\|\nabla
C\|_{L^2}^2) \\
&\le \frac{3N_1k}{2}\|\Delta C_1\|_{L^2}^{4/3}\|\nabla
C\|_{L^2}^{2/3}\|u\|_{L^2}^{2/3}\|\nabla u\|_{L^2}^{2/3} \\
&\quad +\frac{N_2k}{2}\|\nabla C\|_{L^2}^{2}\|u_1\|_{L^2}^{2}\|\nabla
u_1\|_{L^2}^{2}+k g_0\|\nabla C\|_{L^2}^{2} \\
& \le (\|u\|_{L^2}^2+k\|\nabla C\|_{L^2}^2)
\Big(\frac{N_2}{2}\|\nabla C\|_{L^2}^{2}\|u_1\|_{L^2}^{2}\|\nabla
u_1\|_{L^2}^{2}\\
&\quad +\frac{3N_1k}{2}\|\Delta C_1\|_{L^2}^{4/3}\|\nabla
C\|_{L^2}^{2/3}\|u\|_{L^2}^{-4/3}\|\nabla u\|_{L^2}^{2/3} +k g_0
\Big).
\end{align*}
Denote
\begin{gather*}
\phi(t)=\|\nabla C\|_{L^2}^{2}\|u_1\|_{L^2}^{2}\|\nabla
u_1\|_{L^2}^{2}+\|\Delta C_1\|_{L^2}^{4/3}\|\nabla
C\|_{L^2}^{2/3}\|u\|_{L^2}^{-4/3}\|\nabla u\|_{L^2}^{2/3}+1, \\
M=\max\big(\frac{N_2}{2},\frac{3N_1k}{2},kg_0\big) .
\end{gather*}
Then we have the  estimate
$$
\frac{d}{dt}(\exp\Big(M\int_0^t \phi(s)ds\Big)(\|u\|_{L^2}^2+k\|\nabla
C\|_{L^2}^2))\le 0.
$$
for all  $t \ge 0$. From this we deduce that
$$
\exp\Big(M\int_0^t \phi(s)ds\Big)(\|u\|_{L^2}^2+k\|\nabla C\|_{L^2}^2)\le
\|u(0)\|_{L^2}^2+k\|\nabla C(0)\|_{L^2}^2.
$$
Since $u(0)=C(0)=0$, we conclude the uniqueness of solution. We
can now state the theorem on the uniqueness of solution.

\begin{theorem} \label{thm4.1}
Problem \eqref{prob1}-\eqref{prob2} admits a unique solution.
\end{theorem}


\section{Numerical simulations}

For numerical simulations, we will consider the 2D problem without reaction
term and external forces. We will introduce the stream function defined
by the equalities
$$
u_1 = \frac{\partial \psi}{\partial x_2}, \quad
u_2 = -\frac{\partial \psi}{\partial x_1},
$$
and the vorticity
$\omega  = rot(u)$.
The problem becomes
\begin{gather}\label{omega}
\frac{\partial \omega}{\partial t} + \mu_p \omega =
\frac{\partial}{\partial x_1}(\frac{\partial T_{21}}{\partial
x_1}+\frac{\partial T_{22}}{\partial
x_2})-\frac{\partial}{\partial x_2}(\frac{\partial
T_{11}}{\partial x_1}+\frac{\partial T_{12}}{\partial x_2}) , \\
\label{C} \frac{\partial C}{\partial t}+(\frac{\partial \psi}{\partial
x_2},-\frac{\partial \psi}{\partial x_1}).\nabla C=d\Delta C , \\
\label{psi}
\omega = -\Delta \psi\,.
\end{gather}


\subsection*{Numerical method}

We begin with equation \eqref{C}. It  is solved by the alternative direction
implicit finite difference method with Thomas algorithm:
\begin{align*}
&\frac{C_{i,j}^{n+\frac{1}{2}}-C_{i,j}^{n}}{h_t/2}\\
&= d \Big(
\frac{C_{i-1,j}^{n+\frac{1}{2}}-2C_{i,j}^{n+\frac{1}{2}}
+C_{i+1,j}^{n+\frac{1}{2}}}{h_x^2}+\frac{C_{i,j-1}^{n}-2C_{i,j}^{n}
+C_{i,j+1}^{n}}{h_y^2} \Big)\\
&\quad
-\frac{\psi_{i,j+1}^{n}-\psi_{i,j-1}^{n}}{2h_y}
 \frac{C_{i+1,j}^{n+\frac{1}{2}}-C_{i-1,j}^{n+\frac{1}{2}}}{2h_x}
+ \frac{\psi_{i+1,j}^{n}-\psi_{i-1,j}^{n}}{2h_x}\frac{C_{i,j+1}^{n}
 -C_{i,j-1}^{n}}{2h_y} \, ,
\end{align*}
\begin{align*}
&\frac{C_{i,j}^{n+1}-C_{i,j}^{n+\frac{1}{2}}}{h_t/2}\\
&= d \Big(
\frac{C_{i-1,j}^{n+\frac{1}{2}}-2C_{i,j}^{n+\frac{1}{2}}
 +C_{i+1,j}^{n+\frac{1}{2}}}{h_x^2}+\frac{C_{i,j-1}^{n+1}-2C_{i,j}^{n+1}
 +C_{i,j+1}^{n+1}}{h_y^2} \Big) \\
&\quad -\frac{\psi_{i,j+1}^{n}-\psi_{i,j-1}^{n}}{2h_y}\frac{C_{i+1,j}^{n+\frac{1}{2}}
-C_{i-1,j}^{n+\frac{1}{2}}}{2h_x}
+ \frac{\psi_{i+1,j}^{n}-\psi_{i-1,j}^{n}}{2h_x}\frac{C_{i,j+1}^{n+1}
 -C_{i,j-1}^{n+1}}{2h_y} \, .
\end{align*}

In  \eqref{omega} we replace $T_{ij}$ by their expressions through
 the concentrations
\begin{equation}
 \frac{\partial \omega}{\partial t} + \mu_p \omega =
k \Big( \frac{\partial C}{\partial x_1}\Big(\frac{\partial^3
C}{\partial x_2^3}+\frac{\partial^3 C}{\partial x_1^2\partial
x_2}\Big)- \frac{\partial C}{\partial x_2}\Big(\frac{\partial^3
C}{\partial x_1^3}+\frac{\partial^3 C}{\partial x_1\partial
x_2^2}\Big)\Big)
\end{equation}
We use the finite difference scheme
\begin{align*}
\omega_{i,j}^{n+1}
&= \frac{1}{1+h_t\mu_p}\omega_{i,j}^{n}+k\frac{h_t}{1+h_t\mu_p}
\bigg(\frac{C_{i+1,j}^{n+1}-C_{i-1,j}^{n+1}}{2h_x}\\
&\quad\times \Big(\frac{C_{i,j+2}^{n+1}-2C_{i,j+1}^{n+1}+2C_{i,j-1}^{n+1}-C_{i,j-2}^{n+1}}{2h_y^3}\\
&\quad +\frac{(C_{i+1,j+1}^{n+1}-C_{i+1,j-1}^{n+1})-2(C_{i,j+1}^{n+1}-C_{i,j-1}^{n+1})+(C_{i-1,j+1}^{n+1}-C_{i-1,j-1}^{n+1})}{2h_x^2
h_y}\Big)\\
&\quad -\frac{C_{i,j+1}^{n+1}-C_{i,j-1}^{n+1}}{2h_y}
\Big(\frac{C_{i+2,j}^{n+1}-2C_{i+1,j}^{n+1}
+2C_{i-1,j}^{n+1}-C_{i-2,j}^{n+1}}{2h_x^3}\\
&\quad + \Big((C_{i+1,j+1}^{n+1}-C_{i-1,j+1}^{n+1})
 -2(C_{i+1,j}^{n+1}-C_{i-1,j}^{n+1})
+(C_{i+1,j-1}^{n+1}\\
&\quad  -C_{i-1,j-1}^{n+1})\Big)\big/
\big(2h_y^2 h_x\big) \Big)\bigg) \,.
\end{align*}
Equation \eqref{psi} is solved by the fast Fourier transform method.



\paragraph{Numerical results}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=6cm,height=5cm]{fig1a} % drop1(d3_10-3k10-7mup10max1T0)}
\includegraphics[width=6cm,height=5cm]{fig1b} % drop2(d3_10-3k10-7mup10max_682T100)}
\end{center}
\caption{Evolution of the concentration during $100$ seconds for
$d = 3\times 10^{-3}$, $k = 10^{-7}$ and $\mu_p = 100$.}
\label{fig1}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=6cm,height=5cm]{fig2a} % max_psi
\includegraphics[width=6cm,height=5cm]{fig2b} % contour
\end{center}
\caption{The maximum of stream function as function of time for $d
= 3\times 10^{-3}$, $k = 10^{-7}$ and $\mu_p = 10$ (left). The
stream lines for the same parameters and after $100$ s(right)}
\label{fig2}
\end{figure}

An example of numerical simulations is shown in Figures \ref{fig1} and
\ref{fig2}.
Figure \ref{fig1} shows the evolution of the miscible drop in time. The
transient interfacial tension affects the geometry of the drop and
its shape becomes more spherical. At the same time, the maximum of
the concentration decreases due to diffusion (Figure \ref{fig2}, left). The
stream lines are shown in Figure \ref{fig2} (right). Though transient
interfacial phenomena are sufficiently weak, they provoke the
motion of fluid which is initially quiescent.

\begin{thebibliography}{99}

\bibitem{Kar} K. Allali;
\emph{Existence and uniqueness of solutions for
miscible liquids model in porous media. Electronic Journal of
Differential Equations}, 2013, vol. 2013, no 254, p. 1-7.

\bibitem{Bes1} N. Bessonov, V. A. Volpert, J. A. Pojman, B. D. Zoltowski;
\emph{Numerical simulations of convection induced by Korteweg stresses
in miscible polymermonomer systems}, Microgravity-Science and
Technology \textbf{17} (2005) 8--12.

\bibitem{Bes2}  N. Bessonov, J. Pojman, G. Viner, V. Volpert, B. Zoltowski
Instabilities of diffuse interfaces. Math. Model. Nat. Phenom., 3
(2008), No. 1, 108-125.

\bibitem{Bes4}   N. Bessonov, J. Pojman, V. Volpert;
\emph{Modelling of diffuse interfaces with temperature gradients}.
 J. Engrg. Math.  49  (2004),  no. 4, 321--338.

\bibitem{CH} J. W. Cahn, J. E. Hilliard;
\emph{Free Energy of a Nonuniform System. I. Interfacial Free Energy,}
J. Chem. Phys. 28 (1958), 258-267.

\bibitem{Hut} C. A. Hutchinson, P. H. Braun;
 \emph{Phase relations of miscible displacement in oil recovery},
 AIChE Journal \textbf{7} (1961), 64--72.

\bibitem{Jos} D. D. Joseph;
\emph{Fluid dynamics of two miscible liquids with
diffusion and gradient stresses}, European journal of mechanics.
B, Fluids \textbf{9} (1990) 565--596.

\bibitem{Kor} D. Korteweg;
\emph{Sur la forme que prennent les équations du
mouvement des fluides si l'on tient compte des forces capillaires
causées par des variations de densité}, Arch. Néerl Sci. Exa.
Nat., Series II. \textbf{6} (1901) 1--24.

\bibitem{Kos} I. Kostin, M. Marion, R. Texier-Picard, V. Volpert;
\emph{Modelling of miscible liquids with the Korteweg stress},
ESAIM: Mathematical Modelling and Numerical Analysis \textbf{37}(2003) 741--753.

\bibitem{lio} J. L. Lions;
\emph{Quelques méthodes de résolution des problèmes aux
limites non linéaires}, Gauthier-Villars, Paris (1969).

\bibitem{Off} J. Offeringa,  C. Van der Poel;
\emph{Displacement of oil from
porous media by miscible liquids}, Journal of Petroleum Technology
\textbf{5} (1954), 37--43.

\bibitem{PetMax} P. Petitjeans, T. Maxworthy;
\emph{Miscible displacements in
capillary tubes. Part 1. Experiments}, Journal of Fluid Mechanics
\textbf{326} (1996) 37--56.

\bibitem{Poj}  J. Pojman, N. Bessonov, V. Volpert, M. S. Paley;
\emph{Miscible Fluids in Microgravity (Mfmg):  A Zero-Upmass
 Investigation on the International Space Station},
 Microgravity Sci. Tech. 2007, XIX,
 Volume 19, issue 1, pages 33-41.

\bibitem{RN} I. Rousar, E. B. Nauman;
\emph{A Continuum Analysis of Surface Tension in Nonequilibrium Systems,}
Chem. Eng. Comm. 129 (1994), 19-28.

\bibitem{Sch} F. Schwille;
\emph{Groundwater pollution in porous media by fluids
immiscible with water}, Science of the Total Environment \textbf{21}
(1981), 173--185.

\bibitem{Sch1} F. Schwille;
\emph{Migration of organic fluids immiscible with
water in the unsaturated zone}, Pollutants in porous media.
Springer Berlin Heidelberg, (1984), 27--48.

\bibitem{tem} R. Temam;
\emph{Navier-Stokes equations. Theory and numerical analysis,}
North-Holland Publishing Co., Amsterdam, New York,
Stud. Math. Appl. 2 (1979).

\bibitem{Zel} Y. B. Zeldovich;
\emph{About Surface Tension of a Boundary between
two Mutually Soluble Liquids}, Zhur. Fiz. Khim. 23 (in Russian)
(1949), 931-935.

\end{thebibliography}

\end{document}