Electronic Journal of Differential Equations, Vol. 2022 (2022), No. 15, pp. 1–33.

ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu

A NONLINEAR MATHEMATICAL MODEL FOR

TWO-PHASE FLOW IN NANOPOROUS MEDIA

IMANE MELZI, YOUCEF ATIK

Abstract. We propose a mathematical model for the two-phase flow nano-

porous media. Unlike classical models, our model suppose that the rock per-
meability depends on the gradient of pressure. Using usual laws of flows in

porous media, we obtain a system of two nonlinear partial differential equa-

tions: the first is elliptic and the second is parabolic degenerate. We study a
regularized version of our model, obtained by adding a “vanishing” term to the

coefficient causing the degeneracy. We prove the existence of a weak solution

of the regularized model. Our approach consists essentially to use the Rothe’s
method coupled with Galerkin’s method.

1. Introduction

Modeling flow (of shale gas for instance) in nanoporous rocks is becoming an in-
teresting and challenging point for many researchers. A nanoporous media is char-
acterized by an extremely low permeability on the order of a nanodarcy (≈ 10−21

m2) or less. During the exploitation of those kind of porous medium (rocks), there
appears very large pressure gradient at the boundaries of pores causing their ex-
tension or completely their destruction, this phenomena generates a big increase of
the rock permeability. In 2012, Barenblatt et al. [13] proposed a one dimensional
mathematical model describing fluid and gas flow in nanoporous media using a new
formulation of permeability of the rock supposing that it depends on the pressure
gradient (see also [4]). Inspired by the previous work, we propose a three dimen-
sional mathematical model for two-phase flow in nanoporous media. Supposing the
rock permeability depending on the gradient of pressure, using mass conservation,
Darcy’s law, capillary pressure, introducing the concept of global pressure, some
functional coefficients (mobilities, fractional fluxes) and using total velocity u of the
phases; we obtain the following system describing the flow of two incompressible,
immiscible fluids in nanoporous media:

−div
(
λ(s)K(∇p)∇p

)
= q,

φ
∂s

∂t
− div

(
λw(s)K(∇p)∇p+ Λ(s)p′c(s)K(∇p)∇s

)
= qw,

(1.1)

2020 Mathematics Subject Classification. 35M32, 76T99, 74F10.
Key words and phrases. Nonlinear system; nanoporous media; Rothe’s method;

Galerkin’s method.
©2022. This work is licensed under a CC BY 4.0 license.

Submitted January 31, 2021. Published February 28, 2022.

1



2 I. MELZI, Y. ATIK EJDE-2022/15

where s, the saturation of the wetting phase, and p, the global pressure, are the
unknowns. This system is degenerate because the coefficient Λ(s) vanishes for s = 0
and s = 1.

In this article, we prove the existence of weak solution for a regularized version
of the above system to which we associate boundary and initial conditions.

This work is organized as follows: In this section, we complete this Introduction
by presenting the derivation of our model. In Sections 2 and 3, we precise the
hypotheses on the data, regularize the system by adding a term guarantying the
coerciveness of the parabolic equation, extend outside [0, 1] the functional coeffi-
cient depending on s and give the definition of a weak solution of the regularized
system. In Section 4, we discretize in time our system and give the Definition of its
weak discrete time solution. In Section 5 we give Galerkin’s approximations of this
weak solution and prove its existence using a monotonicity method for the pressure
and Brouwder Fixed Point Theorem for the saturation. Section 6 is devoted to
uniform estimates that allow us to pass to the limit on Galerkin’s approximations
in Section 7. We give in Section 8, different uniform estimates on discrete time
solutions which permit us to prove their compactness in Section 9 and to pass to
the limit in Section 10, making the step time goes to zero to obtain our main re-
sult, Theorem 3.2. This work finishes by Section 11, proving a maximum principle
showing that the solution s obtained is a “true” saturation.

1.1. Flow equations. The mass balance equation for each of the fluid phases is

φ(x)
∂(ραsα)

∂t
+ div(ραuα) = ραqα, α = w, n, (1.2)

where α = w denotes the wetting phase (e.g. water), α = n indicates the non
wetting phase (e.g. oil or air), φ is the porosity of the medium Ω which depends
only on x; ρα, sα,uα and qα are respectively the density, (reduced) saturation,
volumetric velocity and external volumetric flow of the α phase.

The Darcy-Muskat’s law is

uα = −KKrα

µα
(∇pα − ραg), α = w, n, (1.3)

where K is the absolute permeability (of the nanoporous medium), pα, µα and Krα

are the pressure, the viscosity and relative permeability of the α phase, respectively.
Several discussions with petroleum engineers and porous media specialists show

us that assuming the absolute permeability a function of pressure gradient seems
to be a good choice.

In this work, we suppose that the absolute permeability is a function of the
pressure gradient (of the wetting phase pw), more precisely, in order to control
that dependency (to control the deformation at the edge of pores), we adopt the
following new formulation of the rock permeability

K(∇pw) = κ1
|∇pw|

1 + η|∇pw|
+ κ2 (1.4)

with κ1 > 0, κ2 > 0, η > 0 three constants. Here η is a positive control constant.
The constant κ2 ensures the coerciveness of our model. Concerning the choice of
K, see the Remark 11.2 at the end of this paper.



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 3

In addition to the above equations, we suppose the customary property of satu-
rations

sw + sn = 1, (1.5)

and introduce the capillary pressure function

pn − pw = pc. (1.6)

To separate the pressure and saturation equations, we introduce the phase mo-
bility functions

λα(x, sα) =
Krα(x, sα)

µα
, α = w, n,

and the total mobility

λ(x, sw) = λw(x, sw) + λn(x, sw).

The fractional flow functions are defined by

fα(x, sw) =
λα(x, sw)

λ(x, sw)
, α = w, n;

finally, we define the total velocity

u = uw + un. (1.7)

In what follows, we re-write the equations in term of primary variables, the
total velocity u, the pressure of wetting phase pw and the saturation of the wetting
phase sw. Under the assumptions that fluids are incompressible (ρα is constant),
summing up equations (1.2) for α = w, n, we obtain

φ
∂

∂t
(sw + sn) + div(uw + un) = qw + qn,

using (1.5) and (1.7), we obtain

div u = q = qw + qn. (1.8)

Also, concerning the total velocity, we have

u = uw + un

= −K(∇pw)λw(sw)(∇pw − ρwg)−Kw(∇p)λn(sw)(∇pn − ρng)

= −K(∇pw)λ(sw)
[λw(sw)

λ(sw)
(∇pw − ρwg) +

λn(sw)

λ(sw)
(∇pn − ρng)

]
,

since pc = pn − pw, we have ∇pn = ∇pw +∇pc and

u = −K(∇pw)λ(sw)
[λw(sw)

λ(sw)
(∇pw − ρwg) +

λn(sw)

λ(sw)
(∇pn − ρng)

]
= −K(∇pw)λ(sw)

[λw(sw)

λ(sw)
(∇pw − ρwg) +

λn(sw)

λ(sw)
(∇pw +∇pc − ρng)

]
= −K(∇pw)λ(sw)

[
∇pw + fn(sw)∇pc − g{fw(sw)ρw + fn(sw)ρn}

]
,

as a result, we have the equation

u = −K(∇pw)λ(sw)
[
∇pw + fn(sw)∇pc − g{fw(sw)ρw + fn(sw)ρn}

]
. (1.9)



4 I. MELZI, Y. ATIK EJDE-2022/15

Remark 1.1. We can obtain expressions

uw = −K(∇pw)λw(sw)(∇pw − ρwg), (1.10)

un = −K(∇pw)λn(sw)(∇pw +∇pc − ρng). (1.11)

So finally, we obtain the system of equations

div u = qw + qn = q,

u = −K(∇pw)λ(sw)
[
∇pw + fn(sw)∇pc − g{fw(sw)ρw + fn(sw)ρn}

]
,

qw = φ
∂sw
∂t
− div

(
K(∇pw)λw(sw)(∇pw − ρwg)

)
,

(1.12)

where the primary unknowns are pw, sw and u. Taking g = 0, the system is written
as

− div
(
λ(sw)K(∇pw)∇pw

)
− div

(
fn(sw)λ(sw)K(∇pw)∇pc

)
= q,

φ
∂sw
∂t
− div

(
λw(sw)K(∇pw)∇pw

)
= qw.

(1.13)

We introduce as in [8], the global pressure

p = pn −
∫ s

0

(
fw
∂pc
∂s

)
(x, ξ) dξ, (1.14)

with s = sw. Making use of the definition of global and capillary pressure, and the
concept of the Differentiation of Integrals (see for example [18, page 213]) we can
write

∇p = ∇pn−∇
∫ s

0

(
fw
∂pc
∂s

)
(x, ξ) dξ = ∇pn−∇s(x)

(
fw
∂pc
∂s

)
(x, s(x))−γ1(x, s(x)),

with

γ1(x, s) =

∫ s

0

∂

∂x

(
fw
∂pc
∂s

)
(x, ξ) dξ.

Now, we have

∇pw = ∇pn −∇pc = ∇p+∇
∫ s

0

fw(ξ)
∂pc
∂s

(x, ξ) dξ −∇pc

= ∇p+ fw(s)∇pc(s) + γ1(x, s)−∇pc
= ∇p+ (fw(s)− 1)∇pc(s) + γ1(x, s)

= ∇p− fn(s)∇pc(s) + γ1(x, s),

so

K(∇pw) = K
(
∇p− fn(s)∇pc(s) + γ1(x, s)

)
.
= K(∇p, s,∇s).

This leads to the system

−div
(
λ(s)K(∇p, s,∇s)(∇p+ γ1(x, s))

)
= q,

φ
∂s

∂t
− div

(
λw(s)K(∇p, s,∇s)∇p−K(∇p, s,∇s)Λ(s)p′c(s)∇s

+ λw(s)K(∇p, s,∇s)γ1(s)
)

= qw,

where Λ(s) = λw(s)λn(s).λ(s). From a theoretical point of view, and in order to
simplify the model, we are going to neglect the term γ1 and assume that

K(∇p, s,∇s) ≈ K(∇p).



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 5

Consequently, we obtain the system

−div
(
λ(s)K(∇p)∇p

)
= q,

φ
∂s

∂t
− div

(
λw(s)K(∇p)∇p+ Λ(s)p′c(s)K(∇p)∇s

)
= qw,

to which we must add initial and boundary conditions. Therefore, we consider the
system:

Find (p, s) solving the equations

−div
(
λ(s)K(∇p)∇p

)
= q in ΩT

