Global existence for a two-phase flow model with cross diffusion

TL;DR

This paper proves the global existence of weak solutions for a complex two-phase flow model with cross diffusion, derived from thermodynamic principles, using entropy methods and a priori bounds.

Contribution

It introduces a novel analysis for a degenerate pseudo-parabolic system with cross diffusion, establishing global solutions in bounded domains.

Findings

Global existence of weak solutions proven

Entropy inequality used for control of degeneracy

Model derived from thermodynamic principles

Abstract

In this work we study a degenerate pseudo-parabolic system with cross diffusion describing the evolution of the densities of an unsaturated two-phase flow mixture with dynamic capillary pressure in porous medium with saturation-dependent relaxation parameter and hypocoercive diffusion operator modeling cross diffusion. The equations are derived in a thermodynamically correct way from mass conservation laws. Global-in-time existence of weak solutions to the system in a bounded domain with equilibrium boundary conditions is shown. The main tools of the analysis are an entropy inequality and a crucial apriori bound which allows for controlling the degeneracy.

Equations

$${\partial }_t{S}_i=$$

$$i=1,\dots ,n,x\in \Omega ,t>0.$$

$$\displaystyle\mu_{i}=\log\Big{(}\frac{S_{i}}{S}\Big{)}+\int_{1/2}^{S}\frac{\tau(\sigma)}{a(\sigma)}d\sigma\qquad i=1,\ldots,n.$$

$$S_i$$

$$S_i(\cdot ,0)$$

$$\mathcal{D}=\left\{\mathbf{S}\in {\mathbb{R}}^n:S_i>0\text{ for }i=1,\dots ,n,\sum _{j=1}^n{S}_j<1.\right\}.$$

$${\mu }_i$$

$$\mathcal{F}(\mathbf{S})$$

$$p^{th}(\mathbf{S})$$

$$\sum _i{S}_i\nabla {\mu }_i=\nabla p^{th}=\frac{S\nabla \beta (S)}{a(S)}=\frac{S\tau (S)}{a(S)}\nabla S=\nabla {\int }_{1/2}^S\frac{\sigma \tau (\sigma )}{a(\sigma )}d\sigma .$$

$$\frac{d}{dt}{\int }_{\Omega }\left(\tilde{\mathcal{F}}(\mathbf{S})+\frac{1}{2}|\nabla \beta (S){|}^2\right)dx$$

$$=-{\int }_{\Omega }{\beta }^{\prime }(S)p_c^{\prime }(S)|\nabla S{|}^{2}dx-{\int }_{\Omega }\sum _{i,j=1}^n{D}_{ij}(\mathbf{S})\nabla {\mu }_i\cdot \nabla {\mu }_{j}dx\le 0,$$

$$\tilde{\mathcal{F}}(\mathbf{S})=\mathcal{F}(\mathbf{S})-\mathcal{F}({\mathbf{S}}^{\Gamma })-\mathbf{\mu }({\mathbf{S}}^{\Gamma })\cdot (\mathbf{S}-{\mathbf{S}}^{\Gamma }).$$

$$\exists C\in \mathbb{R}:\sum _{i=1}^n{S}_i\frac{\partial \mathcal{F}}{\partial S_i}(\mathbf{S})-\mathcal{F}(\mathbf{S})=C+{\int }_{1/2}^S\frac{\sigma \tau (\sigma )}{a(\sigma )}d\sigma .$$

$$\sum _{i=1}^n{S}_i\frac{\partial {\mathcal{F}}_0}{\partial S_i}(\mathbf{S})-{\mathcal{F}}_0(\mathbf{S})=0.$$

$${\mu }_i^{\ast }={\mu }_i-\frac{1}{n}\sum _{j=1}^n{\mu }_j,i=1,\dots ,n.$$

$$S_i=S\frac{e^{{\mu }_i^{\ast }}}{{\sum }_{j=1}^n{e}^{{\mu }_j^{\ast }}},i=1,\dots ,n.$$

$$D_0|\Pi v{|}^2\le \sum _{i,j=1}^n{D}_{ij}(\mathbf{S})v_i{v}_j\le D_1|\Pi v{|}^2\text{ for all }v\in {\mathbb{R}}^n,\mathbf{S}\in \mathcal{D},$$

$$a(S)=\frac{{\lambda }_o(S){\lambda }_w(S)}{{\lambda }_o(S)+{\lambda }_w(S)}=\frac{(1-S{)}^{\lambda }{S}^{\gamma }}{(1-S{)}^{\lambda }+S^{\gamma }},$$

$$p_c^{\prime }(S)=\frac{1}{S^{{\beta }_1}}+\frac{1}{(1-S{)}^{{\beta }_2}},$$

$$\displaystyle\tau(S)=\beta^{\prime}(S)=\frac{S^{\gamma}}{S^{\gamma}+(1-S)^{\lambda}}\Big{[}1+\frac{(1-S)^{\lambda}}{S^{\gamma_{1}}}\Big{]},$$

$$\displaystyle 5<\beta_{1}\leq\gamma_{1}<\gamma<\frac{1}{2}\beta_{1}+\frac{5}{6}(\gamma_{1}-2\big{)},\quad 5<\beta_{2}\leq\lambda<3\beta_{2}-10.$$

$$\sum _{i,j=1}^n{D}_{ij}(\mathbf{S})v_i{v}_j\le D_1|\Pi v{|}^2\text{mbox}forv\in {\mathbb{R}}^n,\mathbf{S}\in \mathcal{D},$$

$$D_0|\Pi v{|}^2\le \sum _{i,j=1}^n{D}_{ij}(\mathbf{S})v_i{v}_j\text{mbox}forv\in {\mathbb{R}}^n,\mathbf{S}\in \mathcal{D},$$

$$\nabla \beta (S),{\partial }_t\beta (S),\sqrt{a(S)p_c^{\prime }(S)},\sqrt{a(S)}\nabla {\partial }_t\beta (S)\in L_{loc}^2(0,\infty ;L^2(\Omega )),$$

$$\text{mbox}fori=1,\dots ,n:(\Pi \mu {)}_i\in L_{loc}^2(0,\infty ;H^1(\Omega )),$$

$$\text{mbox}fori=1,\dots ,n:{\partial }_t{S}_i\in L_{loc}^2(0,\infty ;H^{-1}(\Omega )),$$

