The unsaturated flow in porous media with dynamic capillary pressure

TL;DR

This paper proves the existence of weak solutions for a degenerate pseudoparabolic PDE modeling unsaturated flow in porous media with dynamic capillary pressure, using Galerkin approximation and entropy methods.

Contribution

It introduces a rigorous mathematical framework for a complex PDE with saturation-dependent relaxation, extending previous models with new existence results.

Findings

Existence of weak solutions established

A priori estimates derived using entropy-inspired test functions

Handling of mixed derivatives in the PDE achieved

Abstract

In this paper we consider a degenerate pseudoparabolic equation for the wetting saturation of an unsaturated two-phase flow in porous media with dynamic capillary pressure-saturation relationship where the relaxation parameter depends on the saturation. Following the approach given in [12] the existence of a weak solution is proved using Galerkin approximation and regularization techniques. A priori estimates needed for passing to the limit when the regularization parameter goes to zero are obtained by using appropriate test-functions, motivated by the fact that considered PDE allows a natural generalization of the classical Kullback entropy. Finally, a special care was given in obtaining an estimate of the mixed derivative term by combining the information from the capillary pressure with obtained a priori estimates on the saturation.

Equations

$$\displaystyle\partial_{t}S=\operatorname{div}\Big{(}a(S)\nabla\big{(}\partial_{t}\beta(S)-P_{c}(S)\big{)}\Big{)},$$

$$a(0)$$

$$-P_c^{\prime }(S)\to +\infty$$

$$\beta (S)$$

$$\tau (0)=0,\tau (S)\text{ is bounded when }$$

$$\displaystyle\partial_{t}S=\operatorname{div}\Big{(}k(S)\nabla\big{(}\tau\partial_{t}S-P_{c}(S)\big{)}\Big{)}.$$

$$\displaystyle\frac{d}{dt}\int_{\Omega}\Big{(}\mathcal{E}_{M}(S)+\frac{\tau}{2}|\nabla S|^{2}\Big{)}dx=\int_{\Omega}P_{c}^{\prime}(S)|\nabla S|^{2}dx+\textrm{ low order terms. }$$

$$\displaystyle\partial_{t}S=\operatorname{div}\Big{(}a(S)\nabla\big{(}\tau(S)\partial_{t}S-P_{c}(S)\big{)}\Big{)}.$$

$$\mathcal{E}(S)=S{\int }_{S_D}^S\frac{\tau (\xi )}{a(\xi )}d\xi -{\int }_{S_D}^S\xi \frac{\tau (\xi )}{a(\xi )}d\xi ,$$

$${\mathcal{E}}^{\prime }(S)={\int }_{S_D}^S\frac{\tau (\xi )}{a(\xi )}d\xi ,$$

$${\mathbf{u}}_w=-{\lambda }_w(S)k\nabla p_w,{\mathbf{u}}_o=-{\lambda }_o(S)k\nabla p_o.$$

$$\varphi \frac{\partial S}{\partial t}+\mathrm{div}{\mathbf{u}}_w=0,\varphi \frac{\partial (1-S)}{\partial t}+\mathrm{div}{\mathbf{u}}_o=0,$$

$$p_o-p_w=p_c,$$

$$p_c=P_c(S)-\tau (S)\frac{\partial S}{\partial t}.$$

$$\displaystyle\frac{\partial S}{\partial t}+\operatorname{div}\big{(}b(S)\mathbf{u}\big{)}=\operatorname{div}\big{(}a(S)\nabla(-p_{c})\big{)},$$

$$b(S)=\frac{{\lambda }_w(S)}{{\lambda }_w(S)+{\lambda }_o(S)},a(S)=\frac{{\lambda }_w(S){\lambda }_o(S)}{{\lambda }_w(S)+{\lambda }_o(S)}.$$

$$\mathbf{u}={\mathbf{u}}_o+{\mathbf{u}}_w=0.$$

$$\displaystyle\partial_{t}S+\operatorname{div}\Big{(}a(S)P_{c}^{\prime}(S)\nabla S\Big{)}=\operatorname{div}\Big{(}a(S)\nabla\big{(}\tau(S)\partial_{t}S\big{)}\Big{)}.$$

$$S$$

$$\displaystyle a(S)\nabla\big{(}P_{c}(S)-\partial_{t}\beta(S)\big{)}\cdot\nu$$

$$S(x,0)$$

$$a(s)=\frac{s^{\mu }\cdot (1-s{)}^{\lambda }}{s^{\mu }+(1-s{)}^{\lambda }},$$

$$\displaystyle\tau(s)=\frac{s^{\mu}}{s^{\mu}+(1-s)^{\lambda}}\Big{[}{\mathcal{T}}_{M}+\frac{(1-s)^{\lambda}}{s^{\gamma}}\Big{]}.$$

$$0<{\beta }_1\le \gamma <\mu \text{ and }0<{\beta }_2\le \lambda$$

$$-P_c^{\prime }(s)=\frac{g(s)}{s^{{\beta }_1}}+\frac{h(s)}{(1-s{)}^{{\beta }_2}},$$

$$-P_c^{\prime }(s)\to +\infty \text{ when }s\to 0^{+},s\to 1^{-}.$$

$$\mathcal{R}(x,t,s)=R_0(x,t)\sigma (s),\text{ where }{R}_0\in C^1({\overline{\Gamma }}_N),\sigma \in C_c^{\infty }(0,1).$$

$$\displaystyle\gamma>2,\;\lambda>2\;\textrm{ and }\;\int_{\Omega}S_{i}^{2-\gamma}(x)dx+\int_{\Omega}\big{(}1-S_{i}(x)\big{)}^{2-\lambda}dx<+\infty.$$