.
= Ω×]I[, (1.15)

φ
∂s

∂t
− div

(
λw(s)K(∇p)∇p+ Λ(s)p′c(s)K(∇p)∇s

)
= qw in ΩT , (1.16)

p(x, t) = 0 et s(x, t) = 0 on ∂Ω× [I], (1.17)

s(x, 0) = s0(x) in Ω, (1.18)

where Ω ⊂ R3 represents the nanoporous medium, supposed to be bounded, con-
nected and Lipschitz domain, I

.
=]0, T [ is the time interval, and with the following

expression of the absolute permeability given in page 2:

K(∇p) = κ1
|∇p|

1 + η|∇p|
+ κ2, with κ1 > 0, κ2 > 0, η > 0 three constants,

and q = q(x, t), qw = qw(x, t), s0 = s0(x) three given functions.
In all that follows, we will denote by (S) the system of equations (1.15) and

(1.16), with boundary conditions (1.17) and the initial condition (1.18).

1.2. Functional setting. We denote by V the Sobolev space H1
0 (Ω), equipped

with the inner product (u, v)V =
∫

Ω
∇u · ∇v dx and the gradient norm ‖u‖V =[ ∫

Ω
|∇u|2 dx

]1/2
, its dual is indicated by V ?. For 1 ≤ p < ∞ and B a Banach

space, we denote Lp(I;B) the Bochner space (of classes with respect to equiva-

lence a.e.) of Bochner integrable functions u : I −→ B satisfying
∫ T

0
‖u(t)‖pB dt <

+∞. This space is a Banach space if endowed with the norm ‖u‖Lp(0,T ;B) =( ∫ T
0
‖u(t)‖pB dt

)1/p
. For p = ∞, this norm is ‖u‖L∞(I;B) = ess supt∈I ‖u(t)‖B .

Following [16], we denote W (0, T ), the Sobolev-Bochner space

W (0, T ) = W 2,2(0, T ;V, V ?) =
{
u ∈ L2(0, T, V )

∣∣u′ =
du

dt
∈ L2(0, T ;V ?)

}
·

Equipped by the norm ‖u‖W =
(
‖u‖2L2(I;V ) +‖u′‖2L2(I;V ?)

)1/2
, W (0, T ) is a Hilbert

space which is continuously embedded in C([0, T ];L2(Ω)), equipped with the norm
of uniform convergence. Proofs of the above facts can be found in [10, 14, 16].

2. Hypotheses

(H1) The porosity φ belongs to W 1,+∞(Ω), and for two constants, φ∗ and φ∗ we
have

0 < φ∗ ≤ φ(x) ≤ φ∗ < +∞, a.e. x ∈ Ω. (2.1)

(H2) λα(x, s), α = n,w are measurable in x ∈ Ω and continuous in s ∈ [0, 1],
and satisfies λn(x, 1) = 0, λn(x, s) > 0 for s < 1, λw(x, 0) = 0, λw(x, s) > 0
for s > 0; and there exist two constants λ∗, λ

∗ such that

0 < λ∗ ≤ λ(x, s) ≤ λ∗ < +∞, x ∈ Ω, s ∈ [0, 1], (2.2)



6 I. MELZI, Y. ATIK EJDE-2022/15

where λ(s) = λw(s) + λn(s).
(H3) pc ∈ W 1,∞([0, 1]) and p′c is a continuous function, also there exits two

constants p′c∗, and p′
∗
c such that

0 < p′c∗ ≤ p′c(s) ≤ p′
∗
c< +∞.

(H4) The initial saturation s0 is in L2(Ω), the functions q and qw are positive
functions in L2(ΩT ).

In what follows, we put

Λε(x, s) = Λ(x, s) + ε with Λ(x, s)=
λw(x, s)λn(x, s)

λ(x, s)
and ε > 0.

Remark 2.1. Hypothesis (H1) permits us, among other things, to put 〈φ∂ts, v〉 :=
〈∂ts, φv〉 for v ∈ V = H1

0 (Ω) to give sense to φ∂ts knowing that ∂ts ∈ V ?. This
because φ ∈ W 1,∞(Ω) implies that φv ∈ V , for all v ∈ V . Also, we should note
that, during the entire work, inequality (2.1) is used to obtain different estimations
on the equation of saturation.

3. Regularization: system (Sε)

We extend the coefficients of identities (1.15), (1.16) outside [0, 1] as continuous
functions in s by putting

λw?(x, s) =


λw(x, s), x ∈ Ω, s ∈ [0, 1],

λw(x, 1), x ∈ Ω, s ≥ 1,

λw(x, 0), x ∈ Ω, s ≤ 0,

λn?(x, s) =


λn(x, s), x ∈ Ω, s ∈ [0, 1],

λn(x, 1), x ∈ Ω, s ≥ 1,

λn(x, 0), x ∈ Ω, s ≤ 0.

The capillary pressure pc is extended outside [0, 1] in the same way. Also, we put

λ?(x, s) = λw?(x, s) + λn?(x, s), Λ?(x, s) =
λw?(x, s)λn?(x, s)

λw?(x, s) + λn?(x, s)
,

Λε(x, s) = Λ?(x, s) + ε, ε > 0.

Substituting these functions in the system (S) (equations (1.15)–(1.18)), we obtain
the system (Sε).

3.1. Weak solution of system (Sε).

Definition 3.1. A weak solution of system (Sε) is a couple (p, s) such that

(p, s) ∈ L2(I;V )× L2(I;V ), ∂ts ∈ L2(I;V ?), (3.1)∫
Ω

λ?(s)K(∇p)∇p · ∇ϕdx = (q, ϕ), ∀ϕ ∈ V, a.e. t ∈ I, (3.2)∫ T

0

〈φ∂s
∂t

, ψ
〉
dt+

∫
ΩT

λw?(s)K(∇p)∇p · ∇ψ dx dt

+

∫
ΩT

Λε(s)p
′
c(s)K(∇p)∇s · ∇ψ dx dt

=

∫ T

0

(qw, ψ) dt, ∀ψ ∈ L2(I;V ),

(3.3)

s(x, 0) = s0(x). (3.4)

Theorem 3.2. Under the hypothesis (H1)–(H4), Problem (Sε) has at least one
weak solution in the sense of Definition 3.1.



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 7

For simplicity, we omit subindex ?, standing for extension of λ, λn, and λw
outside [0, 1].

4. Time discretization of system (Sε)

To show the existence of a weak solution of the system (Sε) in the sense of
definition 3.1, we use the method of Rothe (semi-discretization in time) coupled
with Galerkin’s method. To do this, for each positive integer n, we divide the
interval I = ]0, T [ into N = 2n subintervals and we set α = T

N = 2−nT and put
tj = jα and Ij = (tj−1, tj ] for any integer j, j = 1, . . . , N . We approach the time

derivative ∂s
∂t by the time difference operator

∂αt s(x, t) =
s(x, t+ α)− s(x, t)

α
.

If w = w(x, t) is a function, the average in time over Ij is

wα(x, t) =
1

α

∫
Ij

w(x, τ) dτ, t ∈ Ij . (4.1)

The value of wα(·) on the interval Ij is denoted by wαj(·). Also, for any linear
space H, we define

`α(I;H) = {v ∈ L∞(I;H) : v is constant in time on each subinterval Ij ⊂ I}.
The value of a function vα(·) from the space `α(I;H) on the interval Ij is constant
and it is equal to vα(tj)(·) which will be denoted by vjα(·), i.e.

vα(x, t) =

N∑
j=1

vα(x, tj)χ]tj−1,tj ](t) =

N∑
j=1

vjα(x)χ]tj−1,tj ](t).

We define also the function ṽα by

ṽα(x, t) =

N∑
j=1