$$\sum _{i=1}^n{\int }_0^T⟨{\partial }_t{S}_i,{\varphi }_i⟩dt$$

$$+{\int }_0^T{\int }_{\Omega }\sum _{i,j=1}^n{D}_{ij}(\mathbf{S})\nabla {\mu }_j\cdot \nabla {\varphi }_{i}dxdt=0$$

$$\forall {\varphi }_1,\dots ,{\varphi }_n\in L_{loc}^2(0,\infty ;H_0^1(\Omega )),$$

$$\displaystyle S_{i}(\cdot,t)\to S_{i}^{0}\qquad\mbox{strongly in $H^{-1}(\Omega)$ as $t\to 0$.}$$

$$\underset{i=1,\dots ,n}{\min }\underset{\Omega }{\inf }{S}_i^0>0,\underset{\Omega }{\max }{S}^0<1,\beta (S^0)\in H^1(\Omega ).$$

$$\displaystyle\boldsymbol{w}=\boldsymbol{\mu}(\boldsymbol{S})-\boldsymbol{\mu}(\boldsymbol{S}^{\Gamma})+\Big{(}\frac{\beta(S)-\beta(S^{k-1})}{\kappa}\Big{)}\mathbf{1},$$

$$[{\mathcal{F}}^{\prime \prime }{]}_{ij}=\frac{{\partial }^2\mathcal{F}}{\partial S_j\partial S_i}=\{\begin{array}{ll}\frac{1}{1-S},& i\ne j,\\ \frac{1}{S_i}+\frac{1}{(1-S)},& i=j,\end{array}$$

Full text

The first and the third author acknowledge partial support from the Austrian Science Fund (FWF), grants P22108, P24304, W1245, P27352 and P30000. All three authors were partially supported by the bilaterial project No. HR 04/2018 of the Austrian Exchange Sevice OeAD together with the Ministry of Science and Education of the Republic of Croatia MZO

1. Introduction

The problem of describing the transport of chemical mixtures in porous media is very important in many industrial applications. For a general overview on the modeling of multicomponent multiphase flows in porous media, we refer to [2]. In this paper we consider a two-phase flow model with wetting and non-wetting phase (e.g. water and oil), where the non-wetting phase consists of a mixture of $n$ chemical components, including nonequilibrium effects concerning capillary pressure and cross-diffusion effects. The main result of this work is to provide an existence analysis of the proposed model. From a mathematical viewpoint, the transport equations for the mass densities form a degenerate pseudo-parabolic system of PDEs with cross-diffusion terms. The presence of the mixed-derivative third-order term, coming from the nonequilibrium capillary pressure law, in form of a time derivative inside the diffusion operator, as well as the cross-diffusion terms, involving the chemical potentials, make the analysis very demanding. Furthermore, the compactness of an approximate regularized system is obtained by applying the nonstandard compactness results of Dreyer et al. [5].

The modeling of nonequilibrium capillary effects in problems of enhancing oil and gas recovery from rocks was proposed by Barenblatt, Entov and Ryzhik in the classical book [1], and later investigated by many scientists up to nowadays. In our work we follow the approach given by Hassanizadeh and Grey [10], where the nonequilibrium capillary effects are given by a constitutive relationship between the non-wetting phase saturation and the capillary pressure. This relationship is characterized by the presence of the relaxation parameter which depends on the non-wetting saturation as well.

Concerning the mathematical analysis, the global-in-time existence of weak solutions for the Richards’ equation with dynamic capillary pressure and constant relaxation parameter was shown by Mikelić [14]. The first existence result for the two-phase flow model with dynamic capillary pressure and saturation dependent relaxation parameter was obtained by Cao and Pop in [7]. We note that the existence theorem can be proved under certain relations between the orders of the zeros of the relative permeabilities and the relaxation parameter and the order of the singularities of the capillary pressure function. In comparison to [7], here we follow the approach given in [15], where it was shown that it is enough to analyze the case of the countercurrent imbibition flow instead of the full two-phase flow system.

On the other side, the analysis of a model describing the transport of a single-phase fluid mixture in porous media taking into account also certain cross-diffusion effects was studied in [13]. The equations are derived in a thermodynamically consistent way, and global-in-time existence of weak solutions in a bounded domain with equilibrium boundary conditions as well as long-time behaviour was proved with the help of the boundedness-by-entropy method [3, 11, 12]. The mathematical novelties rely on the complex structure of the equations and on the observation that the solution of the binary model satisfies an unexpected integral inequality leading to a minimum principle for this system.

Our goal in this work is to combine the strategies of [15] and [13], leading to a global-in-time existence of weak solutions result for a two-phase flow model with cross diffusion.

Finally, up to our knowledge, the uniqueness and the long-time behaviour of a weak solution for a two-phase flow model with saturation-dependent relaxation parameter and cross diffusion are still open problems. For a uniqueness result of a two-phase flow model with saturation-dependent relaxation parameter but without cross diffusion, we mention the result in [6].

2. Model equations

We consider an incompressible, isothermal fluid mixture with $n$ components in a domain $\Omega\subset{\mathbb{R}}^{3}$. We note that the fact that we work in ${\mathbb{R}}^{3}$ is for convenience only, and can be easily adapted to an arbitrary space dimension $\Omega\subset{\mathbb{R}}^{d}$ with $d\geq 1$. The evolution of this fluid mixture is governed by the transport equations for the single component mass densities $S_{1}(x,t),\ldots,S_{n}(x,t)$ in the following way

[TABLE]

Here $S=\sum_{i=1}^{n}S_{i}$ is the total mass density, $\boldsymbol{S}=(S_{1},\ldots,S_{n})$ is the vector of the single component mass densities, $a(S)$ is the diffusion mobility, $p_{c}(S)$ represents the stationary capillary pressure, $\tau(S)\equiv\beta^{\prime}(S)$ plays the role of a relaxation parameter, $D=(D_{ij}(\boldsymbol{S}))_{i,j=1,\ldots,n}$ is the diffusion matrix, and the quantities $\mu_{1},\ldots,\mu_{n}$, called chemical potentials, are defined in terms of $S_{1},\ldots,S_{n}$ as follows

[TABLE]

The sum $p_{c}^{\textrm{dyn}}(S)\equiv p_{c}(S)+\partial_{t}\beta(S)$ is referred to as dynamic capillary pressure [10]. The quantities $a(S)$, $\tau(S)$, $p_{c}^{\prime}(S)$ are assumed to be positive for $0<S<1$, while the diffusion matrix $D(\boldsymbol{S})$ is assumed to be positive semidefinite.

Following the approach in [13], we impose equilibrium boundary conditions

[TABLE]

where $S_{1}^{\Gamma},\ldots,S_{n}^{\Gamma}>0$ are generic constants, as well as general initial conditions

[TABLE]

For consistency of (1), (2) with the physics, we require the single component concentrations $S_{1},\ldots,S_{n}$ to be positive and the total concentration $S$ to be smaller than 1; that is, we seek for solutions $\boldsymbol{S}$ to (1), (2) which take values in the set

[TABLE]

The chemical potentials $\mu_{1},\ldots,\mu_{n}$ are the partial derivatives with respect to the species concentrations $S_{1},\ldots,S_{n}$ of a free energy density function ${\mathcal{F}}$ satisfying

[TABLE]

The thermodynamic pressure $p^{th}$ is given by the Gibbs-Duhem equation

[TABLE]

The gradient of the thermodynamic pressure $p^{th}$ satisfies the simple relation

[TABLE]

As a consequence of (7), by employing $\mu_{i}-\mu_{i}(\boldsymbol{S}^{\Gamma})$ as a test function in (1), one obtains the following entropy balance equation:

[TABLE]

where the relative entropy density $\tilde{\mathcal{F}}$ is defined as

[TABLE]

Remark 1**.**

Relations (6), (7) easily imply

[TABLE]

Equation (10) constitutes a necessary condition in order for the entropy balance equation (8) to hold; without (10) it is unclear how to handle the contribution of the nonstationary term $\partial_{t}\beta(S)$ in the dynamic capillary pressure $p^{dyn}$. In other words, (10) is a constraint on the possible choices of free energies ${\mathcal{F}}$ which ensure that (1) possesses an entropy structure.**

Since (10) is a linear nonhomogeneous equation, we can write any solution ${\mathcal{F}}$ to (10) as ${\mathcal{F}}={\mathcal{F}}_{0}+{\mathcal{F}}_{1}$, where ${\mathcal{F}}_{1}$ is a specific solution to (10), while ${\mathcal{F}}_{0}$ is a generic solution to the corresponding linear homogeneous equation:

[TABLE]

A simple ansatz ${\mathcal{F}}_{1}(\boldsymbol{S})=\tilde{\mathcal{F}}_{1}(S)$ yields ${\mathcal{F}}(\boldsymbol{S})=\mathcal{E}(S)$ (up to additive constants). On the other hand, Euler’s theorem on homogeneous functions implies that (11) is equivalent to the condition that ${\mathcal{F}}_{0}$ should be homogeneous of degree 1, i.e. ${\mathcal{F}}_{0}(\lambda\boldsymbol{S})=\lambda{\mathcal{F}}_{0}(\boldsymbol{S})$ for every $\boldsymbol{S}\in{\mathcal{D}}$, $\lambda>0$. This condition has to be put together with the requirement that ${\mathcal{F}}$ has to be convex and the mapping $\boldsymbol{S}\in{\mathcal{D}}\mapsto(\mu_{1},\ldots,\mu_{n})\in{\mathbb{R}}^{n}$ globally invertible. A natural choice of ${\mathcal{F}}_{0}$ which fulfills all these requirements is ${\mathcal{F}}_{0}(\boldsymbol{S})=\sum_{i=1}^{n}S_{i}\log(S_{i}/S)$.**

Other quantities that will play a role in the analysis of (1) are the relative chemical potentials:

[TABLE]

The concentrations $S_{1},\ldots,S_{n}$ can be easily written in terms of the total concentration and the relative chemical potentials:

[TABLE]

The structure of the paper is as follows. In Section 3 the main result of the paper is stated and the state of the art for systems of the form (1) is described. In Section 4 some auxiliary results are stated and proved. In Section 5 Theorem 4 is proved. In the Appendix the derivation of the model is shown.

3. Main result

Throughout the paper we make the following assumptions:

1. (H1)

The diffusion matrix $D=(D_{ij}(\boldsymbol{S}))_{i,j=1}^{n}$ is symmetric and positive semidefinite (Onsager’s principle of thermodynamics). Moreover, constants $D_{0}$, $D_{1}>0$ exist such that

[TABLE]

where $\Pi=I-l\otimes l$ is the orthogonal projection on the subspace of ${\mathbb{R}}^{n}$ orthogonal to $l=(1,\ldots,1)/\sqrt{n}$. 2. (H2)

The diffusion mobility $a(S)$ is given by

[TABLE]

for some constants $\lambda,\gamma>0$. 3. (H3)

The stationary capillary pressure $p_{c}(S)$ has the form

[TABLE]

for some constants $\beta_{1},\beta_{2}>0$. 4. (H4)

We assume that the relaxation parameter $\tau(S)$ is given by

[TABLE]

for some constant $\gamma_{1}>0$. 5. (H5)

The following algebraic relations are satisfied:

[TABLE]

Remark 2**.**

In order to avoid technical difficulties, we use explicit forms for $a$, $p_{c}$ and $\tau$ like in [15].

Remark 3**.**

We point out that the upper bound

[TABLE]

is consistent with the fact that the diffusion fluxes $\mathbf{J}_{i}=-\sum_{j=1}^{n}D_{ij}(\boldsymbol{S})\nabla\mu_{j}$ ($i=1,\ldots,n$) sum up to zero: $\sum_{i=1}^{N}\mathbf{J}_{i}=0$. On the other hand, the lower bound

[TABLE]

often referred to as hypocoercivity, is the strongest coercivity property that $D$ can satisfy under the constraint $\sum_{i=1}^{N}\mathbf{J}_{i}=0$. As a consequence of this assumption, the diffusion fluxes $\mathbf{J}_{i}=-\sum_{j=1}^{n}D_{ij}(\boldsymbol{S})\nabla\mu_{j}$ only depend on the gradients of the relative chemical potentials: $\mathbf{J}_{i}=-\sum_{j=1}^{n}D_{ij}(\boldsymbol{S})\nabla\mu_{j}^{*}$.**

We now present our definition of weak solution to (1)–(4). In the following, the symbol $\langle\cdot,\cdot\rangle$ represents the duality product between $H^{-1}(\Omega)$ and $H^{1}_{0}(\Omega)$.

Definition 1** (Weak solution).**

A function $\boldsymbol{S}:\Omega\times(0,\infty)\to\mathcal{D}$ is called a global-in-time weak solution to (1)–(4) if and only if the following properties are fulfilled:

[TABLE]

as well as the weak formulation of (1):

[TABLE]

relation (2), the boundary conditions (3)111We point out that if $\beta(S)$ and $\mu_{i}^{*}=(\Pi\mu)_{i}$ belong to $L^{2}_{loc}(0,\infty;H^{1}(\Omega))$ for $i=1,\ldots,b$, then they admit trace on $\partial\Omega$, therefore also $S_{1},\ldots,S_{n}$ admit trace on $\partial\Omega$ thanks to the invertibility of $S\mapsto\beta(S)$ and relation (12)., and the initial condition (4):

[TABLE]

The result we present in this paper is concerned with the global existence of weak solutions to (1)–(4).

Theorem 4** (Existence of global weak solutions).**

Let $S_{1}^{0},\ldots,S_{n}^{0}:\Omega\to{\mathbb{R}}_{+}$ be measurable functions satisfying

[TABLE]

Assume that Assumptions **(H1)–(H5)*** hold. Then there exists a global-in-time weak solution $\boldsymbol{S}:\Omega\times(0,\infty)\to\mathcal{D}$ to (1)–(4).*

Key idea of the proof

The proof of Thr. 4 is based on the entropy method [3, 11, 12]. The starting point of the argument is the formulation of a time-discretized and regularized version of (1). Such approximate equation is stated in terms of the variables $w_{i}=\mu_{i}+\partial_{t}\beta(S)$, $i=1,\ldots,n$ (or rather a discretized version of it). One of the key ingredients of the proof is the entropy balance equation (8), which yields crucial gradient estimates. The other key tool employed in the proof is a result shown in [5], which allows to prove compactness for the densities $S_{1},\ldots,S_{n}$ if some bounds for the gradient of the relative chemical potentials $\nabla\mu_{1}^{*},\ldots,\nabla\mu_{n}^{*}$ are known, together with compactness of the total density $S$. We point out that in the standard entropy method the approximate problem is formulated in terms of the “entropy variables” defined as partial derivatives of the mathematical entropy (or energy) density, which in the case here considered would be the functions $\mu_{1},\ldots,\mu_{n}$ given by (2). However, this standard approach does not work in this setting: in fact, in order to obtain a crucial estimate for the dynamic capillary pressure, $\partial_{t}\beta(S)$ must be used as a test function in the weak formulation of (1), which would clash with the regularizing terms in case these latter were written in terms of just $\mu_{1},\ldots,\mu_{n}$.

4. Auxiliary results

We present here some results which will be used in the proof of Thr. 4. Define the variable $\boldsymbol{w}$ as follows:

[TABLE]

where we denoted $\mathbf{1}=(1,\ldots,1)$ and $\boldsymbol{S}=(S_{1},\ldots,S_{n})$.

Lemma 5**.**

*(Invertibility of $\boldsymbol{S}\mapsto\boldsymbol{\mu}$ and $\boldsymbol{S}\mapsto\boldsymbol{w}$)

The mappings $\Phi:\boldsymbol{S}\in{\mathcal{D}}\mapsto\boldsymbol{\mu}\in{\mathbb{R}}^{n}$, and $\Phi_{\kappa}:\boldsymbol{S}\in{\mathcal{D}}\mapsto\boldsymbol{w}\in{\mathbb{R}}^{n}$ are invertible, and their Jacobians $\Phi^{\prime}$, $\Phi_{\kappa}^{\prime}$ are uniformely positive definite in ${\mathcal{D}}$.*

Proof.

We note that $\displaystyle\frac{\partial\mu_{i}}{\partial S_{j}}=\frac{\partial^{2}{\mathcal{F}}}{\partial S_{j}\partial S_{i}}$. Direct calculation gives that

[TABLE]

from where it follows that ${\mathcal{F}}^{\prime\prime}$ is uniformly positive definite in ${\mathcal{D}}$, i.e. ${\mathcal{F}}:{\mathcal{D}}\to{\mathbb{R}}$ is a differentiable, strictly convex mapping. As a consequence, its gradient $\Phi={\mathcal{F}}^{\prime}:{\mathcal{D}}\to{\mathbb{R}}^{n}$ is a monotone (and therefore injective) mapping. Its inverse $\Phi^{-1}$ can be explicitly computed: $\Phi^{-1}(\boldsymbol{\mu})_{i}=\frac{e^{\mu_{i}}}{1+\sum_{j=1}^{n}e^{\mu_{j}}}$, $i=1,\ldots,n$, $\boldsymbol{\mu}\in{\mathbb{R}}^{n}$. Therefore $\Phi:{\mathcal{D}}\to{\mathbb{R}}^{n}$ is invertible. Moreover, since $\beta^{\prime}\geq 0$, then $\frac{\partial\boldsymbol{w}}{\partial\boldsymbol{S}}$ is symmetric and positive definite. Furthermore, $\lim_{\boldsymbol{S}\to\partial{\mathcal{D}}}|w(\boldsymbol{S})|=\infty$. Using the Hadamard global inverse theorem, [16, Thm. 2.2], we conclude that $\Phi_{\kappa}:{\mathcal{D}}\to{\mathbb{R}}^{n}$ is invertible. ∎

Lemma 6**.**

Let $f:[0,1]\to{\mathbb{R}}$ be a continuous function with $f^{\prime}(S)>0$ for $S\in(0,1)$. Given any $\boldsymbol{w}\in{\mathbb{R}}^{n}$, we denote by $\boldsymbol{S}=\boldsymbol{S}(\boldsymbol{w})\in\{\boldsymbol{S}\in(0,\infty)^{n}~{}:~{}\sum_{i=1}^{n}S_{i}<1\}$ the only solution to

[TABLE]

Then the matrix $(M_{ij}(\boldsymbol{w}))_{i,j=1,\ldots,n}=(S_{i}(\boldsymbol{w})\frac{\partial S(\boldsymbol{w})}{\partial w_{j}})_{i,j=1,\ldots,n}$ is symmetric and positive semidefinite for every $\boldsymbol{w}\in{\mathbb{R}}^{n}$.

Proof.

The definition of $\boldsymbol{S}$ implies

[TABLE]

Differentiating the above identity with respect to $w_{j}$ leads to

[TABLE]

Since $f$ is strictly increasing, then $F>0$ in $(0,1)$. It follows

[TABLE]

which means that $M(\boldsymbol{w})$ is symmetric and positive semidefinite for every $\boldsymbol{w}\in{\mathbb{R}}^{n}$. This finishes the proof. ∎

Lemma 7**.**

The following bound holds

[TABLE]

Proof.

Through simple calculations using Assumptions (H2)–(H5), (15) can be written as

[TABLE]

Since $S\in(0,1)$, the claim follows from the fact that $\beta_{1}\leq\gamma_{1}$, $\beta_{2}\leq\lambda$. ∎

The next result has been proved in [13, Lemma 5]:

Lemma 8**.**

Let $\boldsymbol{\alpha}$, $\boldsymbol{\beta}\in{\mathbb{R}}^{n}$ be such that $|\boldsymbol{\alpha}|=|\boldsymbol{\beta}|=1$. Then, for any $\boldsymbol{v}\in{\mathbb{R}}^{n}$ it holds that

[TABLE]

Notation. Let ${\mathbb{R}}_{+}\equiv[0,\infty)$. For $x\in{\mathbb{R}}_{+}\times{\mathbb{R}}^{N-1}$, we denote $x=(x_{0},\overline{x})$.

Lemma 9**.**

Let $\mathcal{R}:{\mathbb{R}}_{+}\times{\mathbb{R}}^{N-1}\to{\mathbb{R}}_{+}^{N}$ be a continuous and bounded mapping. Let $K\subset L^{2}(\Omega)$ be relatively compact. Let $\{\phi_{i}\in C_{c}^{\infty}(\Omega;{\mathbb{R}}^{N}):i\in{\mathbb{N}}\}$ be dense in $L^{2}(\Omega;{\mathbb{R}}^{N})$. Then, for every $\delta>0$, there are $C(\delta)>0$, $m(\delta)\in{\mathbb{N}}$ such that, for all $\boldsymbol{w}^{1},\boldsymbol{w}^{2}\in K\times H^{1}(\Omega;{\mathbb{R}}^{N-1})$ it holds

[TABLE]

Proof.

Assume by contradiction that there exists $\delta_{0}>0$ such that, for every $m\in{\mathbb{N}}$, there exist $\boldsymbol{w}^{1,m},\boldsymbol{w}^{2,m}\in K\times H^{1}(\Omega;{\mathbb{R}}^{N-1})$ such that

[TABLE]

Since $\mathcal{R}({\mathbb{R}}_{+}\times{\mathbb{R}}^{N-1})$ is bounded, then $(\overline{\boldsymbol{w}}^{i,m})_{m\in{\mathbb{N}}}$ is bounded in $H^{1}(\Omega;{\mathbb{R}}^{N-1})$ and thus $\overline{\boldsymbol{w}}^{i,m}\rightharpoonup\overline{\boldsymbol{w}}^{i}$ weakly in $H^{1}(\Omega;{\mathbb{R}}^{N-1})$ (as $m\to\infty$), for $i=1,2$. By a compact Sobolev embedding it holds that $\overline{\boldsymbol{w}}^{i,m}\to\overline{\boldsymbol{w}}^{i}$ strongly in $L^{2}(\Omega)$ and a.e. in $\Omega$ (up to a subsequence), for $i=1,2$. Moreover, the compactness of $K$ implies that $\boldsymbol{w}_{0}^{i,m}\to\boldsymbol{w}_{0}^{i}$ strongly in $L^{2}(\Omega)$ (up to a subsequence), for $i=1,2$. Therefore, $\boldsymbol{w}^{i,m}\to\boldsymbol{w}^{i}$ strongly in $L^{2}(\Omega;{\mathbb{R}}^{N})$ and a.e. in $\Omega$. It follows that $\mathcal{R}(\boldsymbol{w}^{i,m})\to\mathcal{R}(\boldsymbol{w}^{i})$ strongly in $L^{2}(\Omega;{\mathbb{R}}^{N})$. On the other hand,

[TABLE]

and so

[TABLE]

Being $(\phi_{i})_{i\in{\mathbb{N}}}$ dense in $L^{2}(\Omega)$, this implies that $\mathcal{R}(w^{1})=\mathcal{R}(w^{2})$. But

[TABLE]

which is a contradiction. This finishes the proof. ∎

We recall the following remark, see [5]. For completeness and clarity, we give a full proof.

Lemma 10**.**

*If a subset $\{u_{\varepsilon}\}_{\varepsilon\in(0,1]}$ of $C([0,T];L^{2}(\Omega))$ is relatively compact in

$C([0,T];L^{2}(\Omega))$, then the set $\mathcal{F}\equiv\cup_{\varepsilon\in(0,1]}\cup_{t\in[0,T]}\{u_{\varepsilon}(t)\}$ is relatively compact in $L^{2}(\Omega)$. In this case, given any $f\in C^{0}({\mathbb{R}})$, the set $\mathcal{F}_{f}\equiv\cup_{\varepsilon\in(0,1]}\cup_{t\in[0,T]}\{f(u_{\varepsilon}(t))\}$ is relatively compact in $L^{2}(\Omega)$.*

Proof.

Let $(u_{\varepsilon_{n}}(t_{n}))_{n\in{\mathbb{N}}}$ be an arbitrary sequence of points of $\mathcal{F}$. The sequence $(u_{\varepsilon_{n}})_{n\in{\mathbb{N}}}\subset C([0,T];L^{2}(\Omega))$ is relatively compact in $C([0,T];L^{2}(\Omega))$, therefore is convergent up to a subsequence. Moreover, the sequence $(t_{n})_{n\in{\mathbb{N}}}\subset[0,T]$ is convergent up to a subsequence, so w.l.o.g. we can write $u_{\varepsilon_{n}}\to u$ strongly in $C([0,T];L^{2}(\Omega))$ and $t_{n}\to t\in[0,T]$. It follows that

[TABLE]

Therefore $\mathcal{F}$ is relatively compact in $L^{2}(\Omega)$. In this case, given any $f\in C^{0}({\mathbb{R}})$, the relative compactness of $\mathcal{F}_{f}$ in $L^{2}(\Omega)$ is straightforward. This finishes the proof of the Lemma. ∎

Our main compactness tool is given in the following lemma (see Corollary 3.7. in [5]).

Lemma 11**.**

For $n\in{\mathbb{N}}$, let $w^{n}:[0,T]\to L^{2}(\Omega;{\mathbb{R}}_{+}\times{\mathbb{R}}^{N-1})$ be continuous. Assume that $K:=\{w_{0}^{n}(\cdot,t)\in L^{2}(\Omega):n\in{\mathbb{N}},t\in[0,T]\}$ is relatively compact in $L^{2}(\Omega)$, and that $\overline{w}^{n}$ is bounded in $L^{1}((0,T);H^{1}(\Omega))$. Furthermore, let $\mathcal{R}:{\mathbb{R}}_{+}\times{\mathbb{R}}^{N-1}\to{\mathbb{R}}_{+}^{N}$ be continuous and bounded. Then, $\mathcal{R}(w^{n})$ is (up to subsequence) strongly convergent in $L^{1}(\Omega\times(0,T))$.

Proof.

Apply Lemma 9. For every $\delta>0$ there exist $C(\delta)>0$, $m(\delta)\in{\mathbb{N}}$ such that, for every $n,n^{\prime}\in{\mathbb{N}}$ it holds that

[TABLE]

By integrating the above estimate in time and exploiting the boundedness of $\overline{w}^{n}$ in $L^{1}((0,T);H^{1}(\Omega))$, we deduce

[TABLE]

The boundedness of the mapping $\mathcal{R}$ implies that, up to subsequences, $\mathcal{R}(w^{n}(t))$ is weakly convergent in $L^{2}(\Omega)$ for a.e. $t\in[0,T]$, and so

[TABLE]

Moreover,

[TABLE]

The dominated convergence theorem yields

[TABLE]

It follows that $\nu\in{\mathbb{N}}$ exists such that, for $n,n^{\prime}\geq\nu$,

[TABLE]

As a consequence, it holds that

[TABLE]

In particular, $\mathcal{R}(w^{n})$ is Cauchy (and therefore convergent) in $L^{1}(\Omega\times(0,T))$. This finishes the proof. ∎

5. Existence proof

The proof is divided into several steps.

Step 1: discretization and regularization. Fix $T>0$. For $N\in{\mathbb{N}}$ we define $\kappa=T/N$, $t_{k}=\kappa k$ ($k=0,\ldots,N$), $S_{i}^{0}=S_{i,0}$ $(i=1,\ldots,n)$.

Consider the implicit Euler discretization:

[TABLE]

for all $\phi_{1},\ldots,\phi_{n}\in H^{1}_{0}(\Omega)$, where $\boldsymbol{S}^{k}:\Omega\times(0,T)\to{\mathbb{R}}^{n}$ is (implicitly) defined by

[TABLE]

and we denoted $S^{k}=\sum_{i=1}^{n}S_{i}^{k}$. Here we assume that $S^{k-1}\in H^{1}(\Omega)$.

Step 2: linearized approximated problem. Using the fact that

[TABLE]

equation (17) can be simply rewritten as

[TABLE]

Now, the linearized problem has the following form:

[TABLE]

for all $\phi_{i}\in H_{0}^{1}(\Omega)$, where $S_{i}^{*}$ is defined by

[TABLE]

and we denoted $S^{*}=\sum_{i=1}^{n}S_{i}^{*}$. The above problem can be summarized as

[TABLE]

where

[TABLE]

[TABLE]

It is easy to see that the functional $F$ is continuous, i.e. it holds

[TABLE]

The bilinear form (21) can be written as:

[TABLE]

with

[TABLE]

where

[TABLE]

Thanks to Lemma 6 and the nonnegativity of $G(S^{*})$:

[TABLE]

From Assumption (H1) we obtain

[TABLE]

Now we apply Lemma 8 and deduce

[TABLE]

Next, since \Big{(}\sum_{i=1}^{n}\sqrt{\frac{S_{i}^{*}}{S^{*}}}\Big{)}^{2}\geq n, we conclude that the bilinear form $a(\boldsymbol{w},\boldsymbol{w})$ is coercive in $H^{1}_{0}(\Omega)$, i.e.

[TABLE]

the last inequality being a consequence of Poincaré’s Lemma. Therefore we can deduce by Lax-Milgram lemma the existence of a unique solution $\boldsymbol{w}\in H^{1}_{0}(\Omega;{\mathbb{R}}^{n})$ to (19).

Remark 12**.**

We note that from the coercivity of the bilinear form $a(\boldsymbol{w},\boldsymbol{w})$ it directly follows that the solution $\boldsymbol{w}\in H^{1}_{0}(\Omega)$ to the linearized problem satisfies $\|\nabla\boldsymbol{w}\|_{L^{2}(\Omega,{\mathbb{R}}^{n})}\leq C(\varepsilon)$.

Step 3: solution of the nonlinear approximated problem. We reformulate (17) as a fixed-point problem for a suitable operator and we solve it via Leray-Schauder fixed point theorem. The Step 2 allows us to define an operator $T:L^{2}(\Omega;{\mathbb{R}}^{n})\times[0,1]\to L^{2}(\Omega;{\mathbb{R}}^{n})$ in the following way: for $\boldsymbol{w}^{*}\in L^{2}(\Omega;{\mathbb{R}}^{n})$, $\sigma\in[0,1]$, it holds that $\boldsymbol{w}=T(\boldsymbol{w}^{*},\sigma)\in H^{1}_{0}(\Omega;{\mathbb{R}}^{n})$ is the solution to (19). In a standard way we can show that the mapping $T$ is continuous. Moreover, $T:L^{2}(\Omega;{\mathbb{R}}^{n})\times[0,1]\to L^{2}(\Omega;{\mathbb{R}}^{n})$ is compact due to the compact Sobolev embedding $H^{1}(\Omega;{\mathbb{R}}^{n})\hookrightarrow L^{2}(\Omega;{\mathbb{R}}^{n})$. Furthermore, it holds that $T(\cdot,0)\equiv 0$. It remains to prove a uniform bound (with respect to $\sigma$) for all fixed points of $T(\cdot,\sigma)$ in $L^{2}(\Omega,{\mathbb{R}}^{n})$. Let $\boldsymbol{w}\in L^{2}(\Omega,{\mathbb{R}}^{n})$ be such a fixed point. Then $\boldsymbol{w}$ solves (20) with a test-function $\phi$ replaced by $\boldsymbol{w}$. We have

[TABLE]

yielding an $H^{1}$ bound for $\boldsymbol{w}$, uniform in $\sigma$. Thanks to Leray-Schauder’s fixed point theorem we get the existence of a solution $\boldsymbol{w}\in H^{1}_{0}(\Omega;{\mathbb{R}}^{n})$ to (20) for $\sigma=1$. In this way we proved the solution to (17).

Step 4: uniform in $\kappa$ a-priori estimates. Let us choose

[TABLE]

in (17). Since $\mu_{i}^{k}=\partial_{i}{\mathcal{F}}(\boldsymbol{S}^{k})$ and ${\mathcal{F}}(\boldsymbol{S})$ is convex, it follows that

[TABLE]

where $\tilde{\mathcal{F}}$ is the relative entropy density defined in (9). Moreover, the nonnegativity and boundedness of $\beta^{\prime}$ allows us to write

[TABLE]

where $\displaystyle C=\frac{1}{\max_{0\leq S\leq 1}\beta^{\prime}(S)}$. In this way we obtain

[TABLE]

Taking into account (23) and Assumption (H1), one gets

[TABLE]

Using the relation (7), we obtain

[TABLE]

In this way we get:

[TABLE]

Next, using the fact that $(a-b)a\geq\frac{1}{2}(a^{2}-b^{2}),$ we obtain

[TABLE]

We have:

[TABLE]

Next, Young inequality gives:

[TABLE]

In this way we get:

[TABLE]

Thanks to Lemma 15, we can estimate the second integral on the right-hand side of (29) by means of the third integral of the left-hand side of (29). In this way we get:

[TABLE]

Let us now introduce a new notation. Let us define the piecewise constant-in-time functions:

[TABLE]

and let $S^{(\kappa)}=\sum_{i=1}^{n}S^{(\kappa)}_{i}$. We also define the discrete backward time derivative operator $D_{\kappa}$ as follows: for every function $f:Q_{T}\to{\mathbb{R}}$,

[TABLE]

The discretized-regularized system (17) can be rewritten, in the new notation, as

[TABLE]

In the new notation, the entropy inequality (30) reads as

[TABLE]

By using the lower bounded $\displaystyle\sum_{i=1}^{n}S_{i}\log\frac{S_{i}}{S}\geq-mnS\geq-mn$, we obtain the following apriori estimates:

Proposition 13**.**

There is a constant $C$, independent of $\kappa$ and $\varepsilon$, such that

[TABLE]

for $i=1,\ldots,n$.

By using the bound (33) on the entropy function we obtain the following bounds:

Lemma 14**.**

There is a constant $C$ independent of $\kappa$ and $\varepsilon$, such that

[TABLE]

Proof.

By using simple calculations, we get

[TABLE]

The bound (40) now follows from (33). ∎

By using (36) we get the following bounds.

Lemma 15**.**

Define the exponents $\alpha_{1}$ and $\alpha_{2}$ as follows:

[TABLE]

Then, there is a constant $C$, independent of $\kappa$ and $\varepsilon$, such that:

[TABLE]

Proof.

Let us denote, for notational simplicity, $S=S^{(\kappa)}$. Then, from (36) and Assumptions (H3), (H4) we get

[TABLE]

with the constant $C$ independent of $\kappa$ and $\varepsilon$. As a consequence

[TABLE]

The inequality stated above implies the following bound for the functions $Z\equiv\min(S,1/2)$, $W\equiv\max(S,1/2)$:

[TABLE]

which can be written as

[TABLE]

The function $Z^{\alpha_{1}}$ and $(1-W)^{\alpha_{2}}$ are in $L^{2}(0,T;L^{2}(\Omega))$ due to Assumption (H5) and Lemma 14. Indeed, Assumption (H5) implies that $2\alpha_{1}\geq 2-\gamma_{1}$ and $2\alpha_{2}\geq 2-\lambda$. We can then use the Sobolev embedding theorem to get the bound:

[TABLE]

Due to (41) these bounds hold also for the function $S$ instead of $Z$ and $W$. ∎

Lemma 16**.**

There exists $p>1$ such that

[TABLE]

where the constant $C>0$ is independent of $\kappa$ and $\varepsilon$.

Proof.

For simplifying the notation we will write $S=S^{(\kappa)}$. We first notice that

[TABLE]

So it is sufficient to prove that $S^{-p\gamma}$ and $(1-S)^{-p\lambda}$ are uniformely bounded in $L^{1}(\Omega\times(0,T))$ for some $p>1$.

It is clear that integrability given by Lemma 14 is not sufficient to prove the estimate (43). Therefore, we will combine estimates from Lemmas 15 and 14 in order to obtain the integrability with requested exponents. Assumptions (H5) on the parameters $\beta_{1}$, $\beta_{2}$, $\gamma$, $\gamma_{1}$ and $\lambda$ imply

[TABLE]

We rewrite the expression $\int_{\Omega}S^{-\gamma p}dx$ using $-\gamma p=\alpha_{1}\Theta+(2-\gamma_{1})\Theta_{1}$ and Hölder’s inequality:

[TABLE]

We take $p_{1}=6/\Theta$ and $p_{2}=6/(6-\Theta)$, $\Theta=2$, $\Theta_{1}=2/3$ and we get

[TABLE]

Because of (42) and (40), the right hand side is uniformly bounded. Condition

[TABLE]

is equivalent to

[TABLE]

Now it is easy to see that (45) and the first inequality in (44) are equivalent to the first inequality in Assumption (H5).

The second inequality in Assumption (H5) in treated in the same way. The calculations are given here for completeness. We rewrite the expression $\int_{\Omega}(1-S)^{-\lambda p}dx$ using $-\lambda p=\alpha_{2}\Theta+(2-\lambda)\Theta_{1}$ and Hölder’s inequality:

[TABLE]

We take $p_{1}=6/\Theta$ and $p_{2}=6/(6-\Theta)$, $\Theta=2$, $\Theta_{1}=2/3$ and obtain

[TABLE]

Because of (42) and (40), the right hand side is uniformly bounded. Condition

[TABLE]

is equivalent to $\lambda<3\beta_{2}-10$. It is now easy to see that this inequality together with the second inequality in (44) are equivalent to the second inequality in Assumption (H5). This concludes the proof of Lemma 16. ∎

Proposition 17**.**

There is an exponent $1<q<2$ such that

[TABLE]

where $C$ is a constant independent of $\kappa$ and $\varepsilon$.

Proof.

Let $p>1$ as in Lemma 16. By choosing $q=2p/(1+p)\in(1,2)$, we get

[TABLE]

By using Lemma 16 and bound (37) we conclude the proof. ∎

Finally, from equation (31) together with the the bounds (36)–(39), we get the following uniform bound for the discrete time derivative:

[TABLE]

5.1. Passing to the limit when $\kappa\to 0$

From (38) and Lemma 5 we get that

[TABLE]

From this and the bound from the discrete time derivative (47) we get by using the nonlinear version of the Aubin-Lions lemma [4] that

[TABLE]

This strong convergence holds also in $L^{q}(0,T;L^{q}(\Omega))$ for any $q<\infty$.

By using the bounds in Proposition 13 and Proposition 17, we obtain that the solution of (31) satisfies

[TABLE]

Thus, after taking the limit $\kappa\to 0$, (43) holds with $S^{(\varepsilon)}$ in place of $S^{(\kappa)}$, i.e.

[TABLE]

Also, estimates (33)-(39), (46), (47) hold with $S^{(\varepsilon)}$ in place of $S^{(\kappa)}$, i.e.

[TABLE]

for $i=1,\ldots,n$.

5.2. Passing to the limit $\varepsilon\to 0$

Now we define the continuous mapping $\mathcal{R}:{\mathbb{R}}_{+}\times{\mathbb{R}}^{n}\to{\mathbb{R}}^{n+1}$ as

[TABLE]

It follows from (12) that $S_{i}^{(\varepsilon)}=\mathcal{R}_{i}(S^{(\varepsilon)},(\mu^{*})^{(\varepsilon)})$ for $i=1,\ldots,n$. Lemma 11 implies that $S_{i}^{(\varepsilon)}$ has a subsequence that is strongly convergent in $L^{1}(\Omega\times(0,T))$, for $i=1,\ldots,n$. From the $L^{\infty}$ bounds for $S_{i}^{(\varepsilon)}$ we conclude that, up to a subsequence, it holds that

[TABLE]

By using this convergence property as well as the bounds (49)–(58), we are able to take the limit $\varepsilon\to 0$ in (48) and obtain that $\boldsymbol{S}=(S_{1},\ldots,S_{n})$ is a weak solution to (1)–(4). This finishes the proof of the theorem.

Appendix

Derivation of the model

We consider an isothermal, immiscible and incompressible two-phase flow of water and oil in a porous media, where oil consists of $n$ chemical components. Let us denote by $\mathcal{V}$ the representative volume (REV), which consists of the solid part $\mathcal{V}_{s}$ and the pore space $\mathcal{V}_{p}$. The flow occurs in a porous domain $\mathcal{V}_{p}$ of volume $\Delta V_{p}$, where the porosity (the relative volume occupied by the pores) is denoted by $\displaystyle\Phi=\Delta V_{p}/\Delta V$. The saturations of the oil and water phase are given by $S_{\alpha}=\Delta V_{\alpha}/\Delta V_{p}$, where $\Delta V_{\alpha}$ is the volume of the phase $\alpha$ with $\alpha=w,o$. Following [2], a generalized Darcy law gives

[TABLE]

Here the subscripts $w$ and $o$ correspond, respectively, to the water (wetting) and the oil (non-wetting) fluids, $\mathbf{u}_{\alpha}$ are the fluxes of the phases, $p_{\alpha}$ are their pressures, and $\lambda_{\alpha}$ are the phase mobilities. We assume that $\lambda_{\alpha}$ depend on the nonwetting-phase saturation $S_{o}$. Furthermore, $k$ is the absolute permeability of the porous medium, and the gravity effects are neglected for simplicity. The mass conservation laws for both phases have the form:

[TABLE]

where $\Phi$ is the porosity of the medium. The model (59)–(60) has to be completed with the capillary pressure law which has the form

[TABLE]

where, due to [10], the capillary pressure saturation relationship is given by

[TABLE]

Here, $p_{c}(S_{o})$ is the static capillary pressure function and $\tau(S_{o})$ is the relaxation parameter.

We assume that the non-wetting phase (oil) is a heterogeneous mixture of hydrocarbon compounds and we derive the mass conservation equation for each compound. More precisely, in the oil phase there are $n$ components whose mass concentrations $c_{o}^{i}$, i.e. the densities of the $i$-th component in the volume of the phase, are given by $c_{o}^{i}=\Delta m_{o}^{i}/\Delta V_{o}$, where $\Delta m_{o}^{i}$ is the mass of the component $i$ in the oil-phase of the REV. The sum of the mass concentrations of all components is given by

[TABLE]

Noting that

[TABLE]

the mass conservation equation for the component $i$ is given by

[TABLE]

The component velocities $\mathbf{u}_{o,i}$ are related to the phase velocity $\mathbf{u}_{o}$ by the expression

[TABLE]

The flux of the oil-phase components consists of the relative movement of the constituents $i$ spreading due to random collisions between molecules of different types (diffusion) followed by the convection, i.e.

[TABLE]

Note that $\sum_{i=1}^{n}\mathbf{J}_{o}^{i}=\mathbf{0}$. Let us introduce the saturation of the component $i$ in the oil phase as

[TABLE]

It is clear that $\sum_{i=1}^{n}S_{o}^{i}=S_{o}$. Furthermore, we assume that each component $i$ of the mixture in the oil phase is incompressible, i.e.

[TABLE]

Now we have

[TABLE]

Next, we make the assumption that the diffusion fluxes are proportional to the spatial gradients of suitable chemical potentials, i.e.

[TABLE]

where $\mu_{j}$ are given by (2) using the notation $S\equiv S_{o}$. In this way, equation (63) reads:

[TABLE]

Now, a simple calculation gives

[TABLE]

Furthermore, we assume that the total flow equals zero, i.e. $\mathbf{u}_{0}+\mathbf{u}_{w}=\bf{0}$, which gives

[TABLE]

Here the diffusion mobility $a(S_{o})$ is given by

[TABLE]

In this way, we obtain the parabolic system of our interest

[TABLE]

where $i=1,\ldots,n$. Notice that (1) is identical to (67) with $k=1$ and $S_{o}$, $S_{o}^{i}$ replaced by $S$, $S_{i}$, respectively.