$$a(s)=\frac{s^{\mu}\cdot(1-s)^{\lambda}}{f(s)\big{(}s^{\mu}+(1-s)^{\lambda}\big{)}},\;\textrm{ and }\;\tau(s)=\frac{s^{\mu}}{f(s)\big{(}s^{\mu}+(1-s)^{\lambda}\big{)}}\Big{[}{\mathcal{T}}_{M}+\frac{(1-s)^{\lambda}}{s^{\gamma}}\Big{]},$$

$$\beta (S)={\int }_0^S\tau (s)ds,$$

$$\tau (s)=\{\begin{array}{ll}0,& s<0,\\ \tau (1),& s>1.\end{array}$$

$$\beta (s)=\{\begin{array}{ll}0,& s<0,\\ \beta (1)+\tau (1)(s-1),& s>1.\end{array}$$

Full text

1. Introduction

Two-phase flow processes in porous media appear in many real situations, such as unsaturated water flow in the subsurface, remediation of contaminated sites and water-oil displacement during oil recovery. The physical-mathematical model of the two-phase flow consists of two mass balance equations, Darcy relations for the two phases and the capillary condition which is a constitutive relationship between the wetting-phase saturation and capillary pressure $P_{c}$. In this paper we include the dynamic capillary pressure relation proposed by Hassanizadeh and Grey [11] to the standard two-phase flow equations. This results in a transport equation containing higher-order mixed derivatives. Moreover, the obtained model represents degenerate pseudoparabolic equation and our goal is to prove the global-in-time existence of its weak solutions. More specifically, we consider the equation

[TABLE]

where

[TABLE]

[TABLE]

The solution of equation (1) satisfies certain initial-boundary conditions which will be specified later.

We note that the nonequilibrium capillary effects in problems of enhancing oil and gas recovery from rocks were proposed in classical book by Barenblatt, Entov and Ryzhik [1], and then investigated by numerous scientists (see for example [2, 3, 4, 16], [6, 10]). The mathematical analysis of dynamic capillary pressure models for Richards’ equation was intensively studied. A significant contribution has been done by Mikelić [13] where the existence, for any time interval, of an appropriate weak solution of the Richards’ equation was proved. The existence of a weak solution to the Richards’ equation with dynamic capillary pressure and hysteresis is studied by Schweizer [14]. For the two-phase flow model with dynamic capillary pressure and hysteresis, the existence of a weak solution is proved by Koch, Rätz and Schweizer [12] under assumption of non-degenerated mobilities.

The first existence result for the two-phase flow model with dynamic capillary pressure and saturation dependent relaxation parameter is obtained by Cao and Pop in [7]. They considered the system in global pressure formulation with a strictly positive relaxation parameter and a constant Dirichlet boundary condition for the saturation.

In this work we consider the two-phase flow model with dynamic capillary pressure and degenerate saturation-dependent relaxation parameter. Moreover, we impose mixed boundary conditions with a non-constant boundary saturation. As our result shows (see also [13] and [7]), the existence theorem can be proved under certain relations between the orders of zeros of relative permeabilities and relaxation parameter and the order of singularities of capillary pressure function (see Proposition 3). Concerning these relations it is sufficient to analyze the case of the countercurrent imbibition flow instead of the full two-phase flow system. Therefore, we consider the countercurrent imbibition flow in the two-phase system that changes wettability (see [9]).

Key ideas. The motivation for this article comes from the result of Mikelić [13] on a global existence for the following nonlinear degenerate pseudoparabolic equation

[TABLE]

Mikelić observed that equation (2) allows a natural generalization of the classical Kullback entropy given by $\mathcal{E}_{M}^{\prime\prime}(s)=1/k(s)$. Namely, formal calculation gives

[TABLE]

The main difficulties for the mathematical analysis of (2) were the degenerations of the coefficients, i.e. $k(0)=0$ and $-P_{c}^{\prime}(S)\to+\infty$ when $S\to 0$. These degenerate coefficients and presence of the initial and the boundary conditions lead to unbounded non-integrable $\mathcal{E}_{M}^{\prime}$. Because of that reason, in [13] the degenerate coefficients were firstly regularized, by introducing the regularization parameter $\varepsilon$. After that, the entropy estimates necessary for passing to the limit in the weak formulation when $\varepsilon\to 0$ were derived. The crucial step in mentioned existence result was the use of $\mathcal{E}_{M}^{\prime}$ as a test-function for obtaining needed a priori estimates.

In comparison with equation (2) considered in [13], besides the additional degeneration of the coefficient $a(S)$ in $S=1$, and the fact that $-P_{c}^{\prime}(S)\to+\infty$ for $S\to 0$ and $S\to 1$, the main novelty in this article is the dependence of the relaxation parameter $\tau$ upon the variable $S$, i.e. after rewriting, the equation of our interest (1) reads

[TABLE]

Similarly like in [13], we observe that equation (4) allows a natural generalization of the classical Kullback entropy. For the equation of our interest (1), the entropy functional is given by formula

[TABLE]

where $S_{D}$ comes from the boundary condition. Easy calculation gives $\mathcal{E}^{\prime\prime}(S)=\tau(S)/a(S)$. In order to obtain needed a priori estimates, adapting the idea given in [13], we use

[TABLE]

as a test-function in a weak formulation given later by (19), leading to an estimate similar to (3).

2. Physical model and basic equations

We consider two-phase water-oil flow (wetting and non-wetting immiscible fluids) in an isotropic and homogeneous porous medium. Following [5], the 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}_{i}$ are the fluxes of the phases, $p_{i}$ are their pressures, and $\lambda_{i}$ are phase mobilities, depending on the wetting-phase saturation $S=S_{w}$. 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 (6)–(7) has to be completed with the capillary pressure law which has the form

[TABLE]

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

[TABLE]

Here, $P_{c}(S)$ is static capillary pressure function and $\tau(S)$ is the relaxation parameter. Expression (8) suggests that a new equilibrium between the dynamic and the static capillary pressure is not attained instantaneously and the relaxation parameter $\tau$ depends on the water saturation. Following the work of Cuesta, van Duijn and Hulshof [8] and the references therein [10, 15], we take that $\tau(S)$ vanishes as $S\to 0$. This assumption can be explained as in [4]. Namely, at low water saturation the water flows to the narrowest pores and the time needed to reconfigure the water distribution tends to zero. We note that this degeneration introduces additional mathematical difficulty into the system in comparison with [13]. For purposes of our mathematical analysis (see inequality (52)), we take that $\tau(S)$ is monotone increasing function.

Hereafter, we normalize the porosity $\phi$ and the permeability $k$ to one. By adding equations (7) we get

[TABLE]

where the coefficients are given by

[TABLE]

For the simplicity reasons we consider here the capillary countercurrent imbibition flow (see [4]) where, due to the incompressibility of both fluids, the bulk fluid flow is identically equal zero:

[TABLE]

More precisely, we study the following equation for the water saturation:

[TABLE]

which is actually the equation of our interest, written in (1).

3. Equation and its weak formulation

Let $\Omega\subset{\mathbb{R}}^{n}$ be an open, bounded domain, with Lipschitz boundary $\partial\Omega$. The boundary $\partial\Omega$ is decomposed into the Dirichlet part $\partial_{D}\Omega$ and the Neumann part $\partial_{N}\Omega$, where $\partial\Omega=\overline{\partial_{D}\Omega}\cup\overline{\partial_{N}\Omega}$ and $\partial_{D}\Omega\cap\partial_{N}\Omega=\emptyset$. We suppose that the Dirichlet boundary $\partial_{D}\Omega$ is closed relative to $\partial\Omega$ and has a positive $(n-1)$-dimensional Hausdorff measure. The Neumann boundary $\partial_{N}\Omega$ may be empty. Let $Q_{T}=\Omega\times(0,T)$, $T>0$. Moreover, we use the notation $\Gamma_{D}=\partial_{D}\Omega\times(0,T)$ and $\Gamma_{N}=\partial_{N}\Omega\times(0,T)$. Here $(0,T)$, $T>0$, is the time interval. We are looking for a solution to equation (1) which satisfies the following boundary and initial conditions:

[TABLE]

where $\nu$ is the unit outer normal on $\partial\Omega$, $\mathcal{R}$ is a given flux, $S_{D}$ is the given saturation at $\Gamma_{D}$ and $S_{i}$ is the initial saturation of the wetting phase. We impose following hypotheses on the coefficients and on the data, for $s\in[0,1]$:

(H1) There are constants $\mu>0$ and $\lambda>0$ such that

[TABLE]

giving $a(0)=a(1)=0$.

(H2) For the relaxation parameter we suppose that there exist constant $\gamma>0$ and ${\mathcal{T}}_{M}>1$, such that

[TABLE]

Furthermore, we take $\mu>\gamma\geq 0$ so we get $\tau(0)=0$ and $\tau(1)={\mathcal{T}}_{M}$. The condition ${\mathcal{T}}_{M}>1$ together with $\mu>\gamma$ ensures that the relaxation parameter $\tau$ is increasing function of water saturation $s$.

(H3) For the derivative of the capillary pressure we suppose that there exist constants

[TABLE]

such that

[TABLE]

where $g,h\in C^{\infty}_{c}({\mathbb{R}})$ are arbitrary functions strictly positive on $[0,1]$. Note that

[TABLE]

(H4) We assume that $S_{D}\in C^{1}([0,T];H^{1}(\Omega))$, $0<S_{D_{min}}\leq S_{D}(x,t)\leq S_{D_{max}}<1$ a.e. on $Q_{T}$. Moreover, we assume that $\displaystyle\|\partial_{t}S_{D}\|_{L^{\infty}(Q_{T})}<+\infty$. The initial condition $S_{i}$ belongs to $H^{1}(\Omega)$ and $0\leq S_{i}\leq 1$ a.e. on $\Omega$.

(H5) We assume that the boundary-flux $\mathcal{R}$ satisfies

[TABLE]

(H6) The initial saturation satisfies the finite entropy condition:

[TABLE]

Remark 1**.**

In order to avoid technical difficulties we use explicit forms for $a$ and $\tau$. This can be easily generalized, for example using the strictly positive on $[0,1]$ function $f\in C^{\infty}_{c}({\mathbb{R}})$ such that

[TABLE]

where $f$ should be chosen in a way that the relaxation parameter $\tau(s)$ is monotone increasing function.

The function $\beta(S)$ is defined as

[TABLE]

and it is bounded: $0\leq\beta(S)\leq\beta(1)$. We note that we extend $\tau$ to $(-\infty,0)$ and to $(1,+\infty)$ as follows:

[TABLE]

So we have

[TABLE]

From hypotheses (H1)–(H3) it is easy to see that the following coefficients are bounded:

[TABLE]

[TABLE]

Our goal is to obtain a global existence of a weak solution for problem (1)–(10)–(11)–(12), for any time interval. Firstly, we define the weak formulation of the problem.

Definition 1**.**

Let

[TABLE]

Find $S\in L^{2}(Q_{T})$, $0\leq S(x,t)\leq 1$ a.e. on $Q_{T}$ such that $\beta(S)\in H^{1}(Q_{T})$, $\beta(S)-\beta(S_{D})\in L^{2}(0,T;V)$, $a(S)\nabla\partial_{t}\beta(S)\in L^{2}(Q_{T})$ and that for all $\varphi$ in $H^{1}(0,T;V)$, $\varphi_{|_{t=T}}=0$ it holds

[TABLE]

In the following theorem we state the main result of the paper.

Theorem 1**.**

Global existence of a weak solution.

Let $n\leq 3$ and let us suppose hypotheses (H1)–(H6). Under the assumptions on the exponents $\beta_{1}$, $\gamma$, $\mu$, and $\lambda$ given by Proposition 3, there exists a weak solution $S$ for problem (1)-(10)-(11)-(12) such that $0\leq S\leq 1$.

The rest of the paper is organized as follows. In Section 4 we firstly regularize the coefficients and define a weak formulation of the regularized problem. Then by Theorem 2 we state the global existence of a weak solution for the regularized problem. Section 5 gives proof of Theorem 2. More precisely, we start with the Galerkin method for the regularized problem, then in Theorem 3, Section 6 we give a priori estimates uniform with respect to the regularization parameter $\varepsilon$. Special care was given in obtaining a priori estimate of the mixed-derivative term in Section 7, where we combined informations from the capillary pressure with those obtained from the entropy functional. Finally, in Section 8 we pass to the limit when $\varepsilon\to 0$ finishing the proof of Theorem 1.

4. Regularization of the coefficients

For proving existence of a weak solution given by Definition 1 we need to define a regularized problem. For that purpose we regularize the coefficients $a$, $P_{C}^{\prime}$ and $\tau$ and extend them to the whole ${\mathbb{R}}$ as follows. Let us introduce the cutting-function $Z$ defined for $s\in\mathbb{R}$ as

[TABLE]

and, for $\varepsilon>0$, the function

[TABLE]

Then we define

[TABLE]

where $s\in{\mathbb{R}}$. We also define

[TABLE]

It is easy to see that the regularized coefficients are bounded, i.e.

[TABLE]

We note that constants $m_{a}^{\varepsilon}$, $m_{aP}^{\varepsilon}$, $m_{\tau}^{\varepsilon}$ depend on $\varepsilon$. Moreover, $m_{\tau}^{\varepsilon}=\tau(\varepsilon)$. Furthermore, we have

[TABLE]

Next, we introduce the notion of a weak solution for the regularized version of problem (1)-(10)-(11)-(12).

Definition 2**.**

Find $S_{\varepsilon}\in H^{1}(Q_{T})$ such that $S_{\varepsilon}-S_{D}\in L^{2}(0,T;V)$, $\nabla\partial_{t}\beta(S_{\varepsilon})\in L^{2}(Q_{T})$ and satisfying

[TABLE]

for all $\varphi\in H^{1}(0,T;V)$ such that $\varphi_{|_{t=T}}=0$, and $S_{\varepsilon}(0)=S_{i}$ in $L^{2}(\Omega)$.

Theorem 2**.**

Global existence of a weak solution for the regularized problem.

Under hypotheses (H1)–(H5)* there exists a solution for the regularized problem (29).*

5. Proof of Theorem 2

Step 1: Galerkin approximation.

Let $(e_{j})_{j\in{\mathbb{N}}}$ be a basis of $V$ and let $V_{N}=\operatorname{span}\{e_{1},\ldots,e_{N}\}$. The Galerkin approximation for the regularized problem (29) reads as follows:

find $\displaystyle\beta_{\varepsilon}(S_{\varepsilon,N})=\beta_{\varepsilon}(S_{D})+\sum_{j=1}^{N}\gamma_{j}(t)e_{j}(x)$, $\gamma_{j}\in C^{1}([0,T])$ such that

[TABLE]

for all $l\in\{1,\ldots,N\}$ satisfying the initial condition

[TABLE]

where $\Pi_{N}$ is the orthogonal projector in $L^{2}(\Omega)$ on the space $V_{N}$.

Lemma 1**.**

There exists $T_{N}>0$ such that problem (30) has a unique solution which satisfies $\beta_{\varepsilon}(S_{\varepsilon,N})-\beta_{\varepsilon}(S_{D})\in C^{1}([0,T_{N}];V_{N})$.

Proof. We define the mapping $\Phi_{\varepsilon}=\beta_{\varepsilon}^{-1}$ such that S_{\varepsilon,N}=\Phi_{\varepsilon}\Big{(}\beta_{\varepsilon}(S_{D})+\sum_{i=1}^{N}\gamma_{i}(t)e_{i}\Big{)}. For simplicity of notation, we omit indices $\varepsilon$ and $N$ from $S_{\varepsilon,N}$. We can rewrite (30) as:

[TABLE]

where $A\in{\mathbb{R}}^{N\times N}$, $B\in{\mathbb{R}}^{N\times N}$ and $F\in{\mathbb{R}}^{N}$. For the coefficients in (32) we have

[TABLE]

[TABLE]

The matrix $A=A(\gamma)$ is symmetric, positive definite, depending smoothly on $\gamma$. Namely, one has

[TABLE]

where $u_{\xi}=\sum_{i}\xi_{i}e_{i}$ and regularized coefficient satisfies $a_{\varepsilon}(S)\geq m_{a}^{\varepsilon}>0$ and $\beta_{\varepsilon}^{\prime}(S)=\tau_{\varepsilon}(S)\geq m_{\tau}^{\varepsilon}>0$. Furthermore, since $F=F(\gamma,t)$ and $B=B(\gamma)$ are continuously differentiable functions of $\gamma$ and continuous functions of $t$, the Cauchy-Lipschitz theorem gives the local existence of the solution. $\hfill\Box$

Step 2: a priori estimates uniform with respect to $N$ for the regularized finite-dimensional problem.

Proposition 1**.**

Under assumptions (H1)–(H5), the solution of the approximate problem (30)–(31) exists for all times $T<+\infty$ and \beta_{\varepsilon}(S_{\varepsilon,N})-\beta_{\varepsilon}(S_{D})\in C^{1}\big{(}[0,T];V_{N}\big{)}.

Proof. Let us show boundedness of the solution, uniform with respect to $N$. For the test-function in (30) we take \partial_{t}\big{(}\beta_{\varepsilon}(S_{\varepsilon,N})-\beta_{\varepsilon}(S_{D})\big{)}. In this way we obtain

[TABLE]

Using the assumptions on the boundedness of the regularized coefficients (25)–(26), the Cauchy inequality, the fact that all integrals which contain $S_{D}$ are bounded by some constant, and the trace theorem, we get the following estimate

[TABLE]

where $C_{\varepsilon}$ is independent on $N$. Next, from the expression

[TABLE]

by integrating over the set $\Omega$, one gets:

[TABLE]

where $C_{\varepsilon}$ is independent of $N$. From (34) and (35) we get

[TABLE]

Now, the Gronwall inequality gives

[TABLE]

where $C_{\varepsilon}$ is independent of $N$. Finally, from (36) we get

[TABLE]

We conclude that, since $\partial_{t}\beta_{\varepsilon}(S_{\varepsilon,N})$ is bounded in $L^{\infty}(0,T;H^{1}(\Omega))$, uniformly with respect to $N$, the solution for the approximate problem (30)-(31) exists for all times $T<+\infty$ and one has $\beta_{\varepsilon}(S_{\varepsilon,N})-\beta_{\varepsilon}(S_{D})\in C^{1}([0,T];V_{N})$. $\Box$

From now on we consider problem (30)–(31) on $[0,T]$, for fixed $T>0$. The calculations done in the proof of Proposition 1 give the following a priori estimates uniform with respect to $N$.

Corollary 1**.**

Under hypotheses (H1)–(H5)* there is a constant $C_{\varepsilon}$ independent of $N$ such that*

[TABLE]

Proof. Estimate (41) follows directly from (37), while estimate (40) follows from (35) and (41). Moreover, estimate (37) gives that

[TABLE]

and, since $\tau_{\varepsilon}$ is bounded (see (26)), estimate (39) follows immediately. Note that (42) follows from (41) due to the boundedness of the coefficient $a_{\varepsilon}$. The estimate (38) follows from (40) and boundedness of regularized coefficients given by (26). $\hfill\Box$

Step 3: passing to the limit when $N\to\infty$.

Using the a priori estimates obtained in Corollary 1 we are able to pass to the limit in the weak formulation when $N\to\infty$. More precisely, under hypotheses (H1)–(H5), the solution $S_{\varepsilon,N}$ of the approximate problem (30)-(31) converges to a function $S_{\varepsilon}\in H^{1}(Q_{T})$ when $N\to+\infty$ satisfying equation (29). Moreover, we have $\partial_{t}\beta_{\varepsilon}(S_{\varepsilon})\in L^{2}(Q_{T})$. Now, it is straightforward to prove that the limit function $S_{\varepsilon}$ satisfies the weak formulation (29) as well as the initial and boundary conditions. This completes the proof of Theorem 2. $\hfill\Box$

6. A priori estimates uniform with respect to $\varepsilon$

For clarity and better readability the proof of Theorem 1 is written in few steps. In this section we state a priori estimates uniform with respect to the regularizing parameter $\varepsilon$ needed for passing to the limit when $\varepsilon\to 0$. The needed a priori estimates are obtained using two different test-functions in weak formulation (29). Furthermore, we prove that the integral of the entropy functional given by expression (5) is bounded.

The key-result for proving the global existence of a weak solution for the regularized problem are the a priori estimates given by the following theorem.

Theorem 3**.**

A priori estimates uniform with respect to $\varepsilon$.*

Under the hypotheses (H1)–(H6) there is a constant $C>0$, independent of $\varepsilon$, such that every solution $S_{\varepsilon}$ for the regularized problem (29) satisfies*

[TABLE]

where $\mathcal{E}_{\varepsilon}$ is given below by expression (49).

Proof. We test variational equation (29) by using two different test-functions.

Test-function 1. The first choice for the test-function is given by

[TABLE]

We note that this selection for $\varphi_{\varepsilon}$ is connected to the regularized entropy, corresponding to \big{(}k_{\varepsilon}(S_{\varepsilon})\big{)}^{-1} in [13]. From Section 1, Subsection Key ideas, we recall that $\mathcal{E}_{\varepsilon}^{\prime\prime}(S)=\tau_{\varepsilon}(S)/a_{\varepsilon}(S)$ where

[TABLE]

so we have $\varphi_{\varepsilon}(S)=\mathcal{E}_{\varepsilon}^{\prime}(S)$.

Using the test-function given by (48), from (29) we obtain the following expression

[TABLE]

We denote

[TABLE]

The Cauchy inequality gives

[TABLE]

Using the mean value theorem and the assumption that $\tau_{\varepsilon}(s)$ is monotone increasing function, it is easy to derive the following inequality

[TABLE]

After inserting (52) into (51), using the Poincaré inequality and the assumptions that $S_{D}\geq S_{D_{min}}>0$ one gets:

[TABLE]

where the constant $C_{p}$ comes from the Poincaré inequality. By rewriting the term

[TABLE]

after using the Cauchy-Schwarz and then the Cauchy inequality we get

[TABLE]

The third term on the left-hand side in (50)

[TABLE]

we rewrite as

[TABLE]

Next, we estimate the integral

[TABLE]

After applying the Cauchy inequality we have

[TABLE]

It remains to estimate the integral

[TABLE]

Similarly like before, the Cauchy inequality gives

[TABLE]

Finally, concerning the estimate of the integral in (50) which contains the boundary-flux term $\mathcal{R}$, due to assumption ${\bf(H5)}$, we firstly note that

[TABLE]

where we denoted by $\operatorname{supp}$ the support of the function and we take the set $A\subset(0,1)$ such that $\operatorname{supp}S_{D}\cup\operatorname{supp}(\sigma)\subseteq A$. In this way we conclude that the boundary-flux integral

[TABLE]

is bounded by a positive constant independent of $\varepsilon$.

In this way, from (50) we obtain:

[TABLE]

where $C(\delta_{1},\delta_{2},\delta_{3})$ is independent of $\varepsilon$. We note that, using hypothesis (H6) and explicit calculation of $\mathcal{E}_{\varepsilon}(S)$ given in Appendix, we easily see that

[TABLE]

with $C$ independent of $\varepsilon$. After taking the supremum with the respect to $t\in[0,T]$ and by choosing appropriate $\delta_{i}$, we get

[TABLE]

Test-function 2. Another choice for the test-function is given by

[TABLE]

In this way, from (29) we get:

[TABLE]

where \displaystyle I_{5}=\int_{0}^{T}\int_{\partial_{N}\Omega}\mathcal{R}\big{(}\partial_{t}\beta_{\varepsilon}(S_{\varepsilon})-\tau_{\varepsilon}(S_{\varepsilon})\partial_{t}S_{D}\big{)}d\Gamma dt.

We start by estimating the boundary-flux term given by $I_{5}$. We note that the usage of the trace theorem would give the unpleasant term $\displaystyle\int_{0}^{T}\int_{\Omega}|\nabla\partial_{t}\beta_{\varepsilon}(S_{\varepsilon})|^{2}dxdt$ which would make troubles. For that reason we use assumption (H5). Firstly, we note that

[TABLE]

Now, we get

[TABLE]

Other integrals on the right-hand side in (55) can be estimated standardly, using the boundedness assumptions of the regularized coefficients and the Cauchy inequality. In this way we obtain the following expressions:

[TABLE]

[TABLE]

[TABLE]

Using the Cauchy inequality, we get

[TABLE]

and

[TABLE]

We continue with estimation of the integrals on the right-hand side in (55). We have:

[TABLE]

Again, using the Cauchy inequality we get

[TABLE]

and

[TABLE]

By the boundedness assumption of the regularized coefficient $\tau_{\varepsilon}$ and after choosing appropriate constants $\delta_{i}$ we get:

[TABLE]

where the last constant comes from the estimate for the boundary-flux. Finally, we have:

[TABLE]

We now combine estimate (57) with previously obtained estimate (53) to get

[TABLE]

from where it directly follows that

[TABLE]

The last inequality gives estimates (45), (46) and (47). Furthermore, from (57) it follows

[TABLE]

In this way we get estimate (44) and also

[TABLE]

Finally, estimate (43) follows directly from (58) and the fact that, using assumption (26), one has

[TABLE]

In this way we conclude the proof of Theorem 3. $\hfill\Box$

7. A priori estimate of the mixed-derivative term

The main result of this section is given by Proposition 3 where we estimate the third-order term $\displaystyle\nabla(\partial_{t}\beta_{\varepsilon}(S_{\varepsilon}))$. For that we need two technical results given in Proposition 2 and Lemma 2 stated below. The needed integrability will be achieved by combining between results (68) and (69) with (59) and (60) given by the following proposition.

Proposition 2**.**

Under hypotheses (H1)–(H6)* there is a constant $C>0$, independent of $\varepsilon$, and constants $\gamma>2$ and $\lambda>2$ such that every solution $S_{\varepsilon}$ of the regularized problem (29) satisfies*

[TABLE]

Proof. The result of this proposition follows directly from estimate (47) and Lemma 2 given below.

Lemma 2**.**

For $\gamma>2$ and $\lambda>2$ there exist positive constants $C$ and $D$, independent of $\varepsilon$, such that for any $s\in{\mathbb{R}}$

[TABLE]

The proof of Lemma 2 is technical and it is given in Appendix.

Now we state the result which gives the a priori estimate of the third-order term.

Proposition 3**.**

Let $n\leq 3$ and let us suppose hypotheses (H1)–(H6). Furthermore, let the exponents $\beta_{1}$, $\gamma$, $\mu$, and $\lambda$ satisfy:

[TABLE]

where $\sigma=(p-1)/p$, for arbitrary $p>2$. Then there exist constant $C$, independent of $\varepsilon$, such that any solution $S_{\varepsilon}$ of the regularized problem (29) satisfies

[TABLE]

Proof. We start with the estimate

[TABLE]

for $1<m<2$, where we used the Hölder inequality with $r=2/m$ and $r^{\prime}=2/(2-m)$. Due to a priori estimate (44), the first integral on the right-hand side of (66) is bounded. We need to conclude that

[TABLE]

and uniformly bounded with respect to $\varepsilon$. Let us set $\delta=m/(2-m)$ and note that $1<\delta<+\infty$. We consider the integral

[TABLE]

Now we recall the regularization of the capillary diffusion coefficient $a$ given by (24), i.e. $\displaystyle a_{\varepsilon}(s)=a(Z_{\varepsilon}(s))$ where $\displaystyle Z_{\varepsilon}(s)=(1-2\varepsilon)Z(s)+\varepsilon$. Using assumption (H1), we get

[TABLE]

It follows that it is sufficient to undertake to estimate the integrals:

[TABLE]

In order to get the estimates of the mentioned integrals, we will use the informations about the capillary pressure given by assumption (H3). Moreover, from (H3) we have that $\mu\geq\beta_{1}$ and $\lambda\geq\beta_{2}$. Next, we use a priori estimate (46) and we have

[TABLE]

For simplicity we denote $y=Z_{\varepsilon}(S_{\varepsilon})$. Using assumptions (H2) and (H3), expression (67) can be written as follows

[TABLE]

Since \displaystyle\Big{(}y^{\mu}+(1-y)^{\lambda}\Big{)}^{-1}\geq C we get

[TABLE]

from where, using the boundedness from below of functions $g$ and $h$, we get in particular the following estimates:

1. (i)

$\displaystyle\iint_{Q_{T}}y^{\mu}(1-y)^{-\beta_{2}}|\nabla y|^{2}dxdt\leq C$, 2. (ii)

$\displaystyle\iint_{Q_{T}}y^{\mu-\gamma-\beta_{1}}(1-y)^{\lambda}|\nabla y|^{2}dxdt\leq C$.

We start with expression (ii). By defining $\widehat{y}=\min\{y,1/2\}$ from (ii) it follows

[TABLE]

so we get

[TABLE]

with $\alpha_{1}=1+(\mu-\gamma-\beta_{1})/2$. In this way we conclude that

[TABLE]

is uniformly bounded. Because of the Sobolev embedding theorem, one has uniform boundedness in the following spaces:

[TABLE]

In similar way, from expression (i) we get uniform boundedness in the following spaces:

[TABLE]

where $\alpha_{2}=1-\beta_{2}/2$.

In the following part of the proof of Proposition 3, using (68) and (69) we are going to prove that for some $\delta>1$ the term $(a_{\varepsilon}(S_{\varepsilon}))^{-\delta}$ is uniformly bounded in $L^{1}(Q_{T})$. For simplicity, we continue to use notation $y=Z_{\varepsilon}(S_{\varepsilon})$. Again, we will consider separately cases $n=1,2$ and $3$.

Let us start with the case $n=3$. We rewrite the expression $\int_{\Omega}y^{-\mu\delta}dx$ using $-\mu\delta=\alpha_{1}\Theta+(2-\gamma)\Theta_{1}$, and then after the Hölder inequality one has

[TABLE]

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

[TABLE]

Because of (68), the first integral on the right-hand side is uniformly bounded. On the other side due to (59), the second integral is bounded as well. From the condition that

[TABLE]

we get

[TABLE]

Due to assumption (H3), we have

[TABLE]

Using $\mu>\gamma$, from (71) we have $-10+3\beta_{1}>\gamma$, while $0<\beta_{1}\leq\gamma<\mu$ gives $5\gamma>3\beta_{1}+10$. By adding two obtained inequalities it follows that $\gamma>5$. Furthermore, since $\gamma<\mu$, one has $\mu>5\;\textrm{ and }\;\beta_{1}>5$. Finally, for $n=3$ we obtain the following sequence of inequalities:

[TABLE]

which gives (70).

For $n=2$, we write again $-\delta\mu=\Theta\alpha_{1}+(2-\gamma)\Theta_{1}$. Next, we take $\Theta=2$ and $\Theta_{1}=1-2/p$ and we apply the Hölder inequality with $p_{1}=p/2$ and $p_{2}=p/(p-2)$, where $p>2$ is arbitrary finite. In this way we get:

[TABLE]

Now, taking into account the expression for $\alpha_{1}$, from $\delta>1$ we get the condition

[TABLE]

In this way we obtain the following inequality

[TABLE]

with $\gamma>2$ and $0<\beta_{1}\leq\gamma<\mu$. Next, we denote $\nu=(p-1)/p$ and from (73) we get

[TABLE]

Moreover, the upper bound in (74) is the lowest for $\gamma=\beta_{1}$, so from $\beta_{1}<\nu\beta_{1}+(\beta_{1}-4\nu)/2$, we get $\beta_{1}>4\nu/(2\nu-1)=:\tilde{f}(\nu)$. Since $\tilde{f}(\nu)$ is nonincreasing function for $\nu>1/2$, it is sufficient to have $\beta_{1}>\tilde{f}(1)=4$.

For $n=1$, we have:

[TABLE]

where we chose $\Theta=2$ and $\Theta_{1}=1$. The above calculation gives $-\delta\mu=2\alpha_{1}+(2-\gamma)$, and consequently $\delta=-(4+\mu-2\gamma-\beta_{1})/\mu>1$, from where it follows $\beta_{1}>4+2(\mu-\gamma)$, which gives $\beta_{1}>4$. By taking into account assumption (H3), we conclude that for $n=1$ one has

[TABLE]

Next, we consider the term $\displaystyle(1-Z_{\varepsilon}(S_{\varepsilon}))^{-\lambda\delta}$. Similarly like above, using the Hölder inequality and estimates (60) and (69), for $n=1,2$ and $3$ we get the expressions for parameters $\beta_{2}$ and $\lambda$ as it is given by (62)–(63)–(64).

Finally, we conclude that $a_{\varepsilon}^{-m/(2-m)}\in L^{1}(Q_{T})$ for $m=m(n)$. This ends the proof of Proposition 3. $\hfill\Box$

8. End of proof of Theorem 1

For conducting the passage to the limit in (29) when $\varepsilon\to 0$, we will need almost everywhere convergence of the sequence $S_{\varepsilon}$ towards $S$, as $\varepsilon\to 0$. For obtaining this convergence we need results given by the following proposition which, roughly speaking, comes from estimate (47) using the results given in the proof of Lemma 2.

Proposition 4**.**

Under hypotheses (H1)–(H6)* there exist a constant $C>0$, independent of $\varepsilon$, such that every solution $S_{\varepsilon}$ for the regularized problem (29), satisfies:*

[TABLE]

where $|\cdot|$ represents the measure on $Q_{T}$.

The proof of Proposition 4 is based on the proof of Lemma 2 and it is given in the Appendix, as well.

From derived a priori estimates (43)-(47) we conclude that there exists a subsequence of $S_{\varepsilon}$, denoted by the same subscript, and $q\in H^{1}(Q_{T})$ such that, as $\varepsilon\to 0$ one has:

[TABLE]

Because of the compact embedding $H^{1}(Q_{T})\hookrightarrow L^{2}(Q_{T})$ it follows

[TABLE]

Let us write $q=q^{+}+q^{-}$. We are going to show that $q^{-}=0$. Using the definition of $\beta$ given in (16), its extension to whole ${\mathbb{R}}$ and the mean-value theorem, direct calculation gives

[TABLE]

From (78) and Proposition 4 we get

[TABLE]

Next, we write \displaystyle\beta_{\varepsilon}(S_{\varepsilon})=\big{(}\beta_{\varepsilon}(S_{\varepsilon})\big{)}^{+}+\big{(}\beta_{\varepsilon}(S_{\varepsilon})\big{)}^{-}. It is obvious that

[TABLE]

so we conclude that $q\geq 0$ a.e. on $Q_{T}$. On the other side, one has

[TABLE]

We define $\displaystyle\beta^{-1}(q)=:S$ and we conclude that $S\geq 0$. In this way we get that

[TABLE]

Using convergence (79) we also have that $\displaystyle(S_{\varepsilon}-1)^{+}\to(S-1)^{+}=0\;\textrm{ a.e. on }\;Q_{T}$, where the last equality comes directly from (76). In this way we get $S\leq 1$.

Finally, concerning the passage to the limit in (29) when $\varepsilon\to 0$ we comment only the integral $\displaystyle\int_{0}^{T}\int_{\Omega}a_{\varepsilon}(S_{\varepsilon})\nabla\partial_{t}\beta_{\varepsilon}(S_{\varepsilon})\cdot\nabla\varphi dxdt$. From (44) we conclude that $\displaystyle\sqrt{a_{\varepsilon}(S_{\varepsilon})}\nabla\partial_{t}\beta_{\varepsilon}(S_{\varepsilon})\rightharpoonup k$ weakly in $L^{2}(Q_{T})$, as $\varepsilon\to 0$. Furthermore, from (65) and the fact that $\beta_{\varepsilon}(S_{\varepsilon})\to\beta(S)$ strongly in $L^{2}(Q_{T})$, we get $\displaystyle\nabla\partial_{t}\beta_{\varepsilon}(S_{\varepsilon})\rightharpoonup\nabla\partial_{t}\beta(S)$ weakly in $L^{m_{0}}(Q_{T})$, for $m_{0}\in(1,2)$ given by Proposition 3. Because of (79), we get $\sqrt{a_{\varepsilon}(S_{\varepsilon})}\to\sqrt{a(S)}\;\textrm{ a.e.}$ when $\varepsilon\to 0$. Now, the Lebesgue theorem on dominated convergence for $L^{p}$ spaces gives

[TABLE]

The passage to the limit when $\varepsilon\to 0$ of the integral $\displaystyle\int_{0}^{T}\int_{\Omega}a_{\varepsilon}(S_{\varepsilon})\nabla\partial_{t}\beta_{\varepsilon}(S_{\varepsilon})\cdot\nabla\varphi dxdt$ now follows as a product of weak and strong converging sequences. The passage to the limit in other terms in (29) is straightforward. This ends the proof of Theorem 1. $\hfill\Box$

9. Appendix

Proof of Lemma 2. Firstly, we recall the definition of the entropy functional $\mathcal{E}(s)$ given by expression (5) with assumptions on the coefficients $a(s)$ and $\tau(s)$ given by (H1) and (H2). Let us firstly introduce the notation:

[TABLE]

where $\displaystyle h_{1}^{\gamma}(s)=\frac{1}{s^{\gamma}}$ and $\displaystyle h_{2}^{\lambda}(s)=\frac{1}{(1-s)^{\lambda}}$ with $\gamma>0,\lambda>0$.

In Section 4 we regularized the expression for $h(s)$ in the following way:

[TABLE]

where $Z_{\varepsilon}(s)$ was given by

[TABLE]

Regularized entropy functional reads

[TABLE]

where

[TABLE]

Direct calculation gives the following cases:

Case 1: $0\leq s\leq 1$

[TABLE]

Case 2: $s<0$

[TABLE]

Using $s/(Z_{\varepsilon}(S_{D}))^{\gamma-1}\geq s/\varepsilon^{\gamma-1}$ we can estimate

[TABLE]

[TABLE]

Case 3: $s>1$

[TABLE]

We note that $\mathcal{E}_{2,\varepsilon}^{\lambda}(s)$ can be transformed in the following form:

[TABLE]

where we neglected the first term on the right-hand side since it is positive. From derived estimates on $\mathcal{E}_{1,\varepsilon}^{\gamma}(s)$ and $\mathcal{E}_{2,\varepsilon}^{\lambda}(s)$ Lemma 2 follows directly. $\hfill\Box$

Proof of Proposition 4. We start the proof by rewriting the result (47), i.e. we have the following estimate for a.e. $t\in(0,T)$:

[TABLE]

From previous estimates on $\mathcal{E}_{1,\varepsilon}^{\gamma}(s)$ and $\mathcal{E}_{2,\varepsilon}^{\lambda}(s)$ we get

[TABLE]

and

[TABLE]

This completes the proof of Proposition 4. $\hfill\Box$

Acknowledgement

The author acknowledges support from the Croatian Science Foundation HRZZ-IP-2013-11-3955, and the Austrian-Croatian Project of the Austrian Exchange Service (ÖAD) and the Ministry of Science and Education of the Republic of Croatia (MZOS).

This work was partially done during the author was visiting Laboratoire de Mathématiques et de leurs Applications at the University of Pau and Pays de l’Adour. The author warmly thanks Brahim Amaziane not only for hospitality but also for introducing her to the interesting problems concerning the nonequilibrium two-phase flow through porous media.