[vα(x, tj)− vα(x, tj−1)

α
(t− tj−1) + vα(x, tj−1)

]
χ[tj−1,tj [(t)

=

N∑
j=1

[vjα(x)− v(j−1)α(x)

α
(t− tj−1) + v(j−1)α(x)

]
χ[tj−1,tj [(t),

where we put ṽ0α(x) = ṽ(x, 0) = v0(x), a given function supposed hereafter to play
the role of the initial condition.

Remark 4.1. Performing simple calculations, one can easily see that

‖wα‖2L2(ΩT ) = α

N∑
j=1

‖wαj‖2L2(Ω), ‖wα‖L2(ΩT ) ≤ ‖w‖L2(ΩT ),

‖vα‖2L2(I;V ) = α

N∑
j=1

‖∇vjα‖2L2(Ω), ‖ṽα‖2L2(I;V ) ≤ 5α

N∑
j=1

‖∇vjα‖2L2(Ω),

‖vα − ṽα‖2L2(I;X) =
α

3

N∑
j=1

‖vjα − v(j−1)α‖2X ,

∂ṽα

∂t
=

N∑
j=1

vjα − v(j−1)α

α
χ[tj−1,tj [, a.e.,



8 I. MELZI, Y. ATIK EJDE-2022/15

∥∥∂ṽα
∂t

∥∥2

L2(I;V )
=

1

α

N∑
j=1

‖vjα − v(j−1)α‖2V .

Definition 4.2. A discrete time solution is a couple of functions

(pα, sα) ∈ `α(I;V )× `α(I;V )

which satisfies ∫
Ω

λ(s(j−1)α)K(∇pjα)∇pjα · ∇ϕdx = (qαj , ϕ),

∀ϕ ∈ V, t ∈ Ij , j = 1, . . . , N,
(4.2)

∫ T

0

(φ∂−αt sα, ψ) dt+

∫
ΩT

λw(sα)K(∇pα)∇pα · ∇ψ dx dt

+

∫
ΩT

Λε(s
α)p′c(s

α)K(∇pα)∇sα · ∇ψ dx dt

=

∫ T

0

(qwα, ψ) dt, ∀ψ ∈ `α(I;V ).

(4.3)

Regarding the first term in (4.3), we have
∫ T

0
(φ∂−αt sα, ψ) dt =

∫ T
0
〈∂s̃

α

∂t , ψ〉 dt this

because ∂−αt sα = ∂s̃α

∂t . In fact,

∂−αt sα(x, t) =
sα(x, t)− sα(x, t− α)

α

=

∑N
j=1(sα(x, tj)− sα(x, tj − α)χ[tj−1,tj [(t)

α

=

∑N
j=1(sjα(x)− s(j−1)α(x))χ[tj−1,tj [(t)

α

=
∂s̃α

∂t
(x, t).

Let us re-write the integral identity (4.3) in an equivalent form. By taking the
test function in the form χIj (t)ϕ(x), where χIj is the characteristic function of the
interval Ij = [tj−1, tj [= [(j − 1)α, jα[= [j′α, jα[, and ϕ is a function in the space
V , we then obtain∫

Ij

(
φ
sα(t)− sα(t− α)

α
, ϕ
)
dt+

∫
Ij

(
λw(sα)K(∇pα)∇pα,∇ϕ

)
dt

+

∫
Ij

(
Λε(s

α)p′c(s
α)K(∇pα)∇sα,∇ϕ

)
dt

=

∫
Ij

(qwα, ϕ) dt.

Since sα(·, t) is constant with respect to t on the interval Ij and it is equal to
sα(·, tj), the same thing is true for pα(·, t), so, we obtain the following integral
identity (

φsjα, ϕ
)

+ α
(
λw(sjα)K(∇pjα)∇pjα,∇ϕ

)
+ α

(
Λε(s

jα)p′c(s
jα)K(∇pjα)∇sjα,∇ϕ

)
=
(
φsj

′α, ϕ
)

+ α(qwαj , ϕ), ∀ϕ ∈ V.

(4.4)



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 9

5. Galerkin’s approximations of discrete time solutions

To use the Galerkin procedure in determining the solution at the level t1 = α,
we choose an orthonormal basis {ei}∞i=1 in V . Therefore the subspaces Hd =
〈e1, · · · , ed〉, d ∈ N, spanned by these functions are denses in V , and then we look
for functions, written as

p1α(·) ∈ V and s1α
d (x) =

d∑
i=1

σ1
i ei(x),

where {σ1
i }di=1 are unknowns real coefficients, and satisfying, for all ϕ ∈ V ,∫

Ω

λ(s0)K(∇p1α)∇p1α · ∇ϕdx = (qα1, ϕ), (5.1)

and, for all ψ ∈ Hd,(
φs1α

d , ψ
)

+ α
(
λw(s1α

d )K(∇p1α)∇p1α,∇ψ
)

+ ακ1

(
Λε(s

1α
d )p′c(s

1α
d )K(∇p1α)∇s1α

d ,∇ψ
)

=
(
φs0, ψ

)
+ α(qwα1, ψ).

(5.2)

To be brief, instead of (s1α
d ), we denote (sd).

5.1. Existence of Galerkin’s approximations. Before giving the proof of exis-
tence, we give the following statement.

Remark 5.1. The finite dimensional spaceHd is equipped with the three equivalent

norms defined, for v =
∑d
k=1 αkek ∈ Hd, by

|v|Rd =
[ d∑
k=1

α2
k

]1/2
, |v|2 =

[ ∫
Ω

|v|2 dx
]1/2

, ‖v‖V =
[ ∫

Ω

|∇v|2 dx
]1/2

.

Let us explain the first step of the existence of Galerkin’s approximation. In the
beginning, we have to find p1α solution of

κ1

∫
Ω

λ(s0)
|∇p1α|

1 + η|∇p1α|
∇p1α · ∇ϕdx+ κ2

∫
Ω

λ(s0)∇p1α · ∇ϕdx = (qα1, ϕ),

for all ϕ ∈ V , which is equivalent to∫
Ω

λ(s0)K12(∇p1α) · ∇ϕdx = (qα1, ϕ), ∀ϕ ∈ V,

with K12(x)
.
= K1(x) +K2(x) where (| · | stands for the Euclidean norm)

R3 3 x 7−→ K1(x) = κ1
|x|x

1 + η|x|
∈ R3, R3 3 x 7−→ K2(x) = κ2x ∈ R3. (5.3)

Consider now the operator A from V = H1
0 (Ω) into its dual V ? = H−1(Ω), given

by

〈A(p), v〉 =

∫
Ω

λ(s0)K12(∇p) · ∇v dx, ∀v ∈ V.

Notice that in fact A(p) ∈ L2(Ω), for all p ∈ V .
We have the following results:



10 I. MELZI, Y. ATIK EJDE-2022/15

(1) Operator A is coercive in the sense that 〈A(p),p〉
‖p‖V → +∞ when ‖p‖V → +∞.

This because, K12x · x ≥ κ2|x|2, for all x ∈ R3, and therefore

〈A(p), p〉
‖p‖

=

∫
Ω
λ(s0)K12(∇p) · ∇p dx
‖∇p‖L2(Ω)

≥
κ2λ∗‖∇p‖2L2

‖∇p‖L2(Ω)
= κ2λ∗‖p‖V .

(2) Operator A is monotone, i.e., for all p, q ∈ V : 〈A(p)−A(q), p− q〉 ≥ 0. In fact,
for p, q ∈ V , we have

〈A(p)−A(q), p− q〉
= 〈A(p), p− q〉 − 〈A(q), p− q〉

=

∫
Ω

λ(s0)[K12(∇p)−K12(∇q)] · ∇(p− q) dx

=

∫
Ω

λ(s0)(K1(∇p)−K1(∇q)) · ∇(p− q) dx+ κ2

∫
Ω

λ(s0)|∇(p− q)|2 dx.

It is easy to see that

(K1(∇p)−K1(∇q)) · ∇(p− q) (Cauchy-Schwarz)

= κ1

{ |∇p|∇p
1 + η|∇p|

− |∇q|∇q
1 + η|∇q|

}
· ∇(p− q)

≥ κ1

{ |∇p|3

1 + η|∇p|
− |∇p|

2|∇q|
1 + η|∇p|

+
|∇q|3

1 + η|∇q|
− |∇p‖∇q|

2

1 + η|∇q|

}
= κ1(|∇p| − |∇q|)

{ |∇p|2

1 + η|∇p|
− |∇q|2

1 + η|∇q|

}
.

Let us now consider the real function R+ 3 ξ 7−→ f(ξ) = ξ2

1+ηξ ∈ R+. We have

f ′(ξ) =
2ξ + ηξ2

(1 + ηξ)2
> 0, ∀ξ > 0.

Putting ξ = |∇p| and σ = |∇q|, we see, by using the Mean Value Theorem, that

|∇p|2

1 + η|∇p|
− |∇q|2

1 + η|∇q|
= f(ξ)− f(σ) = (ξ − σ)f ′(cξη)

where cξη is a point between ξ and η. We conclude that

(K1(∇p)−K1(∇q)) · ∇(p− q) ≥ κ1(|∇p| − |∇q|)2f ′(cξη) ≥ 0, ∀p, q ∈ V.
This implies that

〈A(p)−A(q), p− q〉 ≥ κ2λ∗‖p− q‖2V , ∀p, q ∈ V,
showing that A is in fact strongly monotone, see, for instance [12] or [14].

(3) A is bounded. Let p ∈ V with ‖p‖V ≤ M , we have ‖A(p)‖V ? ≤ M ′. In fact, if
p ∈ V with ‖p‖V ≤M , we have

〈A(p), p〉 =

∫
Ω

λ(s0)K12(∇p) · ∇p dx

≤
(κ1

η
+ κ2

)
λ∗‖∇p‖2L2(Ω)

=
(κ1

η
+ κ2

)
λ∗‖p‖2V ≤

(κ1

η
+ κ2

)
λ∗M2 = C,



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 11

and

‖A(p)‖V ? = sup
q∈V
‖q‖V ≤1

|〈A(p), q〉|

≤ sup
q∈V
‖q‖≤1

(κ1

η
+ κ2

)∫
Ω

λ(s0)|∇p| |∇q| dx

≤ sup
q∈V
‖q‖≤1

λ∗
(κ1

η
+ κ2

)
‖∇p‖L2(Ω) ‖∇q‖L2(Ω)

≤ λ∗
(κ1

η
+ κ2

)
‖∇p‖L2(Ω) ≤ λ∗

(κ1

η
+ κ2

)
M = M ′.

Operator A is hemicontinuous. Let p, q, r ∈ V be three functions. Let us prove
that the application defined from R into R by θ 7→ 〈A(p+ θq), r〉 is continuous. To
see that, we consider a real sequence (θn)n converging to θ. First, because of the
continuity of the function x 7−→ K12(x), we have

K12(∇p+ θn∇q) · ∇r −−−−−→
a.e. x∈Ω

K12(∇p+ θ∇q) · ∇r.

Second, since the sequence (θn)n is convergent, there exists a constant M > 0 such
that |θn| ≤M , ∀n ∈ N. Then, we obtain∣∣∣K12(∇p+ θn∇q) · ∇r

∣∣∣ ≤ (κ1

η
+ κ2

)
|∇p+ θn∇q| |∇r|

≤
(κ1

η
+ κ2

)(
|∇p|+M |∇q|

)
|∇r| ∈ L1(Ω).

Using the Lebesgue’s Dominated Convergence Theorem, we see that∫
Ω

λ(s0)K12(∇p+ θn∇q) · ∇r dx −−−−−→
n→+∞

∫
Ω

λ(s0)K12(∇p+ θ∇q) · ∇r dx.

This means 〈A(p+ θnq), r〉 −−−−−→
n→+∞

〈A(p+ θq), r〉, which is the hemicontinuity of A.

As a result, the operator A is bounded, hemicontinuous, monotone and coercive. Con-
sequently, for qα1 ∈ L2(Ω), there exists p1α solution of (5.1), see [12] or [14].

Now, to prove the existence of sd(= s1α
d ) solution to (5.2), we use a variant of Brouwer’s

Fixed Point Theorem which asserts that a continuous mapping P from Rd into itself
satisfying, for some ρ > 0, P (σ) · σ ≥ 0, for all σ, |σ| = ρ, has at least a zero σ0 ∈ Rd
with |σ0| ≤ ρ, see, for instance, [12, page 53].

Let us therefore consider the operator Rd 3 σ 7→ P (σ) = β ∈ Rd where β
.
= (β1, . . . , βd),

defined, for k = 1, . . . , d, by

βk =

∫
Ω

φ
sσd − s0

α
ek dx+

∫
Ω

λw(sσd )K(∇p1α)∇p1α · ∇ek dx

+

∫
Ω

Λε(s
σ
d )p′c(s

σ
d )K(∇p1α)∇sσd · ∇ek dx− (qwα1, ek).

Here sσd =
∑d
l=1 σlel, for σ = (σ1, . . . , σd).

The operator P has the following properties:
(1) P is continuous. Let {σm}∞m=1

.
= {(σm1 , . . . , σmd )}∞m=1 a sequence in Rd converging

in this space to σ
.
= (σ1, . . . , σd). We have to prove that the sequence {P (σm)}∞m=1 is

converging to P (σ). We do have the following convergences

sσ
m

d (x) =

d∑
l=1

σml el(x)
a.e. x−−−−→
m→∞

d∑
l=1

σlel(x) = sσd (x),



12 I. MELZI, Y. ATIK EJDE-2022/15

∇sσ
m

d (x) =

d∑
l=1

σml ∇el(x)
a.e. x−−−−→
m→∞

d∑
l=1

σl∇el(x) = ∇sσd (x).

Using the continuity of the concerned functions

λw(sσ
m

d (x))
a.e. x−−−−−→
m→∞

λw(sσd (x)), Λε(s
σm

d (x))
a.e. x−−−−−→
m→∞

Λε(s
σ
d (x)),

p′c(s
σm

d (x))
a.e. x−−−−−→
m→∞

p′c(s
σ
d (x)).

Which implies that (we omit the variable x)

λw(sσ
m

d )K(∇p1α)∇p1α · ∇ek
a.e. x−−−−→
m→∞

λw(sσd )K(∇p1α)∇p1α · ∇ek,

Λε(s
σm

d )p′c(s
σm

d )K(∇p1α)∇sσ
m

d · ∇ek
a.e. x−−−−→
m→∞

Λε(s
σ
d )p′c(s

σ
d )K(∇p1α)∇sσd · ∇ek,

We have also the estimate

|λw(s)|+ |Λε(s)|+ |p′c(s)|+ |σm|Rd ≤M, ∀s ∈ R, ∀m ∈ N.

The constant M depends, among other things, on λ∗, λ
∗, p′

∗
c , and ε. Also

|λw(sσ
m

d )K(∇p1α)∇p1α · ∇ek| ≤
(κ1

η
+ κ2

)
M |∇p1α · ∇ek|

≤
(κ1

η
+ κ2

)
M |∇p1α| |∇ek| ∈ L1(Ω),

and

|Λε(sσ
m

d )p′c(s
σm

d )K(∇p1α)∇sσ
m

d · ∇ek| ≤
(κ1

η
+ κ2

)
M2|∇sσ

m

d | |∇ek|

≤
(κ1

η
+ κ2

)
M3

d∑
l=1

|∇ϕl| |∇ek| ∈ L1(Ω),

By the Lebesgue Dominated Convergence Theorem, limm→∞ P (σm) = P (σ), which proves
the continuity of P .

(2) There exists ρ > 0 such that P (σ) · σ ≥ 0 for all σ ∈ Rd with |σ| = ρ. Here the
central dot stands for the classical dot (scalar) product in Rd. We have

P (σ) · σ =

d∑
k=1

σk
(∫

Ω

φ
sσd − s0

α
ek dx+

∫
Ω

λw(sσd )K(∇p1α)∇p1α · ∇ek dx

+

∫
Ω

Λε(s
σ
d )p′c(s

σ
d )K(∇p1α)∇sσd · ∇ek dx− (qwα1, ek)

)
and

P (σ) · σ

=

∫
Ω

φ
sσd − s0

α

( d∑
k=1

σkek
)
dx+

∫
Ω

λw(sσd )K(∇p1α)∇p1α ·
( d∑
k=1

σk∇ek
)
dx

+

∫
Ω

Λε(s
σ
d )p′c(s

σ
d )K(∇p1α)∇sσd ·

( d∑
k=1

σk∇ek
)
dx−

(
qwα1,

( d∑
k=1

σkek
))
,

which is equivalent to

P (σ) · σ =

∫
Ω

φ
sσd − s0

α
sσd dx+

∫
Ω

λw(sσd )K(∇p1α)∇p1α · ∇sσd dx

+

∫
Ω

Λε(s
σ
d )p′c(s

σ
d )K(∇p1α)∇sσd · ∇sσd dx− (qwα1, s

σ
d ).



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 13

Let us estimate each one of the four terms of the above equality. Concerning the first and
second terms, we have the estimates∫

Ω

φ
sσd − s0

α
sσd dx =

∫
Ω

φ

α
(sσd )2 dx−

∫
Ω

φ

α
s0s

σ
d dx

≥ φ∗
α
‖sσd‖2L2(Ω) −

φ∗

α
‖s0‖L2(Ω) ‖s

σ
d‖L2(Ω)

≥ φ∗
α
‖sσd‖2L2(Ω) − Cp

φ∗

α
‖s0‖L2(Ω) ‖∇s

σ
d‖L2(Ω),

where Cp is Poincaré’s constant, and∣∣∣ ∫
Ω

λw(sσd )K(∇p1α)∇p1α · ∇sσd dx
∣∣∣ ≤ λ∗(κ1

η
+ κ2

)
‖∇p1α‖L2(Ω) ‖∇s

σ
d‖L2(Ω).

For the third term, we have the estimate∫
Ω

Λε(s
σ
d )p′c(s

σ
d )K(∇p1α)|∇sσd |2 ≥ εp′c∗κ2‖∇sσd‖2L2(Ω).

Concerning the fourth term,∣∣∣ ∫
Ω

qwα1s
σ
d

∣∣∣ ≤ ‖qwα1‖L2(Ω) ‖s
σ
d‖L2(Ω) ≤ Cp‖qwα1‖L2(Ω) ‖∇s

σ
d‖L2(Ω).

Collecting the previous estimates, we see that

P (σ) · σ ≥ φ∗
α
‖sσd‖2L2(Ω) dx− Cp

φ∗

α
‖s0‖L2(Ω) ‖∇s

σ
d‖L2(Ω)

− λ∗
(κ1

η
+ κ2

)
‖∇p1α‖L2(Ω) ‖∇s

σ
d‖L2(Ω) + εp′c∗κ2‖∇sσd‖2L2(Ω)

− Cp‖qwα1‖L2(Ω) ‖∇s
σ
d‖L2(Ω)

= ‖∇sσd‖L2(Ω)

(
εp′c∗κ2‖∇sσd‖L2(Ω) − Cp

φ∗

α
‖s0‖L2(Ω)

− λ∗
(κ1

η
+ κ2

)
‖∇p1α‖L2(Ω) − Cp‖qwα1‖L2(Ω)

)
+
φ∗
α
‖sσd‖2L2(Ω).

Thus P (σ) · σ ≥ 0 if

‖∇sσd‖ ≥
[
Cp

φ∗

α
‖s0‖L2(Ω) + λ∗

(κ1

η
+ κ2

)
‖∇p1α‖L2(Ω)

+ Cp‖q1
wα‖L2(Ω)

]
{εp′c∗κ2}−1 .

= ρ0.

Let us recall here that {ei}∞i=1 being an orthonormal basis in V , implies that, if sσd =∑d
l=1 σlel, one has |σ|Rd = ‖sσd‖V . Therefore, P (σ) · σ ≥ 0 for all σ ∈ Rd with |σ| = ρ for

all ρ ≥ ρ0.
We are now in a position to apply the variant of Brouwer’s Fixed Point Theorem

mentioned above: there exists σ0 = (σ01, . . . , σ0d) ∈ Rd, with |σ0| ≤ ρ, such that P (σ0) =

0. It is now easy to see that s1α
d = sσ0d =

∑d
l=1 σ0lel is a solution of (5.2).

6. Uniform estimates on Galerkin’s approximations

Proposition 6.1. Let (p1α, s1α
d ) a solution to the system (5.1)–(5.2) at the time level

t1 = α. For the functions (s1α
d )d≥1 the following estimate holds,

‖s1α
d ‖V ≤ C, ∀d ≥ 1, (6.1)

where C is a positive constant independent of d.



14 I. MELZI, Y. ATIK EJDE-2022/15

Proof. By writing (5.2) with ei as the test function, then multiplying by σ1
i and summing

i from 1 to d, we have∫
Ω

φ
sd − s0

α
sd dx+

∫
Ω

λw(sd)K(∇p1α)∇p1α · ∇sd dx

+

∫
Ω

Λε(sd)p
′
c(sd)K(∇p1α)|∇sd|2 dx = (qwα1, sd),

(6.2)

concerning the first term, we have∫
Ω

φ(sd − s0)sd dx =

∫
Ω

φ(sd)
2 dx−

∫
Ω

φs0 × sd dx

with ∫
Ω

φ(sd)
2 dx ≥ φ∗

∫
Ω

(sd)
2 dx = φ∗‖sd‖2L2(Ω)

and by Young’s Inequality, for β > 0, we obtain∣∣∣ ∫
Ω

φs0 sd dx
∣∣∣ ≤ φ∗Cp 1

2β
‖s0‖2L2(Ω) + φ∗Cp

β

2
‖∇sd‖2L2(Ω).

For the second term of (6.2), using Young’s Inequality, with β1 > 0, we obtain∣∣∣ ∫
Ω

λw(sd)K(∇p1α)∇p1α · ∇sd dx
∣∣∣

≤ λ∗
(κ1

η
+ κ2

)
‖∇p1α‖L2(Ω) ‖∇sd‖L2(Ω),

≤ 1

2β1

(
λ∗
[κ1

η
+ κ2

])2

‖∇p1α‖2L2(Ω) +
β1

2
‖∇sd‖2L2(Ω).

Concerning the third term of (6.2), we have∫
Ω

Λε(sd)p
′
c(sd)K(∇p1α)|∇sd|2 dx ≥ εp′c∗κ2‖∇sd‖2L2(Ω),

and finally for the last term and using again Young’s Inequality, for β2 > 0, we have∣∣∣ ∫
Ω

qwα1sd dx
∣∣∣ ≤ ‖qwα1‖L2(Ω) ‖sd‖L2(Ω)

≤ Cp‖qwα1‖L2(Ω) ‖∇sd‖L2(Ω)

≤ 1

2β2
(Cp‖q1

wα‖L2(Ω))
2 +

β2

2
‖∇sd‖2L2(Ω).

By taking into account all the previous estimates, we have

αεp′c∗κ2‖∇sd‖2L2(Ω) + φ∗‖sd‖2L2(Ω)

≤ φ∗Cp
1

2β
‖s0‖2L2(Ω) + φ∗Cp

β

2
‖∇sd‖2L2(Ω) +

α

2β1

(
λ∗
[κ1

η
+ κ2

])2

‖∇p1α‖2L2(Ω)

+
αβ1

2
‖∇sd‖2L2(Ω) +

α

2β2
(Cp‖qwα1‖L2(Ω))

2 +
αβ2

2
‖∇sd‖2L2(Ω),

which implies that(
αεp′c∗κ2 − φ∗Cp

β

2
− αβ1

2
− αβ2

2

)
‖∇sd‖2L2(Ω)

≤ φ∗Cp
1

2β
‖s0‖2L2(Ω) +

α

2β1

(
λ∗
[κ1

η
+ κ2

])2

‖∇p1α‖2L2(Ω) +
α

2β2
(Cp‖qwα1‖L2(Ω))

2 .

Taking

β =
αεp′c∗κ2

2Cpφ∗
, β1 = β2 =

εp′c∗κ2

2
,

we obtain
1

4
αεp′c∗κ2‖∇sd‖2L2(Ω)



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 15

≤ φ∗Cp
1

2β
‖s0‖2L2 +

α

2β1

(
λ∗
[κ1

η
+ κ2

])2

‖∇p1α‖2L2(Ω) +
α

2β2
C2
p‖qwα1‖2L2(Ω),

then estimate (6.1) holds. �

7. Passing to the limit (with respect of d) in Galerkin’s approximations

We proved previously that

‖∇s1α
d ‖L2(Ω) ≤ C, ∀d ≥ 1 .

Since the sequence {s1α
d }∞d=1 is bounded in V (associated with the norm of gradient) we

can extract a subsequence (denoted in the same symbol) such that

s1α
d ⇀ s1α weakly in V and s1α

d → s1α a.e. in Ω, (7.1)

more precisely, by the Rellich-Kondrachov theorem s1α
d → s1α strongly in L2(Ω) and by the

inverse of the Dominated Convergence theorem of Lebesgue, we can extract a subsequence
which converge almost everywhere. Let d0 a positive integer, since the sequence of linear
spaces Hd are nested, we have(

φs1α
d , ψ

)
+ α

(
λw(s1α

d )K(∇p1α)∇p1α,∇ψ
)

+ α
(
Λε(s

1α
d )p′c(s

1α
d )K(∇p1α)∇s1α

d ,∇ψ
)

=
(
φs0, ψ

)
+ α(qwα1, ψ)ΓN , ∀ψ ∈ Hd0 , ∀d ≥ d0.

(7.2)

Let us fix ψ in Hd0 . Using the convergences (7.1) and taking into account that the
functions λw, Λε and p′c are bounded and continuous, by making d goes to infinity in the
equation of saturation (7.2) and using the Dominated Convergence Theorem of Lebesgue,
one obtains (

φs1α, ψ
)

+ α
(
λw(s1α)K(∇p1α)∇p1α,∇ψ

)
+ α

(
Λε(s

1α)p′c(s
1α)K(∇p1α)∇s1α,∇ψ

)
=
(
φs0, ψ

)
+ α(qwα1, ψ)ΓN , ∀ψ ∈ Hd0 .

(7.3)

Now, using the density of ∪∞d=1Hd in V , we see that the previous integral identity is satisfied
for all ψ ∈ V . This makes an end to the proof of existence of the couple (p1α, s1α) solution
of the system (S)ε at the time level t1 = α.

Note that the same reasoning permits us to prove inductively the existence of the
discrete time solution (pjα, sjα) at each time level tj = jα for j = 2, . . . , N .

Knowing the functions pjα, sjα at levels j = 1, . . . , N , we construct the Rothe’s func-
tions pα and sα which are in `α(I, V ), see the beginning of Section 4. We construct also
s̃α as explained there, with s̃ 0α(0) = s0, the initial condition.

8. Uniform estimates for discrete time solutions

Lemma 8.1. Let (pα, sα) be a time discrete solution of (Sε) in the sense of Definition 4.2.
Then, there exists a positive constant C (independent of α) such that

‖pα‖L2(I;V ) ≤ C, ∀α > 0, (8.1)

‖sα‖L2(I;V ) ≤ C, ∀α > 0, (8.2)

‖s̃α‖L2(I;V ) ≤ C, ∀α > 0, (8.3)

N∑
j=1

‖sjα(·)− sj
′α(·))‖2L2(Ω) ≤ C. (8.4)

Proof. Let us begin by the equation of pressure. Testing Equation (4.2) with ϕ = pjα, for
j = 1, . . . , N , we obtain∫

Ω

λ(sj
′α)K(∇pjα)∇pjα · ∇pjα dx = (qαj , p

jα), t ∈ Ij ,



16 I. MELZI, Y. ATIK EJDE-2022/15

which implies that

λ∗κ2‖∇pjα‖2L2(Ω) ≤ ‖qαj‖L2(Ω) ‖p
jα‖L2(Ω) ≤ Cp‖qαj‖L2(Ω) ‖∇p

jα‖L2(Ω).

Using Young inequality, we obtain

λ∗κ2‖∇pjα‖2L2(Ω) ≤
1

2β
C2
p‖qαj‖2L2(Ω) +

β

2
‖∇pjα‖2L2(Ω), β > 0,

and choosing β = λ∗κ2, we obtain

λ∗κ2

2
‖∇pjα‖2L2(Ω) ≤

1

2λ∗κ2
C2
p‖qαj‖2L2(Ω),

this shows that

λ∗κ2

2
α

N∑
j=1

‖∇pjα‖2L2(Ω) ≤
C2
p

2λ∗κ2
α

N∑
j=1

‖qαj‖2L2(Ω),

it results that

‖pα‖2L2(I;V ) ≤
( Cp
λ∗κ2

)2

‖q‖2L2(I;L2(Ω)).

Remark 8.2. If q ∈ L∞(I;L2(Ω)), then pα ∈ L∞(I;V ) with

‖pα‖2L∞(I;V ) ≤
( Cp
λ∗κ2

)2

‖q‖2L∞(I;L2(Ω)).

Concerning the equation of saturation, we test Equation (4.4) with ϕ = sjα, for j =
1, . . . , N , and obtain(

φ(sjα − sj
′α), sjα

)
+ α

(
λw(sjα)K(∇pjα)∇pjα,∇sjα

)
+ α

(
Λε(s

jα)p′c(s
jα)K(∇pjα)∇sjα,∇sjα

)
= α(qwαj , s

jα).

For the first term, using the identity a(a− b) = 1
2
[a2 − b2 + (a− b)2], we obtain∫

Ω

φ(sjα(x)− sj
′α(x))sjα(x) dx =

∫
Ω

φ
1

2

[
(sjα(x))2 − (sj

′α(x))2 + (sjα(x)− sj
′α(x))2] dx.

Consequently

N∑
j=1

∫
Ω

φ(sjα(x)− sj
′α(x))sjα(x) dx

=

m∑
j=1

∫
Ω

φ
1

2

[
(sjα(x))2 − (sj

′α(x))2 + (sjα(x)− sj
′α(x))2] dx

=
1

2

∫
Ω

φ(sNα(x))2 dx− 1

2

∫
Ω

φ(s0(x))2 dx+
1

2

N∑
j=1

∫
Ω

φ|(sjα(x)− sj
′α(x))|2 dx.

Concerning the second term, summing j from 1 to N , we have∣∣∣ N∑
j=1

α
(
λw(sjα)K(∇pjα)∇pjα,∇sjα

)∣∣∣ ≤ λ∗(κ1

η
+ κ2

) N∑
j=1

α‖∇pjα‖L2(Ω) ‖∇s
jα‖L2(Ω).

Using Hölder’s Inequality, we obtain

λ∗
(κ1

η
+ κ2

) N∑
j=1

α‖∇pjα‖L2(Ω) ‖∇s
jα‖L2(Ω)

≤ λ∗
(κ1

η
+ κ2

)( N∑
j=1

α‖∇pjα‖2L2(Ω)

)1/2( N∑
j=1

α‖∇sjα‖2L2(Ω)

)1/2



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 17

= λ∗
(κ1

η
+ κ2

)
‖pα‖L2(I;V ) ‖s

α‖L2(I;V ),

then applying Young’s Inequality, for β > 0, we obtain∣∣∣ N∑
j=1

α
(
λw(sjα)K(∇pjα)∇pjα,∇sjα

)∣∣∣ ≤ [λ∗(κ1

η
+ κ2

)]2 1

2β
‖pα‖2L2(V ) +

β

2
‖sα‖2L2(V ).

Summing the third term from j = 1 to N , we obtain

N∑
j=1

α
(

Λε(s
jα)p′c(s

jα)K(∇pjα)∇sjα,∇sjα
)
≥ εκ2p

′
c∗

N∑
j=1

α

∫
Ω

|∇sjα|2 dx

= εκ2p
′
c∗

N∑
j=1

α‖∇sjα‖2L2(Ω)

= εκ2p
′
c∗‖sα‖2L2(I;V ).

Finally, for the second member and using Hölder’s Inequality, we obtain∣∣∣ N∑
j=1

α(qwαj , s
jα)
∣∣∣ ≤ N∑

j=1

α‖qwαj‖L2(Ω) ‖s
jα‖L2(Ω)

≤ Cp
N∑
j=1

α‖qwαj‖L2(Ω) ‖∇s
jα‖L2(Ω)

≤ Cp
( N∑
j=1

α‖qwα(tj)‖2L2(Ω)

)1/2( N∑
j=1

α‖∇sjα‖2L2(Ω)

)1/2

= Cp‖qwα‖L2(I;L2(Ω)) ‖s
α‖L2(I;V ).

Then, using that ‖qwα‖L2(I;L2(Ω)) ≤ ‖qw‖L2(I;L2(Ω)) and Young’s Inequality, for β1 > 0,∣∣∣ N∑
j=1

α(qwjα, s
jα)
∣∣∣ ≤ C2

p

2β1
‖qw‖2L2(I;L2(Ω)) +

β1

2
‖sα‖2L2(I;V ).

Taking into account all the previous estimates, after reorganizing terms, we obtain

1

2

∫
Ω

φ(sNα(x))2 dx+
1

2

N∑
j=1

∫
Ω

φ|(sjα(x)− sj
′α(x))|2 dx+ εκ2p

′
c∗‖sα‖2L2(I;V )

≤
[
λ∗
(κ1

η
+ κ2

)]2 1

2β
‖pα‖2L2(I;V ) +

β

2
‖sα‖2L2(I;V )

+
C2
p

2β1
‖qw‖2L2(I;L2(Ω)) +

β1

2
‖sα‖2L2(I;V ) +

1

2

∫
Ω

φ(s0(x))2 dx .

As a result(
εκ2p

′
c∗ −

β

2
− β1

2

)
‖sα‖2L2(I;V ) +

1

2

N∑
j=1

∫
Ω

φ|(sjα(x)− sj
′α(x))|2 dx

≤
C2
p

2β1
‖qw‖2L2(I;L2(Ω)) +

[
λ∗
(κ1

η
+ κ2

)]2 1

2β
‖pα‖2L2(I;V ) +

1

2
φ∗‖s0‖2L2(Ω) .

Taking β = β1 = εκ2p
′
c∗/2 and using (8.1), we see that the estimates (8.2), (8.3) and (8.4)

are valid. �

Lemma 8.3. There exists a constant C > 0 (independent of α) such that∥∥φ∂s̃α
∂t

∥∥
L2(I;V ?)

≤ C, ∀α > 0. (8.5)



18 I. MELZI, Y. ATIK EJDE-2022/15

Proof. For each j = 1, . . . , N , let Lαj be the linear form (and continuous) on V defined by

Lαj (ψ) =

∫
Ω

φ
sjα − sj

′α

α
ψ dx

= −
∫

Ω

λw(sjα)K(∇pjα)∇pjα · ∇ψ dx

−
∫

Ω

Λε(s
jα)p′c(s

jα)K(∇pjα)∇sjα · ∇ψ dx+

∫
Ω

qwαjψ dx, ∀ψ ∈ V.

It follows that for all ψ ∈ V ,∣∣Lαj (ψ)
∣∣ =

∣∣∣− ∫
Ω

λw(sjα)K(∇pjα)∇pjα · ∇ψ dx

−
∫

Ω

Λε(s
jα)p′c(s

jα)K(∇pjα)∇sjα · ∇ψ dx+

∫
Ω

qwαjψ dx
∣∣∣

≤ λ∗
(κ1

η
+ κ2

)
‖∇pjα‖L2‖∇ψ‖L2

+ (ε+ c)p
′∗
c

(κ1

η
+ κ2

)
‖∇sjα‖L2‖∇ψ‖L2 + ‖qwαj‖L2‖ψ‖L2 ,

where c = (λ∗)2/λ∗ and j = 1, . . . , N . Using Poincaré’s Inequality, we obtain∣∣Lαj (ψ)
∣∣ ≤ λ∗(κ1

η
+ κ2

)
‖∇pjα‖L2‖∇ψ‖L2

+
(
ε+ c)p

′∗
c

(κ1

η
+ κ2

)
‖∇sjα‖L2‖∇ψ‖L2 + Cp‖qwαj‖L2‖∇ψ‖L2

=
[
λ∗
(κ1

η
+ κ2

)
‖∇pjα‖L2 + (ε+ c)p

′∗
c

(κ1

η
+ κ2

)
‖∇sjα‖L2

+ Cp‖qwαj‖L2

]
‖∇ψ‖L2 , for j = 1, . . . , N

In what follows, C denotes a constant which can change from one line to another. The
previous inequality can be rewritten as∣∣Lαj (ψ)

∣∣
‖∇ψ‖L2

≤
(κ1

η
+ κ2

){
λ∗‖∇pjα‖L2 + (ε+ c)p

′∗
c ‖∇sjα‖L2

}
+ Cp‖qwαj‖L2

≤ C
{
‖∇pjα‖L2(Ω) + ‖∇sjα‖L2(Ω) + ‖qwαj‖L2(Ω)

}
, j = 1, . . . , N,

where C is the maximum of constants involved in the preceding inequality. Consequently

‖Lαj ‖V ? ≤ C
{
‖∇pjα‖L2(Ω) + ‖∇sjα‖L2(Ω) + ‖qwαj‖L2(Ω)

}
, j = 1, . . . , N

which leads to

‖Lαj ‖2V ? ≤ 3C2{‖∇pjα‖2L2(Ω) + ‖∇sjα‖2L2(Ω) + ‖qwαj‖2L2(Ω)

}
, j = 1, . . . , N. (8.6)

Before going further, it seems good to notice that for each j = 1, . . . , N , the function

φ sjα−sj
′α

α
is in L2(Ω) and then in V ?, since L2(Ω) ↪→ V ?. We can therefore consider that

[0, T ] 3 t 7−→ φs̃α(·, t) is a path in V ? and we have -in the sense of classical derivatives,

V ? 3 φ(·) ∂s̃
α(·, t)
∂t

= φ(·)s
jα(·)− sj

′α(·)
α

, ∀t ∈]tj′ , tj [, j = 1, . . . , N.

Now let us calculate∫ T

0

∥∥φ∂s̃α(·, t)
∂t

∥∥2

V ?
dt =

N∑
j=0

∫ tj

tj′

∥∥φ∂s̃α(·, t)
∂t

∥∥2

V ?
dt

=

N∑
j=0

∫ tj

tj′

∥∥φsjα − sj′α
α

∥∥2

V ?
dt



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 19

=

N∑
j=0

∫ tj

tj′

‖Lαj ‖2V ? dt.

Using estimate (8.6), we obtain∫ T

0

∥∥φ∂s̃α(·, t)
∂t

∥∥2

V ?
dt

≤ C
N∑
j=0

∫ tj

tj′

{
‖∇pjα‖2L2(Ω) + ‖∇sjα‖2L2(Ω) + ‖qwαj‖2L2(Ω)

}
dt

= C
{
‖∇pα‖2L2(ΩT ) + ‖∇sα‖2L2(ΩT ) + ‖qwα‖2L2(ΩT )

}
≤ C.

(8.7)

Here C does not depend on α, because we proved earlier that the sequences {pα}α>0,
{sα}α>0 are bounded in L2(0, T ;V ) and the sequence {qwα}α>0 is bounded in L2(ΩT ). �

9. Compactness of discrete time solutions

First, we give the following remark, see for instance [3].

Remark 9.1. Let w be a function belonging to L2(Ω) and wα the average function in
time defined by relation (4.1). Then

lim
α→0

wα = w in L2(ΩT ) strongly. (9.1)

Lemma 9.2. Let sα satisfy the saturation equation (4.3). Then, there exists a constant
C such that

1√
ξ

∫ T

ξ

∫
Ω

φ{sα(x, t)− sα(x, t− ξ)}2 dx dt ≤ C, ∀ ξ > 0. (9.2)

Proof. We follow [1] (see also [2] and [8]). Let k be fixed (1 ≤ k ≤ N) and let τ ∈]kα, T ],
so there exists j ≥ k + 1 such that τ ∈ Ij = ]tj−1, tj ]. Let R(τ) = ](j − k)α, jα] and take

ω(x, t) = kαχR(τ)(t)∂
−kα
t sα(x, τ) as a test function in the equation of saturation (4.3).

For the parabolic term, we obtain∫
I

(φ∂−αt sα, ω)Ω dt =

∫
I

(
φ∂−αt sα, kαχR(τ)(t)∂

−kα
t sα(x, τ)

)
Ω
dt

=

∫
I

∫
Ω

φ∂−αt sαkαχR(τ)(t)∂
−kα
t sα(x, τ) dx dt

=

∫
Ω

[
φ(x)kα∂−kαt sα(x, τ)

∫
I

∂−αt sα(x, t)χR(τ)(t) dt
]
dx,

with ∫
I

∂−αt sα(x, t)χR(τ)(t) dt =

∫ jα

(j−k)α

∂−αt sα(x, t) dt

=

∫ jα

(j−k)α

sα(x, t)− sα(x, t− α)

α
dt

=

k∑
r=1

∫ (j−k+r)α

(j−k+r−1)α

sα(x, t)− sα(x, t− α)

α
dt

=

k∑
r=1

sα(x, (j − k + r)α)− sα(x, (j − k + r − 1)α)

= sα(x, jα)− sα(x, (j − k)α)

= sα(x, τ)− sα(x, τ − kα)

= kα∂−kαt sα(x, τ) = kα∂−kαt sα(x, jα),



20 I. MELZI, Y. ATIK EJDE-2022/15

which means that the parabolic term is equal to∫
Ω

φ(x)kα∂−kαt sα(x, τ)kα∂−kαt sα(x, jα) dx =

∫
Ω

φ(x)(kα)2
(
∂−kαt sα(x, τ)

)2

dx

By integrating this equality with respect to τ from kα to T , we obtain∫ T

kα

∫
I

(φ∂−αt sα, ω)Ω dt dτ =

∫ T

kα

∫
Ω

φ(x)(kα)2(∂−kαt sα(x, τ)
)2
dx dτ. (9.3)

We have also∫
ΩT

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)
∇ω(x, t) dx dt

=

∫
ΩT

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)
× kαχR(τ)(t)∇∂−kαt sα(x, τ) dx dt

=

∫
Ω

kα∇∂−kαt sα(x, τ)

∫ T

0

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)

×K(∇pα)∇sα
)
χR(τ)(t) dt dx

=

∫
Ω

kα∇∂−kαt sα(x, τ)

∫ jα

(j−k)α

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)
dt dx.

Let

F (x) =

∫ jα

(j−k)α

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)
dt,

G(x) = kα∇∂−kαt sα(x, τ).

Applying Hölder’s inequality, we obtain∫
ΩT

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)
∇ω(x, t) dx dt

≤
(∫

Ω

F 2(x) dx
)1/2 (∫

Ω

G2(x) dx
)1/2

.

(9.4)

Then we have ∫
Ω

G2(x) dx =

∫
Ω

(
kα

sα(x, τ)− sα(x, τ − kα)

kα

)2

dx

=

∫
Ω

(
∇sα(x, τ)−∇sα(x, τ − kα)

)2

dx.

According to the inequality: (a+ b)2 ≤ 2a2 + 2b2, for all a, b real, we have∫
Ω

G2(x) dx ≤
∫

Ω

(
|∇sα(x, τ)|+ |∇sα(x, τ − kα)|

)2

dx

≤
∫

Ω

2|∇sα(x, τ)|2 + 2|∇sα(x, τ − kα)|2 dx.

Using Hölder’s inequality, we obtain∫
Ω

F 2(x) dx

=

∫
Ω

[ ∫ jα

(j−k)α

1×
(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)
dt
]2
dx

≤
∫

Ω

∫ jα

(j−k)α

12 dt×
∫ jα

(j−k)α

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)2
dt dx



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 21

= kα

∫
Ω

∫ jα

(j−k)α

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)2
dt dx.

Taking into account all the previous estimates, (9.4) implies∫
ΩT

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)
∇ω(x, t) dx dt

≤
√

2kα
(∫

Ω

∫ jα

(j−k)α

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)2

dt dx
)1/2

×
(∫

Ω

|∇sα(x, τ)|2 + |∇sα(x, τ − kα)|2 dx
)1/2

≤
√

2kα
(∫

Ω

∫ jα

(j−k)α

(
λ∗(

κ1

η
+ κ2)∇pα + Cεp

′
c
∗
(
κ1

η
+ κ2)∇sα

)2

dt dx
)1/2

×
(∫

Ω

|∇sα(x, τ)|2 + |∇sα(x, τ − kα)|2 dx
)1/2

.

In what follows, C denotes a constant which can change from one line to another.∫
ΩT

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)
∇ω(x, t) dx dt

≤
√

2kαC
(∫

Ω

∫ jα

(j−k)α

(
∇pα +∇sα

)2
dt dx

)1/2

×
(∫

Ω

|∇sα(x, τ)|2 + |∇sα(x, τ − kα)|2 dx
)1/2

,

where C = max
(
λ∗(κ1

η
+ κ2), Cεp

′
c
∗
(κ1
η

+ κ2), λ
∗2

λ∗
+ ε
)

. Then∫
ΩT

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)
∇ω(x, t) dx dt

≤
√

2kαC
(∫

Ω

∫ jα

(j−k)α

|∇pα|2 + |∇sα|2 dt dx
)1/2

×
(∫

Ω

|∇sα(x, τ)|2 + |∇sα(x, τ − kα)|2 dx
)1/2

.

Now, using the fact that, for 0 < p < 1 and a, b two real positive, (a + b)p ≤ ap + bp, we
obtain ∫

ΩT

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)
∇ω(x, t) dx dt

≤
√

2kαC
[( ∫

ΩT

|∇pα|2 dx dt
)1/2

+
(∫

ΩT

|∇sα|2 dx dt
)1/2]

×
[( ∫

Ω

|∇sα(x, τ)|2 dx
)1/2

+
(∫

Ω

|∇sα(x, τ − kα)|2 dx
)1/2]

.

Using (8.1) and (8.2), we have∫
ΩT

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)
∇ω(x, t) dx dt

≤ 2
√

2kαC
[( ∫

Ω

|∇sα(x, τ)|2 dx
)1/2

+
(∫

Ω

|∇sα(x, τ − kα)|2 dx
)1/2]

.

Integrating this inequality with respect to τ from kα to T , and using estimate (8.2), we
obtain∫ T

kα

∫
ΩT

(
λw(sα)K(∇pα)∇pα + Λε(s

α)p′c(s
α)K(∇pα)∇sα

)
∇ω(x, t) dx dt dτ



22 I. MELZI, Y. ATIK EJDE-2022/15

≤
√

2kαC
[ ∫ T

kα

(∫
Ω

|∇sα(x, τ)|2 dx
)1/2

dτ +

∫ T

kα

(∫
Ω

|∇sα(x, τ − kα)|2 dx
)1/2

d τ
]
,

≤
√

2kαC
[ ∫ T

0

(∫
Ω

|∇sα(x, τ)|2 dx
)1/2

dτ +

∫ T−kα

0

(∫
Ω

|∇sα(x, s)|2 dx
)1/2

d s
]
,

≤
√

2kαC
[ ∫ T

0

(∫
Ω

|∇sα(x, τ)|2 dx
)1/2

dτ +

∫ T

0

(∫
Ω

|∇sα(x, s)|2 dx
)1/2

d s
]
,

≤
√

2kαC.
(9.5)

Finally, concerning the second term in the equation of saturation,∫
ΩT

qwαkαχR(τ)(t)∂
−kα
t sα(x, τ) dx dt =

∫
Ω

∫ jα

(j−k)α

qwαkα∂
−kα
t sα(x, τ) dx dt

=

∫
Ω

kα∂−kαt sα(x, τ)
(∫ jα

(j−k)α

qwα dt
)
dx.

Putting

E(x) =

∫ jα

(j−k)α

qwα dt, G(x) = kα∂−kαt sα(x, τ)

and using Hölder’s inequality, we obtain∫
ΩT

qwαkαχR(τ)(t)∂
−kα
t sα(x, τ) dx dt ≤

(∫
Ω

E2(x) dx
)1/2 (∫

Ω

G2(x) dx
)1/2

.

We have (∫
Ω

G2(x) dx
)1/2

≤ Cp
(∫

Ω

∇G2(x) dx
)1/2

.

Using the same techniques as above, we obtain(∫
Ω

G2(x) dx
)1/2

≤ Cp
(∫

Ω

∇G2(x) dx
)1/2

≤ Cp
√

2
[( ∫

Ω

|∇sα(x, τ)|2 dx
)1/2

+
(∫

Ω

|∇sα(x, τ − kα)|2 dx
)1/2]

.

Also, ∫
Ω

E2(x) dx =

∫
Ω

(∫ jα

(j−k)α

qwα(x, t) dt
)2

dx

≤
∫

Ω

(∫ jα

(j−k)α

12 dt×
∫ jα

(j−k)α

q2
wα(x, t) dt

)
dx

=

∫
Ω

kα

∫ jα

(j−k)α

q2
wα(x, t) dt dx

≤
∫

ΩT

kαq2
wα(x, t) dt dx ≤ kα

∫
ΩT

q2
w(x, t) dt dx.

Consequently, for the second member we have∫
ΩT

qwαkαχR(τ)(t)∂
−kα
t sα(x, τ) dx dt

≤ Cp
√

2kα
(∫

ΩT

q2
w(x, t) dt dx

)1/2

×
[( ∫

Ω

|∇sα(x, τ)|2 dx
)1/2

+
(∫

Ω

|∇sα(x, τ − kα)|2 dx
)1/2]

=
√

2kαC
[( ∫

Ω

|∇sα(x, τ)|2 dx
)1/2

+
(∫

Ω

|∇sα(x, τ − kα)|2 dx
)1/2]



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 23

with C = Cp
( ∫

ΩT
q2
w(x, t) dt dx

)1/2
. Integrating the previous inequality with respect to τ

from kα to T , we obtain∫ T

kα

∫
ΩT

qwαkαχR(τ)(t)∂
−kα
t sα(x, τ) dx dt dτ

=
√

2kαC
[ ∫ T

kα

(∫
Ω

|∇sα(x, τ)|2 dx
)1/2

dτ +

∫ T

kα

(∫
Ω

|∇sα(x, τ − kα)|2 dx
)1/2

dτ
]
.

Using the same arguments as in (9.5), we have∫ T

kα

∫
ΩT

qwαkαχR(τ)(t)∂
−kα
t sα(x, τ) dx dt dτ ≤

√
2kαC. (9.6)

Now, taking into account (9.3), (9.5) and (9.6), one obtains∫ T

kα

∫
Ω

φ(x)(kα)2(∂−kαt sα(x, τ)
)2
dx dτ ≤

√
kαC

which implies that

1√
kα

∫ T

kα

∫
Ω

φ(x)(kα)2(∂−kαt sα(x, τ)
)2
dx dτ ≤ C, (9.7)

consequently, according to [1] (see also [2] and [8]), this concludes the proof of (9.2).
In our opinion, the proof can be completed as follows: For a fixed ξ > 0, there exists

k ≥ 0 such that ξ ∈]kα, (k+1)α]. If k ≥ 1 and since we are integrating a positive function,
we have

1√
ξ

∫ T

ξ

∫
Ω

φ{sα(x, t)− sα(x, t− ξ)}2 dx dt

≤ 1√
kα

∫ T

kα

∫
Ω

φ{sα(x, t)− sα(x, t− ξ)}2 dx dt

≤ 2√
kα

∫ T

kα

∫
Ω

φ{sα(x, t)− sα(x, t− kα)}2 dx dt

+
2√
kα

∫ T

kα

∫
Ω

φ{sα(x, t− kα)− sα(x, t− ξ)}2 dx dt.

(9.8)

Since, by construction, the value of sα on each subinterval Ij , j = 1, . . . , N , is equal to

its value at the end of Ij , we have s(k+1)α = sα(·, (k + 1)α) on ]kα, (k + 1)α], and, for

t ∈ ]lα, (l+1)α], sα(t−kα) = s(l−k+1)α and sα(t− ξ) = s(l−k+1)α or s(l−k)α (because t− ξ
belongs to ](l−k−1)α, (l−k+1)α]). Then sα(t−kα)−sα(t−ξ) = 0 or sα(t−kα)−sα(t−ξ) =

s(l−k+1)α − s(l−k)α. Necessarily we have

2√
kα

∫ T

kα

∫
Ω

φ{sα(x, t− kα)− sα(x, t− ξ)}2 dx dt

=
2√
kα

N−k−1∑
m=0

∫ (k+m+1)α

(k+m)α

∫
Ω

φ{sα(x, t− kα)− sα(x, t− ξ)}2 dx dt

≤ 2α‖φ‖∞√
kα

N−k−1∑
m=0

‖s(m+1)α − smα‖2L2(Ω).

Since k ≥ 1, α ≤ 1 and
∑
l ‖s

(l+1)α − slα‖2L2 ≤ C (according to (8.4)), we obtain

2√
kα

∫ T

kα

∫
Ω

φ{sα(x, t− kα)− sα(x, t− ξ)}2 dx dt ≤ 2C‖φ‖∞. (9.9)

Consequently, using (9.7), (9.8) and (9.9), we obtain

1√
ξ

∫ T

ξ

∫
Ω

φ{sα(x, t)− sα(x, t− ξ)}2 dx dt ≤ C. (9.10)



24 I. MELZI, Y. ATIK EJDE-2022/15

If k = 0 then ξ ≤ α with α ≤ 1. We can write

1√
ξ

∫ T

ξ

∫
Ω

φ{sα(x, t)− sα(x, t− ξ)}2 dx dt

=
1√
ξ

∫ α

ξ

∫
Ω

φ{sα(x, t)− sα(x, t− ξ)}2 dx dt

+
1√
ξ

N−1∑
k=1

∫ (k+1)α

kα

∫
Ω

φ{sα(x, t)− sα(x, t− ξ)}2 dx dt.

For t ∈]ξ, α], we have t−ξ ∈]0, α], and then sα(x, t) = sα(x, t−ξ) = s1α. Thus, the integral
on [ξ, α] is zero. For t ∈]kα, (k+1)α], we distinguish two cases: The first is t ∈]kα, kα+ξ],

then t− ξ ∈](k− 1)α, kα], so sα(x, t) = s(k+1)α and sα(x, t− ξ) = skα. The second case is

t ∈]kα+ξ, (k+1)α], then t−ξ ∈]kα, (k+1)α] and therefore sα(x, t) = sα(x, t−ξ) = s(k+1)α.
So, we obtain (remember the estimate (8.4)):

1√
ξ

∫ T

ξ

∫
Ω

φ{sα(x, t)− sα(x, t− ξ)}2 dx dt

=
1√
ξ

N−1∑
k=1

∫ (k+1)α

kα

∫
Ω

φ{sα(x, t)− sα(x, t− ξ)}2 dx dt

=
1√
ξ

N−1∑
k=1

∫ kα+ξ

kα

∫
Ω

φ{sα(x, t)− sα(x, t− ξ)}2 dx dt

≤ ξ‖φ‖∞√
ξ

N−1∑
k=1

‖s(k+1)α − skα‖2L2(Ω) ≤ C.

This completes the proof of estimate (9.2). �

9.1. Passing to the limit in time discretization.

Lemma 9.3. The sequence (pα)α>0 contains a subsequence converging weakly in L2(I;V )
to a function p as α goes to zero.

Proof. According to estimate (8.1), the sequence (pα)α>0 is bounded in L2(I;V ), therefore
it contains a subsequence (denoted in the same way) such that

pα ⇀ p ∈ L2(I;V ) (weakly). (9.11)

�

Lemma 9.4. The sequence (sα)α>0 contains a subsequence converging strongly in L2(ΩT )
to a function s and a.e. in ΩT as α goes to zero.

Proof. Consider the set F = {sα : α > 0} and let the spaces X = V , B = L2(Ω),
Y = X? = V ? = H−1(Ω) (the dual space). It is well known that X ↪→

cont.
B ↪→

comp.
Y . We

have the following

(1) F is uniformly bounded in L2(I;X), i.e., ‖sα‖L2(I;X) ≤ C, with C independent of α).

(2) limξ→0 ‖τξf − f‖L2(0,T−ξ;X?) = 0 uniformly for f ∈ F . Here (τξf)(t) = f(t + ξ). In
fact, using (9.2), we have

φ∗

∫ T

ξ

∫
Ω

(
sα(x, τ)− sα(x, τ − ξ)

)2
dx dτ ≤ C

√
ξ.

Putting σ = τ − ξ, we obtain∫ T−ξ

0

∫
Ω

(
sα(x, σ + ξ)− sα(x, σ)

)2
dx dσ ≤ C

φ∗

√
ξ,



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 25

meaning that ∫ T−ξ

0

‖τξsα − sα‖2L2(Ω) dσ ≤
C

φ∗

√
ξ, ∀α > 0.

Now, since L2(Ω) is continuously embedded in X?, we have ‖ · ‖X? ≤ c‖ · ‖L2(Ω), showing
that ∫ T−ξ

0

‖τξsα − sα‖2X? dσ ≤
C

φ∗

√
ξ, ∀α > 0,

which implies that

lim
ξ→0
‖τξsα − sα‖L2(0,T−ξ;X?) = 0.

Consequently, using [17, Theorem 5, p. 84], we see that F is relatively compact in
L2(I, L2(Ω)) = L2(ΩT ). Therefore, from (sα)α>0, we can extract a subsequence (denoted
in the same way) converging strongly in L2(ΩT ) and a.e. in ΩT to a function s, s ∈
L2(ΩT ). �

9.2. Consequences of estimates and the initial condition. Since the sequences
(sα)α and (s̃α)α are bounded in L2(I;V ), we have

sα ⇀ s weakly in L2(I;V ) and strongly in L2(ΩT ),

s̃α ⇀ s1 in L2(I;V ) and weakly in L2(ΩT ).

Using estimate (8.4) and Remark 4.1, we see that sα − s̃α → 0 in L2(ΩT ), so that s = s1.

On the one hand, from the estimate (8.5) we have φ ∂s̃
α

∂t
⇀ w in L2(I;V ?). On the other

hand, since s̃α ⇀ s in L2(I;V ), one can deduce that s̃α ⇀ s in L2(I;L2(Ω)) and s̃α ⇀ s in
L2(I;V ?); consequently, s̃α ⇀ s in D′(I;V ?) (the space of distributions on I with values
in V ?) and ∂t(φs̃

α) ⇀ ∂t(φs) in D′(I;V ?). Then w = ∂t(φs) = φ∂ts.
Now, since s̃α ⇀ s in L2(I;V ) and ∂ts̃

α ⇀ ∂ts in L2(I;V ?), we see that s̃α ⇀ s in
W (0, T ). Let us recall that W (0, T ) ↪→cont C([0, T ];L2(Ω)) (see, for instance, [15, 9, 14]),
and, for ξ a fixed element in D(Ω) , consider the linear functional Fξ defined by W (0, T ) 3
u 7→ Fξ(u) =

∫
Ω
u(0)(x)ξ(x) dx ∈ R. If we write∣∣∣ ∫

Ω

u(0)(x)ξ(x) dx
∣∣∣ ≤ ‖u(0)‖L2(Ω)‖ξ‖L2(Ω)

≤ sup
0≤t≤T

‖u(t)‖L2(Ω)‖ξ‖L2(Ω)

≤ C‖ξ‖L2(Ω)‖u‖W (0,T ),

where C is a positive constant, We see that Fξ is continuous with ‖Fξ‖(W (0,T ))? ≤
C‖ξ‖L2(Ω). Therefore,

lim
α↓0

Fξ(s̃
α) = lim

α↓0

∫
Ω

s̃α(0)(x)ξ(x) dx =

∫
Ω

s0(x)ξ(x) dx = Fξ(s) =

∫
Ω

s(0)(x)ξ(x) dx.

As ξ is arbitrary in D(Ω), we conclude that s(0) = s0, see for instance [6, Corollary 4.24].
The initial condition is thus satisfied.

10. Proof of Theorem 3.2

1. Equation of pressure. First, let us remember the approximate equations of pressure∫
Ω

λ(sj
′α)K(∇pjα)∇pjα · ∇ϕdx =

∫
Ω

qαj ϕdx, ∀ϕ ∈ V, ∀j = 1, . . . , N.

Let ψ in D(I;V ), then for all t ∈ [tj−1, tj [, ψ(t) ∈ V , by taking it as a test function in the
equation above, we obtain∫

Ω

λ(sj
′α)K(∇pjα)∇pjα · ∇ψ(t) dx =

∫
Ω

qαj ψ(t) dx, for j = 1, . . . , N.



26 I. MELZI, Y. ATIK EJDE-2022/15

Integrating with respect to t from tj−1 to tj and then by summing j from 1 to N , we have

N∑
j=1

∫ tj

tj−1

∫
Ω

λ(sj
′α)K(∇pjα)∇pjα · ∇ψ(t) dx dt =

N∑
i=1

∫ tj

tj−1

∫
Ω

qαj ψ(t) dx dt,

which is equivalent to∫
ΩT

λ(sα(t− α))K(∇pα)∇pα · ∇ψ dx dt =

∫
ΩT

qα ψ dx dt with ψ ∈ D(I;V ).

Let us make α go to zero. Using hypothesis (H2) on the function λ, Lemma 9.4 and
denoting ζ the weak limit of the sequence K(∇pα)∇pα in L2(ΩT ), we obtain

lim
α→0

∫
ΩT

λ(sα(t− α))K(∇pα)∇pα · ∇ψ dx dt

=

∫
ΩT

λ(s)ζ · ∇ψ dx dt, ∀ψ ∈ D(I;V ).

(10.1)

By Remark 9.1, we have

lim
α→0

∫
ΩT

qα ψ dx dt =

∫
ΩT

q ψ dx dt, ∀ψ ∈ D(I;V ). (10.2)

Combining (10.1)and(10.2), we see that∫
ΩT

λ(s)ζ · ∇ψ dx dt =

∫
ΩT

q ψ dx dt, ∀ψ ∈ D(I;V ), (10.3)

which implies ∫
ΩT

λ(s)ζ · ∇ψ dx dt =

∫
ΩT

q ψ dx dt, ∀ψ ∈ L2(I;V ). (10.4)

Now, taking pjα as a test function in the equation of pressure, then integrating with
respect to t from tj−1 to tj and then by summing j from 1 to N , we have∫

ΩT

λ(sα(t− α))K(∇pα)|∇pα|2 dx dt =

∫
ΩT

qαp
α dx dt, (10.5)

from (8.1), remark 9.1, lemma 9.3 and passing to the limit when α goes to zero, we have

lim
α→0

∫
ΩT

λ(sα(t− α))K(∇pα)|∇pα|2 dx dt =

∫
ΩT

qp dx dt. (10.6)

To justify that limα→0

∫
ΩT

qαp
α dx dt =

∫
ΩT

qp dx dt, we write∣∣∣ ∫
ΩT

qαp
α − qp

∣∣∣ =
∣∣∣ ∫

ΩT

qαp
α − qpα + qpα − qp

∣∣∣
=
∣∣∣ ∫

ΩT

pα(qα − q) +

∫
ΩT

q(pα − p)
∣∣∣

≤ ‖pα‖ ‖qα − q‖+
∣∣∣ ∫

ΩT

q(pα − p)
∣∣∣

= C‖qα − q‖+
∣∣∣ ∫

ΩT

q(pα − p)
∣∣∣.



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 27

For all ϕ in L2(I;V ), we have

0 ≤
∫

ΩT

λ(sα(t− α))(K(∇pα)∇pα −K(∇ϕ)∇ϕ) · (∇pα −∇ϕ) dx dt

=

∫
ΩT

λ(sα(t− α))K(∇pα)|∇pα|2 dx dt

−
∫

ΩT

λ(sα(t− α))K(∇pα)∇pα · ∇ϕdx dt

−
∫

ΩT

λ(sα(t− α))K(∇ϕ)∇ϕ · (∇pα −∇ϕ) dx dt,

(10.7)

using (10.5), the above inequality becames

0 ≤
∫

ΩT

qαp
α dx dt−

∫
ΩT

λ(sα(t− α))K(∇pα)∇pα · ∇ϕdx dt

−
∫

ΩT

λ(sα(t− α))K(∇ϕ) · ∇ϕ(∇pα −∇ϕ) dx dt,

passing to the limit as α goes to zero, we obtain

0 ≤
∫

ΩT

qp dx dt−
∫

ΩT

λ(s)ζ · ∇ϕdx dt−
∫

ΩT

λ(s)K(∇ϕ)∇ϕ · (∇p−∇ϕ) dx dt.

Using (10.4), the previous inequality is equivalent to

0 ≤
∫

ΩT

λ(s)ζ · ∇p dx dt−
∫

ΩT

λ(s)ζ∇ϕdx dt−
∫

ΩT

λ(s)K(∇ϕ)∇ϕ · (∇p−∇ϕ) dx dt.

Consequently,

0 ≤
∫

ΩT

λ(s)(ζ −K(∇ϕ)∇ϕ)(∇p−∇ϕ) dx dt.

Taking ϕ = p− ηϕ̃ with η > 0, we obtain∫
ΩT

λ(s)
(
ζ −K(∇ϕ)(∇p− η∇ϕ̃)

)
· ∇ϕ̃ ≥ 0, (10.8)

when η approaches zero. Using the continuity of the operator A defined in subsection 5.1,
we have ∫

ΩT

λ(s)(ζ −K(∇p)∇p) · ∇ϕ̃ ≥ 0, ∀ϕ̃ ∈ L2(I;V ).

Then, replacing ϕ̃ by −ϕ̃ in (10.8) and making η tend to zero, we deduce the equality∫
ΩT

λ(s)(ζ −K(∇p)∇p) · ∇ϕ̃ = 0,

which is exactly∫
ΩT

λ(s)ζ · ∇ϕ̃ dx dt =

∫
ΩT

λ(s)K(∇p)∇p · ∇ϕ̃ dx dt, ∀ϕ̃ ∈ L2(I;V ).

This shows that ζ = K(∇p)∇p, and as a result, (10.4) becomes∫
ΩT

λ(s)K(∇p)∇p · ∇ψ dx dt =

∫
ΩT

q ψ dx dt, ∀ψ ∈ L2(I;V ); (10.9)

hence, p satisfies the equation of pressure.

Lemma 10.1. The following holds,

lim
α→0

∫
ΩT

λ(sα(t− α))
(
K(∇pα)∇pα −K(∇p)∇p

)(
∇pα −∇p

)
= 0.



28 I. MELZI, Y. ATIK EJDE-2022/15

Proof. From (10.6) and (10.9), we deduce that

lim
α→0

∫
ΩT

λ(sα(t− α))K(∇pα)|∇pα|2 dx dt =

∫
ΩT

qp dx dt

=

∫
ΩT

λ(s)ζ∇p dx dt

=

∫
ΩT

λ(s)K(∇p)|∇p|2 dx dt.

(10.10)

Using (10.4), (10.9), (10.10), and Lemma 9.3, we obtain

0 ≤ lim
α→0

∫
ΩT

λ(sα(t− α))
(
K(∇pα)∇pα −K(∇p)∇p

)
·
(
∇pα −∇p

)
dx dt

= lim
α→0

∫
ΩT

λ(sα(t− α))K(∇pα)|∇pα|2 dx dt

− lim
α→0

∫
ΩT

λ(sα(t− α))K(∇pα)∇pα · ∇p dx dt

− lim
α→0

∫
ΩT

λ(sα(t− α))K(∇p)∇p ·
(
∇pα −∇p

)
=

∫
ΩT

λ(s)K(∇p)|∇p|2 dx dt−
∫

ΩT

λ(s)K(∇p)∇p · ∇p dx dt

−
∫

ΩT

λ(s)K(∇p)∇p
(
∇p−∇p

)
= 0. �

Lemma 10.2. The sequence (∇pα)α>0 converges in measure to ∇p in ΩT and a.e. for a
subsequence.

Proof. Let ε1, δ > 0 be two fixed numbers, and set

D = {|∇pα −∇p| ≥ δ} .= {(x, t) ∈ ΩT : |∇pα(x, t)−∇p(x, t)| ≥ δ}.

Then we consider the function K12 introduced in subsection 5.1: K12(x) = κ1
|x|x

1+η|x|+κ2x,

x ∈ R3. The same method used to prove the monotony of operator A in the mentioned
subsection, shows that

(K12(x)−K12(y)) · (x− y) ≥ κ2|x− y|2, ∀x, y ∈ R3.

Writing ∫
ΩT

λ(sα(t− α))[K12(∇pα)−K12(∇p)] · [∇pα −∇p] dx dt

≥ λ∗
∫
D

[K12(∇pα)−K12(∇p)] · [∇pα −∇p] dx dt

≥ λ∗
∫
D

κ2|∇pα −∇p|2 dx dt ≥ λ∗κ2δ
2 meas(D),

we see that

meas(D) ≤ 1

λ∗κ2δ2

∫
ΩT

λ(sα(t− α))[K12(∇pα)−K12(∇p)] · [∇pα −∇p] dx dt.

Since the right-hand side of the previous inequality tends to zero as α does, this by
Lemma 10.1, meas(D) can be made less than ε1 for α sufficiently small. We conclude
that meas({|∇pα −∇p| ≥ δ}) ≤ ε1, for all ε1 > 0. This proves that the sequence (∇pα)
converges in measure to ∇p. Therefore this sequence contains a subsequence, denoted in
the same way, converging a.e. to ∇p in ΩT . �



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 29

Remark 10.3. Using Lemma 10.2, one can have that K(∇pα) converge a.e. in ΩT to
K(∇p) and since K(∇pα) ≤ κ1

η
+ κ2, we deduce K(∇pα) converge strongly to K(∇p) in

L2(ΩT ).

2. Equation of saturation. First we recall a well known result.

Lemma 10.4 (A discrete integration by parts formula). For α > 0, T > 0 two real
numbers and Φ a smooth real function defined on the interval [0, T ], let us put

∂αt Φ(t) =
Φ(t+ α)− Φ(t)

α
, t ∈ [0, T − α].

If α < T and Ψ is a real smooth function defined on the same interval [0, T ], then∫ T

α

Φ(t)∂−αt Ψ(t) dt

=
1

α

∫ T

T−α
(ΦΨ)(t) dt− 1

α

∫ α

0

(ΦΨ)(t) dt−
∫ T−α

0

Ψ(t)∂αt Φ(t) dt.

(10.11)

Proof. For 0 < t < T − α, we can write

(ΦΨ)(t+ α)− (ΦΨ)(t) = αΦ(t+ α)∂αt Ψ(t) + αΨ(t)∂αt Φ(t).

Integrating the left-hand side on [0, T − α], we obtain∫ T−α

0

(ΦΨ)(t+ α) dt−
∫ T−α

0

(ΦΨ)(t) dt =

∫ T

T−α
(ΦΨ)(s) ds−

∫ α

0

(ΦΨ)(t) dt.

Now, integrating on the right-hand side, we have

α

∫ T−α

0

Φ(t+ α)∂αt Ψ(t) dt+ α

∫ T−α

0

Ψ(t)∂αt Φ(t)

= α

∫ T

α

Φ(t)∂−αt Ψ(t) dt+ α

∫ T−α

0

Ψ(t)∂αt Φ(t).

Putting the results together, we obtain formula (10.11). �

Now, let us remember the approximate equation of saturation∫ T

0

(φ∂−αt sα, ψ) dt+

∫
ΩT

λw(sα)K(∇pα)∇pα · ∇ψ dx dt

+

∫
ΩT

Λε(s
α)p′c(s

α)K(∇pα)∇sα · ∇ψ dx dt

=

∫ T

0

(qwα, ψ) dt, ∀ψ ∈ `α(I;V ).

Following [3], the pressure equation in the weak sense of Definition 3.1 can be seen to
hold since ∪∞n=1`

α(I;V ) (remember that α = T
N

= T
2n

) is dense in L2(I;V ). Also, for the
equation of saturation, and for all ψ ∈ ∪∞n=1`

α(I;V ), we have for the second term, using
the same technique when passing to the limit in the equation of pressure, we obtain

lim
α→0

∫
ΩT

λw(sα)K(∇pα)∇pα · ∇ψ dx dt =

∫
ΩT

λw(s)K(∇p)∇p · ∇ψ dx dt,

for the third term, using Remark 10.3, we obtain

lim
α→0

∫
ΩT

Λε(s
α)p′c(s

α)K(∇pα)∇sα · ∇ψ dx dt =

∫
ΩT

Λε(s)p
′
c(s)K(∇p)∇s · ∇ψ dx dt.

For the last term, using Remark 9.1, we have

lim
α→0

∫ T

0

(qwα, ψ) dt =

∫ T

0

(qw, ψ) dt.



30 I. MELZI, Y. ATIK EJDE-2022/15

It follows from (4.3) that

lim
α→0

∫ T

0

(φ∂−αt sα, ψ) dt+

∫
ΩT

λw(s)K(∇p)∇p · ∇ψ dx dt

+

∫
ΩT

Λε(s)p
′
c(s)K(∇p)∇s · ∇ψ dx dt

=

∫ T

0

(qw, ψ), ∀ψ ∈ ∪∞n=1`
α(I;V ).

For any ψ ∈ L2(I;V ), ψα ∈ `α(I;V ), and because sα(·, t) is constant over each interval
Ij = (tj−1, tj ], we observe that∫ T

0

(φ∂−αt sα, ψ) dt =

∫ T

0

(φ∂−αt sα, ψα) dt, (10.12)

then, identity (4.3) can be written as∫ T

0

(φ∂−αt sα, ψα) dt = −
∫

ΩT

λw(s)K(∇p)∇p · ∇ψα dx dt

−
∫

ΩT

Λε(s)p
′
c(s)K(∇p)∇s · ∇ψα dx dt+

∫ T

0

(qw, ψα).

This implies that ∣∣ ∫ T

0

(φ∂−αt sα, ψ) dt
∣∣ ≤ C‖ψ‖L2(I;V ), ∀ψ ∈ L2(I;V ).

The sequence (φ∂−αt sα) is thus bounded in L2(I;V ?). Consequently, for a subsequence,
(φ∂−αt sα) converges weakly in L2(I;V ?). For ψ ∈ D(I;V ) and α > 0 small enough, using
Formula (10.11), we have∫ T

0

(φ∂−αt sα(·, t), ψ(·, t)) dt = −
∫ T−α

0

(φsα, ∂αt ψ) dt

→ −
∫ T

0

(φs, ∂tψ) dt =

∫ T

0

〈φ∂ts, ψ〉 dt,

as a distribution. Therefore, φ∂−αt sα ⇀ φ∂ts weakly in L2(I;V ?). Combining these
results, the saturation equation holds in the weak sense of Definition 3.1 since ∪∞n=1`

α(I;V )
is dense in L2(I;V ). Thus the proof of Theorem 3.2 is complete.

11. Maximum principles about weak solutions

Theorem 11.1. If (p, s) is a weak solution of system (Sε), then 0 ≤ s(x, t) ≤ 1 a.e. x in
Ω and for all t in [0, T ].

Proof. To show that s(x, t) ≥ 0, we prove that its negative part s− is zero on ΩT . Let us
first remark that for (p, s) a weak solution of system (Sε), the equation of saturation (3.3)
implies that

〈φ∂ts, v〉 −
∫

Ω

λw(s)K(∇p)∇p · ∇v dx−
∫

Ω

Λε(s)p
′
c(s)K(∇p)∇s · ∇v dx

=

∫
Ω

qwv dx, ∀v ∈ V, a.e. in (0, T ).

(11.1)

Let us fix a number t ∈ (0, T ]. Since the function R 3 r 7−→ 1
2
(|r| − r) .

= r− ∈ R is

Lipschitz, it is licit to take v = −s−(σ), with σ ∈ (0, t) non exceptional, as a test function
in the Equation (11.1) written at the time point σ. We obtain

〈φ∂ts(σ),−s−(σ)〉 −
∫

Ω

λw(s)K(∇p)∇p · ∇s−dx−
∫

Ω

Λε(s)p
′
c(s)K(∇p)∇s · ∇s− dx



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 31

= −
∫

Ω

qws
−dx.

Since λw(s) = 0 for s ≤ 0, we obtain
∫

Ω
λw?(s)K(∇p)∇p ·∇s− dx dt = 0. Also, using that

−
∫

Ω

Λε(s)p
′
c(s)K(∇p)∇s · ∇s− dx =

∫
Ω

Λε(s)p
′
c(s)K(∇p)|∇s−|2 dx

≥ εp′c∗κ2

∫
Ω

|∇s−|2 dx,

and the positivity of function qw (hypothesis (H4)), we obtain

〈φ∂ts(σ),−s−(σ)〉+ εp′c∗κ2

∫
Ω

|∇s−|2 dx ≤ −
∫

Ω

qws
−dx ≤ 0.

We deduce that
〈φ∂ts(σ),−s−(σ)〉 ≤ 0 a.e. σ ∈ (0, t). (11.2)

Note that

〈φ∂ts(σ),−s−(σ)〉 =
1

2

∂

∂t

∫
Ω

φ(x)|s−(x, σ)|2 dx a.e. σ ∈ (0, T ). (11.3)

Let us now suppose that s ∈ D([0, T ];V ), the space of restrictions to [0, T ] of functions
indefinitely differentiable with compact support and values in V . In this case, we can
write 〈

φ
∂s

∂t
(σ),−s−(σ)

〉
= −

∫
Ω

φ(x)
∂s

∂t
(x, σ)s−(x, σ) dx

= −
∫

Ω∩{s(x,σ)<0}
φ(x)

∂s

∂t
(x, σ)s−(x, σ) dx

=

∫
Ω∩{s(x,σ)<0}

φ(x)
∂s−

∂t
(x, σ)s−(x, σ) dx

=

∫
{s<0}

φ(x)
1

2

∂

∂t
|s−(x, σ)|2 dx

=
1

2

∂

∂t

∫
Ω

φ(x)|s−(x, σ)|2 dx

Adopting the techniques used by Chipot [9, Lemma 11.2, page 203], we can prove that
the (11.3) remains true for s ∈ W (0, T ) = H1(0, T ;V, V ?). Integrating inequality (11.2),
we obtain∫ t

0

〈
φ
∂s

∂t
(σ),−s−(σ)

〉
dσ =

1

2

∫
Ω

φ(x)|s−(x, t)|2 dx− 1

2

∫
Ω

φ(x)|s−(x, 0)|2 dx

=
1

2

∫
Ω

φ(x)|s−(x, t)|2 dx ≤ 0.

This because s(0) = s0(x) ≥ 0 (hypothesis (H4)), giving
∫

Ω
φ(x)|s−(x, 0)|2 dx = 0, and

the Inequality (11.2). Now, using (H1), we see that
∫

Ω
|s−(x, t)|2 dx = 0, a.e. in t ∈]0, T ].

This proves that s(x, t) ≥ 0 a.e. in ΩT .
To show that s(x, t) ≤ 1 a.e., we prove that (s − 1)+, the positive part of s − 1, is

zero on ΩT . Using the same techniques as before, we fix a number t ∈ (0, T ], and take
v = (s − 1)+(σ), with σ ∈ (0, t) non exceptional, as a test function in (11.1) written at
the time point σ. We obtain

〈φ∂ts(σ), (s− 1)+(σ)〉+

∫
Ω

λw(s)K(∇p)∇p · ∇(s− 1)+dx

+

∫
Ω

Λε(s)p
′
c(s)K(∇p)∇s · ∇(s− 1)+dx

=

∫
Ω

qw(s− 1)+dx.

(11.4)



32 I. MELZI, Y. ATIK EJDE-2022/15

Putting {s > 1} = Ω ∩ {s(x, σ) > 1}, by the hypothesis (H2) and the definition of
extensions of the coefficients in equations, see Section 3, the equation of pressure leads to∫

Ω

λw(s)K(∇p)∇p · ∇(s− 1)+dx =

∫
{s>1}

λw(s)K(∇p)∇p · ∇(s− 1)+dx

=

∫
{s>1}

λ(s)K(∇p)∇p · ∇(s− 1)+dx

=

∫
{s>1}

q(s− 1)+dx =

∫
Ω

q(s− 1)+dx.

Now, using the inequality∫
Ω

Λε(s)p
′
c(s)K(∇p)∇s · ∇(s− 1)+dx ≥ εp′c∗κ2

∫
Ω

|∇(s− 1)+|2 dx,

and equation (11.4), we obtain

〈φ∂ts(σ), (s− 1)+(σ)〉+

∫
Ω

q(s− 1)+dx+ εp′c∗κ2

∫
Ω

|∇(s− 1)+|2 dx ≤
∫

Ω

qw(s− 1)+dx

Therefore, since q − qw = qn, which is a positive function, we obtain

〈φ∂ts(σ), (s− 1)+(σ)〉+ εp′c∗κ2

∫
Ω

|∇(s− 1)+|2 dx ≤ −
∫

Ω

qn(s− 1)+dx ≤ 0.

Consequently,

〈φ∂ts(σ), (s− 1)+(σ)〉 ≤ 0 a.e. σ ∈ (0, t). (11.5)

To go further, we note as above that

〈φ∂ts(σ), (s− 1)+(σ)〉 = 〈φ∂t(s− 1)(σ), (s− 1)+(σ)〉

=
1

2

∂

∂t

∫
Ω

φ(x)|(s− 1)+(x, σ)|2 dx a.e. σ ∈ (0, T ).
(11.6)

This can be seen using the denseness for s ∈ D([0, T ];V ) in W (0, T ).
Integrating the previous (11.5), we obtain∫ t

0

〈
φ
∂s

∂t
(σ), (s− 1)+(σ)

〉
dσ

=
1

2

∫
Ω

φ(x)|(s− 1)+(x, t)|2 dx− 1

2

∫
Ω

φ(x)|(s− 1)+(x, 0)|2 dx

=
1

2

∫
Ω

φ(x)|(s− 1)+(x, t)|2 dx ≤ 0.

This because s(0) = s0(x) ≤ 1 (hypothesis (H4)), giving
∫

Ω
φ(x)|(s − 1)+(x, 0)|2 dx = 0,

and the Inequality (11.5). Now, using Hypothesis (H1), we see that
∫

Ω
|(s−1)+(x, t)|2 dx =

0, a.e. in t ∈]0, T ]. This proves that s(x, t) ≤ 1 a.e. in ΩT . �

Remark 11.2. To finish, we mention that all results of this paper are in fact true for
a family of absolute permeability, in the sense that our results remain true if we replace
the expression of absolute rock permeability given in page 2 by any continuous function
K : R3 −→ R, bounded from below and above by positive constants with K(x)x monotone.

Acknowledgements. The authors would like to thank Professor M. Saad, of École Cen-
trale de Nantes, for his help in carrying out this work. In particular for suggesting the
expression of absolute permeability given on page 2.



EJDE-2022/15 TWO-PHASE FLOW IN NANOPOROUS MEDIA 33

References

[1] H. W. Alt, S. Luckhaus; Quasilinear elliptic-parabolic differential equations, Math. Z. 183,
(1983), 311–341.

[2] J. Arbogast; The existence of weak solutions to single porosity and simple dual-porosity

models of two-phase incompressible flow, Journal of Nonlinear Analysis: Theory, Methods,
and Applications, 19 (1992), 1009–1031.

[3] Y. Atik; An existence result for heat and mass transfer in porous media, ENS-Kouba, January

(2016), Unpublished.
[4] G. I. Barenblatt, P. J. M. Monteiro, C. H. Rycroft; On a boundary layer problem related to

the gas flow in shales, J Eng Math, 84 (2013), 11–18.
[5] L. Boccardo, T. Gallouët; Nonlinear elliptic equations with right hand side measures, Com-

munications in Partial Differential Equations, 17 (3-4) (1992), 189–258.

[6] H. Brezis; Functional Analysis, Sobolev Spaces and Partial Differential Equations, Springer,
New York, 2011.

[7] Z. Chen, G. Huan, Y. Ma; Computational Methods for Multiphase Flows in Porous Media,

SIAM, Philadelphia, 2006.
[8] Z. Chen; Degenerate two-phase incompressible flow I, existence, uniqueness and regularity of

a weak solution, Journal of Differential Equations, 171 (2001), 203–232.

[9] M. Chipot; Elements of Nonlinear Analysis, Birkhäuser-Verlag, Basel, 2000.
[10] R. Dautray, J. L. Lions; Mathematical Analysis and Numerical Methods for Science and

Technology, Vol. 5: Evolution Problems I, Springer-Verlag, Berlin, 2000.

[11] H. Kabir; Étude théorique de couplages entre écoulements Forchheimer et déformations

mécaniques dans l’extraction d’hydrocarbures, PhD thesis, ENS B. El Ibrahimi, Kouba, Al-

giers, 2020.
[12] J. L. Lions; Quelques Méthodes de Résolution des Problèmes aux Limites Non Linéaires,

Dunod & Gauthier-Villars, Paris, 1969.

[13] P. J. M. Monteiro, C. H. Rycroft, G. I. Barenblatt; A Mathematical model of fluid and gas
flow in nanoporous media, PNAS Early Edition, (2012), 1–5.

[14] T. Roubicek; Nonlinear Partial Differential Equations with Applications, Vol. 153. Springer,

2013.
[15] S. Salsa; Partial Differential Equations in Action, From Modeling to Theory, Springer, Italy,

2008.
[16] R. E. Showalter; Monotone Operators in Banach Space and Nonlinear Partial Differential

Equations, AMS, Providence, Rhode Island, 1997.

[17] J. Simon; Compact sets in the space Lp(0, T ;B), Annali di Matematica pura ed applicata
(IV), Vol. CXLVI (1987), 65–96.

[18] R. S. Strichartz; The Way of Analysis, Jones and Bartlett Publishers International, London,

2000.

Imane Melzi
Laboratory of nonlinear partial differential equations & history of mathematics, de-
partment of mathematics, Ecole Normale Supérieure, Kouba, Algiers, Algeria

Email address: imane.melzi@g.ens-kouba.dz

Youcef Atik

Laboratory of nonlinear partial differential equations & history of mathematics, de-
partment of mathematics, Ecole Normale Supérieure, Kouba, Algiers, Algeria

Email address: youcefatiq@gmail.com


	1. Introduction
	1.1. Flow equations
	1.2. Functional setting

	2. Hypotheses
	3. Regularization: system (S)
	3.1. Weak solution of system (S)

	4. Time discretization of system (S)
	5. Galerkin's approximations of discrete time solutions
	5.1. Existence of Galerkin's approximations

	6. Uniform estimates on Galerkin's approximations
	7. Passing to the limit (with respect of d) in Galerkin's approximations
	8. Uniform estimates for discrete time solutions
	9. Compactness of discrete time solutions
	9.1. Passing to the limit in time discretization
	9.2. Consequences of estimates and the initial condition

	10. Proof of Theorem ??
	11. Maximum principles about weak solutions
	Acknowledgements

	References