Electronic Journal of Differential Equations, Vol. 2016 (2016), No. 52, pp. 1–28.

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

HOMOGENIZATION OF IMMISCIBLE COMPRESSIBLE TWO-PHASE FLOW IN DOUBLE POROSITY MEDIA

LATIFA AIT MAHIOUT, BRAHIM AMAZIANE, ABDELHAFID MOKRANE, LEONID PANKRATOV

Abstract. A double porosity model of multidimensional immiscible com- pressible two-phase flow in fractured reservoirs is derived by the mathematical theory of homogenization. Special attention is paid to developing a general approach to incorporating compressibility of both phases. The model is writ- ten in terms of the phase formulation, i.e. the saturation of one phase and the pressure of the second phase are primary unknowns. This formulation leads to a coupled system consisting of a doubly nonlinear degenerate para- bolic equation for the pressure and a doubly nonlinear degenerate parabolic diffusion-convection equation for the saturation, subject to appropriate bound- ary and initial conditions. The major difficulties related to this model are in the doubly nonlinear degenerate structure of the equations, as well as in the coupling in the system. Furthermore, a new nonlinearity appears in the tem- poral term of the saturation equation. The aim of this paper is to extend the results of [9] to this more general case. With the help of a new compactness re- sult and uniform a priori bounds for the modulus of continuity with respect to the space and time variables, we provide a rigorous mathematical derivation of the upscaled model by means of the two-scale convergence and the dilatation technique.

1. Introduction

The modeling of displacement process involving two immiscible fluids is of con- siderable importance in groundwater hydrology and reservoir engineering such as petroleum and environmental problems. More recently, modeling multiphase flow received an increasing attention in connection with gas migration in a nuclear waste repository and sequestration of CO2. Furthermore, fractured rock domains corre- sponding to the so-called Excavation Damaged Zone (EDZ) receives increasing at- tention in connection with the behaviour of geological isolation of radioactive waste after the drilling of the wells or shafts, see, e.g., [50].

A fissured medium is a structure consisting of a porous and permeable matrix which is interlaced on a fine scale by a system of highly permeable fissures. The majority of fluid transport will occur along flow paths through the fissure system, and the relative volume and storage capacity of the porous matrix is much larger

2010Mathematics Subject Classification. 35B27, 35K65, 76S05, 76T10.

Key words and phrases. Compressible immiscible; double porous media; two-phase flow;

fractured media homogenization; two-scale convergence.

c

2016 Texas State University.

Submitted February 2, 2016. Published February 18, 2016.

1

than that of the fissure system. When the system of fissures is so well developed that the matrix is broken into individual blocks or cells that are isolated from each other, there is consequently no flow directly from cell to cell, but only an exchange of fluid between each cell and the surrounding fissure system. Therefore the large- scale description will have to incorporate the two different flow mechanisms. For some permeability ratios and fissure widths, the large-scale description is achieved by introducing the so-called double porosity model. It was introduced first for describing the global behaviour of fractured porous media by Barenblatt et al.

[16]. It has been since used in a wide range of engineering specialties related to geohydrology, petroleum reservoir engineering, civil engineering or soil science. For more details on the physical formulation of such problems see, e.g., [17, 49, 51].

During recent decades mathematical analysis and numerical simulation of mul- tiphase flows in porous media have been the subject of investigation of many re- searchers owing to important applications in reservoir simulation. There is an extensive literature on this subject. We will not attempt a literature review here but will merely mention a few references. Here we restrict ourselves to the mathe- matical analysis of such models. We refer, for instance, to the books [13, 27, 31, 36, 43, 45, 52] and the references therein. The mathematical analysis and the homog- enization of the system describing the flow of two incompressible immiscible fluids in porous media is quite understood. Existence, uniqueness of weak solutions to these equations, and their regularity has been shown under various assumptions on physical data; see for instance [3, 13, 14, 25, 27, 28, 29, 36, 48] and the references therein. A recent review of the mathematical homogenization methods developed for incompressible immiscible two-phase flow in porous media and compressible miscible flow in porous media can be viewed in [4, 44, 45]. We refer for instance to [18, 19, 20, 21, 22, 41, 42] for more information on the homogenization of incom- pressible, single phase flow through heterogeneous porous media in the framework of the geological disposal of radioactive waste.

The double porosity problem was first studied in [15], and was then revisited in the mathematical literature by many other authors. Here we restrict ourself to the mathematical homogenization method as described in [45] for flow and transport in porous media. For a recent review of the methods developed for flow through double porosity media, we refer for instance to [12, 15, 23, 30, 32, 34, 53] and the references therein.

However, as reported in [9], the situation is quite different for immiscible com- pressible two-phase flow in porous media, where, only recently few results have been obtained. In the case of immiscible two-phase flows with one (or more) compressible fluids without any exchange between the phases, some approximate models were studied in [37, 38, 39]. Namely, in [37] certain terms related to the compressibility are neglected, and in [38, 39] the mass densities are assumed not to depend on the physical pressure, but on Chavent’s global pressure. In the articles [26, 40, 46, 47], a more general immiscible compressible two-phase flow model in porous media is considered for fields with a single rock type and [10] treated the case with several types of rocks. In [4, 11] homogenization results were obtained for water-gas flow in porous media using the phase formulation, i.e. where the phase pressures and the phase saturations are primary unknowns.

Let us also mention that, recently, a different approach based on a new global pressure concept was introduced in [5, 7] for modeling immiscible, compressible

two-phase flow in porous media without any simplifying assumptions. The resulting equations are written in a fractional flow formulation and lead to a coupled system which consists of a nonlinear parabolic equation (the global pressure equation) and a nonlinear diffusion-convection one (the saturation equation). This new formulation is fully equivalent to the original phase equations formulation, i.e. where the phase pressures and the phase saturations are primary unknowns. For this model, an existence result is obtained in [8] and homogenization results in [6].

Let us note that all the aforementioned homogenization works are restricted to the case where the wetting phase (water) is incompressible while the non-wetting phase (gas) is compressible, contrarily to the present work. In this paper we extend our previous results obtained in [9] to the more complex case where both phases are compressible which is more reasonable in gas reservoir engineering. The major difficulties related to this model are in the nonlinear degenerate structure of the equations, as well as in the coupling in the system. In this case a new nonlinearity appears in the temporal term of the saturation equation. The compactness result used in [9] is no longer valid. To obtain these results we elaborated a new approach based on the ideas from [24, 35] to establish a new compactness result and uniform a priori bounds for the modulus of continuity with respect to the space and time variables.

In this paper, we will be concerned with a degenerate nonlinear system of diffusion-convection equations in a periodic domain modeling the flow and trans- port of immiscible compressible fluids through heterogeneous porous media, taking into account capillary and gravity effects. We consider double porosity media, i.e.

we consider a porous medium made up of a set of porous blocks with permeability
of order ε^{2} surrounded by a system of connected fissures, ε, is a small parameter
which characterizes the periodicity of the blocks. There are two kinds of degener-
acy in the studied system. The first one is the classical degeneracy of the capillary
diffusion term and the second one represents the evolution terms degeneracy. In
both cases the presence of degeneracy weakens the energy estimates and makes a
proof of compactness results more involved.

The outline of the rest of the paper is as follows. In Section 2 we describe the physical model and formulate the corresponding mathematical problem. We also provide the assumptions on the data and a weak formulation of the problem in terms of the global pressure and the saturation. Section 3 is devoted to the presentation of some a priori estimates for the solutions of the problem. They are essentially based on an energy equality. In Section 4, firstly we construct the extensions of the saturation and the global pressure functions defined in the fissures system and secondly we prove a compactness result adapted to our model. It’s based on the compactness criterion of Kolmogorov–Riesz–Fr´echet (see, e.g., [24, 35]). Finally, we formulate the corresponding two–scale convergence results. In Section 5 we are dealing with the dilations of the functions defined in the matrix part. Firstly, we introduce the notion of the dilation operator and describe its properties. Secondly, we derive the system of equations for the dilated functions and obtain the corresponding uniform estimates for them. Finally, we formulate the convergence results for the dilated functions. The main result of the paper is formulated in Section 6 and its proof is given in Section 7. The proof is based on the two-scale convergence and the dilation techniques.

2. Formulation of the problem

The outline of the section is as follows. First, in subsection 2.1 we present the model equations which are valid in fractures and rock matrix. A fractional flow formulation using the notion of the global pressure is discussed in subsection 2.2.

Then in the last subsection 2.3, we give the definition of a weak solution to our system.

2.1. Microscopic model. We consider a reservoir Ω ⊂ R^{d} (d = 2,3) which is
assumed to be a bounded, connected Lipschitz domain with a periodic microstruc-
ture. More precisely, we will scale this periodic structure by a parameter εwhich
represents the ratio of the cell size to the whole region Ω and we assume that
0< ε 1 is a small parameter tending to zero. Let Y := (0,1)^{d} be a periodicity
cell; we paveR^{d} withY. We assume thatYmis an open set with piecewise smooth
boundary ∂Ym such that Ym bY and we reproduceYm by periodicity, obtaining
a periodic open setM inR^{d}. We denote byFthe periodic setF:=R^{d}\M, which
is obtained from the set Yf := Y \Ym. Thus Y = Ym∪Yf ∪Γf m, where Γf m

denotes the interface between the two media. Finally, we denote byχf andχmthe
characteristic functions of the setsFandM. Thenχm(^{x}_{ε}) is the periodic function
of periodεY which takes the value 1 in the setM^{ε}, union of the sets obtained from
εYmby translations of vectors εPn

i=1ki~ei, whereki ∈Zand~ei, 1≤i≤d, is the
canonical basis ofR^{d}, and which takes the value 0 in the setF^{ε}, complementary in
R^{d}of this union. In other words,χ_{m}(^{x}_{ε}) is the characteristic function of the setM^{ε},
whileχ_{f}(^{x}_{ε}) is the characteristic function ofF^{ε}. Now we can define the subdomains
Ω^{ε}_{r}with r= “f^{00} or “m^{00} corresponding to the porous medium with the index “r^{00}.
We set:

Ω^{ε}_{m}:={x∈Ω :χ^{ε}_{m}(x) = 1} and Ω^{ε}_{f} :=

x∈Ω :χ^{ε}_{f}(x) = 1 .

Then Ω = Ω^{ε}_{m}∪Γ^{ε}_{f m}∪Ω^{ε}_{f}, where Γ^{ε}_{f m} := ∂Ω^{ε}_{f} ∩∂Ω^{ε}_{m}∩Ω and the subscript m
andf refer to the matrix and fracture, respectively. For the sake of simplicity, we
assume that Ω^{ε}_{m}∩∂Ω =∅. We also set:

Ω_{T} := Ω×(0, T), Ω^{ε}_{r,T} := Ω^{ε}_{r}×(0, T), and Γ^{ε,T}_{f,m}:= Γ^{ε}_{f,m}×(0, T), (2.1)
whereT >0 is fixed.

We consider an immiscible compressible two-phase flow system in a porous
medium which fills the domain Ω. We focus here on the general case where both
phases are compressible, the phases being ` and g. Let Φ^{ε}(x) be the porosity of
the reservoir Ω; K^{ε}(x) be the absolute permeability tensor of Ω; S_{`}^{ε} = S_{`}^{ε}(x, t),
S_{g}^{ε}=S_{g}^{ε}(x, t) be the phase saturations; k_{r,`}=k_{r,`}(S_{`}^{ε}),k_{r,g}=k_{r,g}(S_{g}^{ε}) be the rela-
tive permeabilities of the phases;p^{ε}_{`} =p^{ε}_{`}(x, t),p^{ε}_{g}=p^{ε}_{g}(x, t) be the phase pressures;

ρ`, ρg be the phase densities andPc the capillary pressure.

In what follows, for the sake of presentation simplicity we neglect the source
terms, and we denote S^{ε}=S_{`}^{ε}. The model for the two-phase flow is described by

(see, e.g., [27, 31, 43]):

06S^{ε}61 in Ω_{T};
Φ^{ε}(x)∂Ξ^{ε}_{`}

∂t −divn

K^{ε}(x)λ`(S^{ε})ρ`(p^{ε}_{`}) ∇p^{ε}_{`}−ρ`(p^{ε}_{`})~go

= 0 in ΩT;
Φ^{ε}(x)∂Ξ^{ε}_{g}

∂t −divn

K^{ε}(x)λg(S^{ε})ρg(p^{ε}_{g}) ∇p^{ε}_{g}−ρg(p^{ε}_{g})~go

= 0 in ΩT;
P_{c}(S^{ε}) =p^{ε}_{g}−p^{ε}_{`} in Ω_{T},

(2.2)

whereλ_{g}(S^{ε}) =eλ_{g}(1−S^{ε}); Ξ^{ε}_{`} :=S^{ε}ρ_{`}(p^{ε}_{`}) and Ξ^{ε}_{g}:= (1−S^{ε})ρ_{g}(p^{ε}_{g}); each function
γ^{ε}:=S^{ε}, p_{`}, p_{g}, Ξ^{ε}_{`}, and Ξ^{ε}_{g} is defined as:

γ^{ε}(x, t) =χ^{ε}_{f}(x)γ_{f}^{ε}(x, t) +χ^{ε}_{m}(x)γ^{ε}_{m}(x, t). (2.3)
The velocities of the phases~q_{`}^{ε},~q^{ε}_{g} are defined by Darcy–Muskat’s law:

~q_{`}^{ε}:=−K^{ε}(x)λ`(S^{ε}_{`})

∇p^{ε}_{`}−ρ`(p^{ε}_{`})~g

withλ`(S_{`}^{ε}) := k_{r,`}

µ`

(S_{`}^{ε}); (2.4)

−

→q^{ε}_{g} :=−K^{ε}(x)eλ_{g}(S_{g}^{ε})

∇p^{ε}_{g}−ρ_{g}(p^{ε}_{g})~g

witheλ_{g}(S_{g}^{ε}) :=kr,g

µ_{g} (S^{ε}_{g}) (2.5)
with~g,µ`,µg being the gravity vector and the viscosities, respectively.

Now we specify the boundary and initial conditions. We suppose that the bound-
ary∂Ω consists of two parts Γ_{1} and Γ_{2} such that Γ_{1}∩Γ_{2}=∅,∂Ω = Γ_{1}∪Γ_{2}. The
boundary conditions are given by

p^{ε}_{g}(x, t) = 0 =p^{ε}_{`}(x, t) on Γ1×(0, T); (2.6)

~

q_{`}^{ε}· −→ν =−→

q_{g}^{ε}· −→ν = 0 on Γ_{2}×(0, T). (2.7)
Finally, the initial conditions read

S^{ε}(x,0) =S^{0}(x) and p^{ε}_{g}(x,0) =p^{0}_{g}(x) in Ω. (2.8)
2.2. A fractional flow formulation. In the sequel, we use a formulation obtained
after transformation using the concept of the global pressure introduced in [13, 27].

The global pressure is defined as follows:

P^{ε}=p^{ε}_{`}−G`(S^{ε}) =p^{ε}_{g}−Gg(S^{ε}), (2.9)
where the functionsG`(S^{ε}),Gg(S^{ε}) are given by:

Gg(S^{ε}) :=Gg(0) +
Z S^{ε}

0

λ`(s)

λ(s)P_{c}^{0}(s)ds, (2.10)
G`(S^{ε}) =Gg(S^{ε})−Pc(S^{ε}), (2.11)
withλ(s) :=λ`(s) +λg(s), the total mobility.

Performing some simple calculations, we obtain the following properties for the global pressure which will be used in the sequel:

λ`(S^{ε})∇p^{ε}_{`}+λg(S^{ε})∇p^{ε}_{g}=λ(S^{ε})∇P^{ε}, (2.12)

∇G`(S^{ε}) =−λ_{g}(S^{ε})

λ(S^{ε})P_{c}^{0}(S^{ε})∇S^{ε}. (2.13)
Notice that from (2.10), (2.13) we obtain

λ_{`}(S^{ε})∇G`(S^{ε}) =∇β(S^{ε}) and λ_{g}(S^{ε})∇Gg(S^{ε}) =−∇β(S^{ε}), (2.14)

where

β(S^{ε}) :=

Z S^{ε}

0

α(u)du withα(s) :=λg(s)λ`(s)

λ(s) |P_{c}^{0}(s)|. (2.15)
Furthermore, we have the important relation

λg(S^{ε})|∇p^{ε}_{g}|^{2}+λ`(S^{ε})|∇p^{ε}_{`}|^{2}=λ(S^{ε})|∇P^{ε}|^{2}+|∇b(S^{ε})|^{2}, (2.16)
where

b(s) :=

Z s

0

a(ξ)dξ witha(s) :=

s

λ_{g}(s)λ_{`}(s)

λ(s) |P_{c}^{0}(s)|. (2.17)
If we use the global pressure and the saturation as new unknown functions, then
problem (2.2) reads

06S^{ε}61 in Ω_{T};
Φ^{ε}(x)∂Θ^{ε}_{`}

∂t −divn

K^{ε}(x)ρe^{ε}_{`}

λ`(S^{ε})∇P^{ε}+∇β(S^{ε})−λ`(S^{ε})ρe^{ε}_{`}~go

= 0 in ΩT;
Φ^{ε}(x)∂Θ^{ε}_{g}

∂t −divn

K^{ε}(x)ρe^{ε}_{g}

λg(S^{ε})∇P^{ε}− ∇β(S^{ε})−λg(S^{ε})ρegε

~ go

= 0 in ΩT, (2.18) we introduced the notation

ρe^{ε}_{`} :=ρ`(P^{ε}+G`(S^{ε})) and ρe^{ε}_{g}:=ρg(P^{ε}+Gg(S^{ε})); (2.19)
Θ^{ε}_{`} = Θ`(S^{ε}, P^{}) :=S^{ε}ρe^{ε}_{`} and Θ^{ε}_{g}= Θg(S^{ε}, P^{ε}) := (1−S^{ε})ρe^{ε}_{g}. (2.20)
The system (2.18) is completed by the following boundary and initial conditions.

P^{ε}= 0 on Γ1×(0, T); (2.21)

Q^{ε}_{`}· −→ν =Q^{ε}_{g}· −→ν on Γ2×(0, T) (2.22)
whereQ^{ε}_{`} andQ^{ε}_{g} are defined by:

Q^{ε}_{`}:=−K^{ε}(x)ρe^{ε}_{`}

λ_{`}(S^{ε})∇P^{ε}+∇β(S^{ε})−λ_{`}(S^{ε})ρe^{ε}_{`}~g
,
Q^{ε}_{g}:=−K^{ε}(x)ρe^{ε}_{g}

λg(S^{ε})∇P^{ε}− ∇β(S^{ε})−λg(S^{ε})ρe^{ε}_{g}~g
.
Finally, the initial conditions read

S^{ε}(x,0) =S^{0}(x) and P^{ε}(x,0) =P^{0}(x) in Ω. (2.23)
2.3. A weak formulation of the problem. Let us begin this section by stating
the following assumptions.

(A1) The porosity Φ = Φ(y) is a Y-periodic function defined by: Φ^{ε}(x) =
χ^{ε}_{f}(x)Φf +χ^{ε}_{m}(x)Φm with 0 < Φf,Φm < 1, where Φf and Φm are con-
stant that do not depend onε.

(A2) The absolute permeability tensorK^{ε} is given by:

K^{ε}(x) :=Kfχ^{ε}_{f}(x)I+ε^{2}Kmχ^{ε}_{m}(x)I,

where Iis the unit tensor andKf,Km are positive constants that do not depend onε.

(A3) The density ρ_{k} =ρ_{k}(p), (k = `, g) is a monotone C^{1}-function in R such
that

ρk(p) =ρmin forp6pmin; ρk(p) =ρmax forp>pmax;

ρ_{min}< ρ_{k}(p)< ρ_{max} forp_{min}< p < p_{max}. (2.24)

ρmin,ρmax,pmin,pmax are constants such that 0< ρmin< ρmax<+∞and 0< pmin< pmax<+∞.

(A4) The capillary pressure functionPc∈C^{1}([0,1];R^{+}). Moreover,P_{c}^{0}(s)<0 in
[0,1] andP_{c}(1) = 0.

(A5) The functions λ_{`}, λ_{g} belong to the space C([0,1];R^{+}) and satisfy the fol-
lowing properties:

(i) 06λ_{`}, λ_{g}61 in [0,1];

(ii) λ_{`}(0) = 0 andλ_{g}(1) = 0;

(iii) there is a positive constantL0such thatλ`(s) =λ`(s)+λg(s)>L0>0 in [0,1].

Moreover,λ`(s)∼s^{κ}^{`} ass→0 andλg(s)∼(1−s)^{κ}^{g} ass→1 (κ`,κg>0).

(A6) The functionαgiven by (2.15) is a continuous function in [0,1].Moreover, α(0) =α(1) = 0 andα >0 in (0,1).

(A7) The function β^{−1}, inverse of β defined in (2.15) is a H¨older function of
order θ with θ ∈ (0,1) on the interval [0, β(1)]. Namely, there exists a
positive constantCβsuch that for alls1, s2∈[0, β(1)],|β^{−1}(s1)−β^{−1}(s2)|6
Cβ|s1−s2|^{θ}.

(A8) The initial data for the global pressure and the saturation defined in (2.23)
are such thatP^{0}∈L^{2}(Ω) and 06S^{0}61.

Assumptions (A1)–(A8) are classical and physically meaningful for existence results and homogenization problems of two-phase flow in porous media. They are similar to the assumptions made in [10] that dealt with the existence of a weak solution of the studied problem.

Next we introduce the Sobolev space

H_{Γ}^{1}_{1} :={u∈H^{1}(Ω) :u= 0 on Γ1},
wich is a Hilbert space when it is equipped with the norm

kuk_{H}^{1}

Γ1(Ω) =k∇uk_{(L}2(Ω))^{d}.

Definition 2.1. We say that the pair of functionshP^{ε}, S^{ε}i is a weak solution to
problem (2.18)–(2.23) if

(i) 06S^{ε}61 a.e in ΩT.
(ii) P^{ε}∈L^{2}(0, T;H_{Γ}^{1}

1(Ω)).

(iii) The boundary conditions (2.21)–(2.22) are satisfied.

(iv) For anyϕ`, ϕg ∈C^{1}([0, T];H_{Γ}^{1}

1(Ω)) satisfyingϕ`(T) =ϕg(T) = 0, we have

− Z

Ω_{T}

Φ^{ε}(x)Θ^{ε}_{`}∂ϕ^{ε}_{`}

∂t dx dt+ Z

Ω

Φ^{ε}(x)Θ^{0}_{`}ϕ^{0}_{`}dx
+

Z

ΩT

K^{ε}(x)ρe^{ε}_{`}n

λ`(S^{ε}) ∇P^{ε}−ρe^{ε}_{`}~g

+∇β(S^{ε})o

· ∇ϕ^{ε}_{`} dx dt= 0;

(2.25)

− Z

ΩT

Φ^{ε}(x)Θ^{ε}_{g}∂ϕ^{ε}_{g}

∂t dx dt+ Z

Ω

Φ^{ε}(x)Θ^{0}_{g}ϕ^{0}_{g}dx
+

Z

Ω_{T}

K^{ε}(x)ρe^{ε}_{g}n

λg(S^{ε}) ∇P^{ε}−ρe^{ε}_{g}~g

− ∇β(S^{ε})o

· ∇ϕ^{ε}_{`} dx dt= 0,

(2.26)

where ρe^{ε}_{g} andρe^{ε}_{`} are defined in (2.19);ϕ^{0}_{`} =ϕ_{`}(0, x), ϕ^{0}_{g} =ϕ_{g}(0, x); Θ^{0}_{`} =
S^{0}ρ`(P^{0}+G`(S^{0})) and Θ^{0}_{g}= (1−S^{0})ρg(P^{0}+Gg(S^{0})).

According to [10], under conditions (A1)–(A8), for eachε >0, problem (2.25)–

(2.26) has at least one weak solution.

In what followsC, C1, . . . denote generic constants that do not depend onε.

3. A priori estimates

To obtain the needed uniform estimates for the solution of problem (2.2) (or the equivalent problem (2.18)), we follow the choice of the test functions as in [9]:

R`(p^{ε}_{`}) =
Z p^{ε}_{`}

0

dξ

ρ`(ξ) and Rg(p^{ε}_{g}) =
Z p^{ε}_{g}

0

dξ

ρg(ξ). (3.1) Then, as in [9], the following results hold.

Lemma 3.1. Let hp^{ε}_{g}, p^{ε}_{`}ibe a solution to (2.2). Then we have the energy equality
d

dt Z

Ω

Φ^{ε}(x)ζ^{ε}(x, t)dx+
Z

Ω

K^{ε}(x)n

λ`(S^{ε})∇p^{ε}_{`}·

∇p^{ε}_{`}−ρ`(p^{ε}_{`})~g
+λg(S^{ε})∇p^{ε}_{g}·

∇p^{ε}_{g}−ρg(p^{ε}_{g})~go
dx= 0

(3.2) in the sense of distributions. Here

ζ^{ε}:=S^{ε}R`(p^{ε}_{`}) + (1−S^{ε})Rg(p^{ε}_{g}) +F(S^{ε}),
where Rk(p) := ρk(p)Rk(p)−p, (k = `, g) and F(s) := Rs

1 Pc(u)du. Moreover,
ζ^{ε}>0 inΩT.

Lemma 3.2. Let hp^{ε}_{g}, p^{ε}_{`}ibe a solution to (2.2). Then
kp

K^{ε}(x)λ`(S^{ε})∇p^{ε}_{`}k_{L}2(ΩT)+kq

K^{ε}(x)λg(S^{ε})∇p^{ε}_{g}k_{L}2(ΩT)6C. (3.3)
Corollary 3.3. Let hp^{ε}_{g}, p^{ε}_{`}ibe a solution to (2.2). Then

kq

λ_{`}(S_{f}^{ε})∇p^{ε}_{`,f}k_{L}2(Ω^{ε}_{f,T})+kq

λ_{g}(S_{f}^{ε})∇p^{ε}_{g,f}k_{L}2(Ω^{ε}_{f,T})

+εkp

λ`(S_{m}^{ε})∇p^{ε}_{`,m}k_{L}2(Ω^{ε}_{m,T})+εkq

λg(S_{m}^{ε})∇p^{ε}_{g,m}k_{L}2(Ω^{ε}_{m,T})6C.

(3.4)
Then we obtain the following uniform a priori estimates for the functions P^{ε}
andβ(S^{ε}).

Lemma 3.4. Let the pair of functionshP^{ε}, S^{ε}ibe a solution to (2.18). Then
k∇β(S_{f}^{ε})k_{L}2(Ω^{ε}_{f,T})+k∇P_{f}^{ε}k_{L}2(Ω^{ε}_{f,T})+εk∇β(S_{m}^{ε})k_{L}2(Ω^{ε}_{m,T})

+εk∇P_{m}^{ε}k_{L}2(Ω^{ε}_{m,T})6C. (3.5)

Moreover,

kP_{f}^{ε}k_{L}2(Ω^{ε}_{f,T})+kβ(S^{ε})k_{L}2(ΩT)6C, (3.6)
kP_{m}^{ε}k_{L}2(Ω^{ε}_{m,T})6C. (3.7)
Now we pass to the uniform bounds for the time derivatives of the functions Θ^{ε}_{g}
andS^{ε}. In a standard way (see, e.g., [4, 9]) we can prove the following lemma.

Lemma 3.5. Let the pair of functions hP^{ε}, S^{ε}i be a solution to (2.18). Then for
r=f, m,

{∂_{t}(Φ_{r}Θ^{ε}_{`,r})}_{ε>0} is uniformly bounded in L^{2}(0, T;H^{−1}(Ω^{ε}_{r})); (3.8)
{∂_{t}(Φ_{r}Θ^{ε}_{g,r})}_{ε>0} is uniformly bounded inL^{2}(0, T;H^{−1}(Ω^{ε}_{r})). (3.9)

4. Convergence of{P_{f}^{ε}}ε>0, {S_{f}^{ε}}ε>0, {Θ^{ε}_{`,f}}ε>0, {Θ^{ε}_{g,f}}ε>0

In this section, we obtain compactness results that will be used in passing to
the limit asεtends to zero in the weak formulation. The compactness result used
in [9] is no longer valid. To obtain these results we elaborated a new approach
based on the ideas from [24, 35] to establish a new compactness result and uniform
a priori bounds for the modulus of continuity with respect to the space and time
variables. It is achieved in several steps. First, in subsection 4.1 we extend the
functionS_{f}^{ε} from the subdomain Ω^{ε}_{f} to the whole Ω and obtain uniform estimates
for the extended function Se_{f}^{ε}. Then in Section 4.2, using the uniform estimates
for the functionPe_{f}^{ε} which follow from Lemma 3.4, the definition of the extension
operator, and the corresponding bounds for Se_{f}^{ε}, we prove the compactness result
for the families {Θe^{ε}_{`,f}}ε>0 and {Θe^{ε}_{g,f}}ε>0, where Θe^{ε}_{`,f}, Θe^{ε}_{g,f} are extensions of the
functions Θ^{ε}_{`,f}, Θ^{ε}_{g,f} from the subdomain Ω^{ε}_{f} to the whole Ω which will be specified
at the end of subsection 4.1. Finally, in subsection 4.3 we formulate the two–scale
convergence which will be used in the derivation of the homogenized system.

4.1. Extensions of the functionsS_{f}^{ε},Θ^{ε}_{`,f}, Θ^{ε}_{g,f}. Consider the functionS_{f}^{ε}. To
extendS_{f}^{ε}, following the ideas of [23], we use the monotone function β defined in
(2.15). Let us introduce the function

β_{f}^{ε}(x, t) :=β(S_{f}^{ε}) =
Z S^{ε}_{f}

0

α(u)du. (4.1)

Then it follows from condition (A6) that 0 6 β^{ε}_{f} 6 max_{s∈[0,1]}α(s) a.e. in Ω^{ε}_{f,T}.
Furthermore, from (3.5) we have

k∇β_{f}^{ε}kL^{2}(Ω^{ε}_{f,T})6C. (4.2)
Let Π^{ε} : H^{1}(Ω^{ε}_{f})→ H^{1}(Ω) be the standard extension operator cf. [1]. Then we
have

06βe^{ε}_{f} := Π^{ε}β^{ε}_{f}6 max

s∈[0,1]α(s) a.e. in ΩT,k∇eβ_{f}^{ε}kL^{2}(Ω_{T})6C.

Now we can extend the function S_{f}^{ε} from the subdomain Ω^{ε}_{f} to the whole Ω. We
denote this extension bySe_{f}^{ε} and define it as:

Se_{f}^{ε}:= (β)^{−1}(eβ_{f}^{ε}).

This implies that Z

ΩT

|∇β(Se_{f}^{ε})|^{2}dx dt6C ,06Se_{f}^{ε}61 a.e. in ΩT. (4.3)
Finally consider the sequences{Θ^{ε}_{`,f}}ε>0and{Θ^{ε}_{g,f}}ε>0. We recall that Θ^{ε}_{`,f} :=

ρ` P_{f}^{ε}+G`(S_{f}^{ε})

S_{f}^{ε} and Θ^{ε}_{g,f} := ρg P_{f}^{ε}+Gg(S_{f}^{ε})

(1−S_{f}^{ε}). Then we define the
extension of the function Θ^{ε}_{f} to the whole Ω by

Θe^{ε}_{`,f} :=ρ_{`}(Pe_{f}^{ε}+G_{`}(Se_{f}^{ε}))eS_{f}^{ε} and Θe^{ε}_{g,f}:=ρ_{g}(Pe_{f}^{ε}+G_{g}(Se_{f}^{ε}))(1−Se_{f}^{ε}), (4.4)
wherePe_{f}^{ε}:= Π^{ε}P_{f}^{ε}is the extension of the functionP_{f}^{ε}.

4.2. Compactness of the sequences{Θe^{ε}_{`,f}}ε>0, {Θe^{ε}_{g,f}}ε>0. The following con-
vergence result is valid.

Proposition 4.1. Under our standing assumptions there exist the functionsL1,L2

such that

Θe^{ε}_{`,f} →L1 strongly inL^{2}(ΩT), (4.5)
Θe^{ε}_{g,f}→L_{2} strongly inL^{2}(Ω_{T}). (4.6)
Proof. We will prove the convergence result for the sequence{Θe^{ε}_{`,f}}ε>0, the corre-
sponding result for the sequence{Θe^{ε}_{g,f}}_{ε>0} can be obtained by similar arguments.

The scheme of the proof is as follows. We apply the compactness criterion of
Kolmogorov-Riesz-Fr´echet (see, e.g., [24, 35]) in the spaceL^{1}(Ω_{T}). To this end we
have to obtain the moduli of continuity with respect to the space and temporal
variables (see Lemmata 4.2, 4.6 below). Finally, the uniform boundedness of the
function Θ^{ε}_{`,f} imply the desired convergence result (4.5) in the spaceL^{2}(ΩT).

We start with the following result which can be proved by arguments similar to those from Lemma 4.2 in [9].

Lemma 4.2 (Modulus of continuity with respect to the space variable). Let con- ditions(A1)–(A8) be fulfilled. Then for|∆x| sufficiently small,

Z T

0

Z

Ω

eΘ^{ε}_{`,f}(x+ ∆x, t)−Θe^{ε}_{`,f}(x, t)

2dx dt6C|∆x|^{θ}, (4.7)
whereθ∈(0,1)is defined in condition (A7).

Now we turn to the derivation of the modulus of continuity with respect to the temporal variable. To do this, for anyδ >0, we introduce the functions

S_{r}^{ε,δ}:= min

1−δ,max(δ, S_{r}^{ε}) withr=f, m.

Let us estimate the norm of the functionS_{r}^{ε,δ} (r=f, m) in L^{2}(0, T;H^{1}(Ω^{ε}_{r,T})).

Lemma 4.3. Under our standing assumptions,

kS_{f}^{ε,δ}k_{L}2(0,T;H^{1}(Ω^{ε}_{f,T}))+εkS_{m}^{ε,δ}k_{L}2(0,T;H^{1}(Ω^{ε}_{m,T}))6C δ^{−η}, (4.8)
whereη:= (κ`+κg)andC is a constant that does not depend on ε, δ.

Proof. We consider the function S_{f}^{ε,δ}, the estimate for the function S_{m}^{ε,δ} can be
obtained in a similar way. It is evident that S_{f}^{ε,δ} (as well as S_{f}^{ε}) belongs to the
spaceL^{2}(Ω^{ε}_{f,T}). Now we are going to estimate ∇S^{ε,δ}_{f} . From Lemma 3.4 we know
that

k∇β(S_{f}^{ε})k_{L}2(Ω^{ε}_{f,T})6C.

Then it is clear that

k∇β(S_{f}^{ε,δ})k_{L}2(Ω^{ε}_{f,T})6k∇β(S_{f}^{ε})k_{L}2(Ω^{ε}_{f,T})6C, (4.9)
where C is a constant that does not depend on ε, δ. Moreover, condition (A5)
implies the inequalities:

λ`(S_{f}^{ε,δ})>C1δ^{κ}^{`} and λg(S_{f}^{ε,δ})>C2δ^{κ}^{g}. (4.10)

Then taking into account the definition of the functions β and α(see (2.15)) and conditions (A4), (A5), from (4.10), we have that

α(S_{f}^{ε,δ})>C1δ^{(}^{κ}^{`}^{+}^{κ}^{g}^{)} withC1:= C_{1}C_{2}

2 min

s∈[0,1]|P_{c}^{0}(s)|. (4.11)
Then from (4.9), (4.11)

C>k∇β(S_{f}^{ε,δ})k^{2}_{L}2(Ω^{ε}_{f,T})=
Z T

0

Z

Ω^{ε}_{f}

α^{2}(S_{f}^{ε,δ})|∇S^{ε,δ}_{f} |^{2}dx dt

>C^{2}1δ^{2(}^{κ}^{`}^{+}^{κ}^{g}^{)}
Z T

0

Z

Ω^{ε}_{f}

|∇S_{f}^{ε,δ}|^{2}dx dt

=C^{2}1δ^{2(}^{κ}^{`}^{+}^{κ}^{g}^{)}k∇S_{f}^{ε,δ}k^{2}_{L}2(Ω^{ε}_{f,T}).

(4.12)

The inequality (4.11) implies that

k∇S_{f}^{ε,δ}k^{2}_{L}2(Ω^{ε}_{f,T})6C δ^{−2(}^{κ}^{`}^{+}^{κ}^{g}^{)} (4.13)
and inequality (4.8) is proved. This completes the proof of Lemma 4.3.

In what follows we use a technical lemma which can be proved using the Fubini theorem.

Lemma 4.4. For h sufficiently small, 0 < h < ^{T}_{2} and for integrable functions
G1(t),G2(t)it holds:

Z T

0

G1(t)Z min(t+h,T) max(t,h)

G2(τ)dτ dt=

Z T

h

G2(t)Z t t−h

G1(τ)dτ dt.

Now, forε >0 and 0< h < ^{T}_{2}, let us introduce the function
ϕ^{ε,δ,h}(x, t) :=

Z min(t+h,T)

max(t,h)

h∂^{−h}Θ^{ε,δ}_{`} (x, τ)dτ,
with ∂^{−h}v(t) := v(t)−v(t−h)

h ,

(4.14)

where

Θ^{ε,δ}_{`} (x, t) :=ρ`(P^{ε}+G`(S^{ε,δ}))S^{ε,δ}. (4.15)
The properties ofϕ^{ε,δ,h}are described by the next lemma.

Lemma 4.5. Let ε > 0, 0 < δ <1, and let h >0 be small enough. There exist a constant C which does not depend on ε, δ, and h such that for the sequence of functions defined by (4.14) it holds

ϕ^{ε,δ,h}∈L^{2}(0, T;H_{Γ}^{1}_{1}(Ω)); (4.16)

ϕ^{ε,δ,h}(x, T) = 0; (4.17)

kϕ^{ε,δ,h}k_{L}2(ΩT)6C h; (4.18)

k∇ϕ^{ε,δ,h}kL^{2}(Ω^{ε}_{f,T})6C h δ^{−η}; (4.19)
εk∇ϕ^{ε,δ,h}k_{L}2(Ω^{ε}_{m,T})6C h δ^{−η}. (4.20)
Here

η:= (κ`+κg). (4.21)

Proof. The regularity property (4.16) follows immediately from Lemma 3.4 and
Lemma 4.3. Moreover, taking into account thatS^{ε} = 1 and P^{ε}= const on Γ1×
(0, T) (see Section 2), according to the definition of the functionS^{ε,δ}, we have that
ϕ^{ε,δ,h}= 0 on Γ_{1}×(0, T).

Result (4.17) follows directly from the definition of the functionϕ^{ε,δ,h}. In fact,
ϕ^{ε,δ,h}(x, T) =

Z min(T+h,T)

max(T ,h)

h∂^{−h}Θ^{ε,δ}_{`} (x, τ)dτ =
Z T

T

h ∂^{−h}Θ^{ε,δ}_{`} (x, τ)dτ = 0.

Bound (4.18) also follows immediately from the definition ofϕ^{ε,δ,h}since min(t+
h, T)−max(t, h)6hand the function Θ^{ε,δ}_{`} is uniformly bounded inL^{∞}(ΩT).

For bound (4.19), we have
k∇ϕ^{ε,δ,h}k^{2}_{L}2(Ω^{ε}_{f,T})

6 Z

Ω^{ε}_{f,T}

hZ min(t+h,T)

max(t,h)

∇Θ^{ε,δ}_{`,f}(x, τ)− ∇Θ^{ε,δ}_{`,f}(x, τ−h)
dτi2

dx dt:=J^{ε,δ}. (4.22)
Let us estimate the right-hand side of (4.22). Since [min(t+h, T)−max(t, h)]6h,
from Cauchy’s inequality, for a.e. (x, t)∈Ω^{ε}_{f,T} we obtain

Z min(t+h,T)

max(t,h)

∇Θ^{ε,δ}_{`,f}(x, τ)− ∇Θ^{ε,δ}_{`,f}(x, τ−h)
dτ

6h^{1/2}hZ min(t+h,T)
max(t,h)

∇Θ^{ε,δ}_{`,f}(x, τ)− ∇Θ^{ε,δ}_{`,f}(x, τ−h)

2dτi1/2

.

Therefore, from this inequality we obtain
J^{ε,δ}6C h

Z T

0

Z

Ω^{ε}_{f}

hZ min(t+h,T)

max(t,h)

∇Θ^{ε,δ}_{`,f}(x, τ)− ∇Θ^{ε,δ}_{`,f}(x, τ−h)

2dτi

dx dt. (4.23)
Now we apply Lemma 4.4 withG1(t) := 1 andG2(t) :=|∇Θ^{ε,δ}_{`,f}(x, τ)− ∇Θ^{ε,δ}_{`,f}(x, τ−
h)|^{2}in the right–hand side of (4.23). We have:

Z T

0

hZ min(t+h,T)

max(t,h)

∇Θ^{ε,δ}_{`,f}(x, τ)− ∇Θ^{ε,δ}_{`,f}(x, τ−h)

2dτi dt

= Z T

h

∇Θ^{ε,δ}_{`,f}(x, t)− ∇Θ^{ε,δ}_{`,f}(x, t−h)

2hZ t

t−h

1dτi dt

=h Z T

h

∇Θ^{ε,δ}_{`,f}(x, t)− ∇Θ^{ε,δ}_{`,f}(x, t−h)

2dt.

Then, by (4.23), we deduce that
J^{ε,δ}

6C h^{2}
Z

Ω^{ε}_{f}

Z T

h

∇Θ^{ε,δ}_{`,f}(x, t)− ∇Θ^{ε,δ}_{`,f}(x, t−h)

2dt dx

6Ch^{2}hZ

Ω^{ε}_{f}

Z T

h

∇Θ^{ε,δ}_{`,f}(x, t)

2dt dx+ Z

Ω^{ε}_{f}

Z T

h

∇Θ^{ε,δ}_{`,f}(x, t−h)

2dt dxi

6C h^{2}k∇Θ^{ε,δ}_{`,f}k^{2}_{L}2(Ω^{ε}_{f,T}).

(4.24)

It remains to estimate the right-hand side of (4.24). To this end we rewrite

∇Θ^{ε,δ}_{`,f} as follows:

∇Θ^{ε,δ}_{`,f} =ρ_{`} P^{ε}+G_{`}(S^{ε,δ})

∇S^{ε,δ}+ρ^{0}_{`}S^{ε,δ}∇P^{ε}+ρ^{0}_{`}S^{ε,δ}∇G_{`}(S^{ε,δ}). (4.25)
Then, from (4.25), conditions (A3), (A5), and the definition of the functionG`, we
have

k∇Θ^{ε,δ}_{`,f}k^{2}_{L}2(Ω^{ε}_{f,T})

6ρ^{2}_{max}k∇S^{ε,δ}k^{2}_{L}2(Ω^{ε}_{f,T})+Ch

k∇P^{ε}k^{2}_{L}2(Ω^{ε}_{f,T})+kP_{c}^{0}∇S^{ε,δ}k^{2}_{L}2(Ω^{ε}_{f,T})

i.

(4.26) To estimate the right-hand side of (4.26), we make use of condition (A4), bound (4.8), and Lemma 4.3. Taking into account thatδis sufficiently small, we obtain

k∇Θ^{ε,δ}_{`,f}k^{2}_{L}2(Ω^{ε}_{f,T})6C δ^{−2(}^{κ}^{`}^{+}^{κ}^{g}^{)}. (4.27)
Now, from (4.22), (4.24), and (4.27), forδ sufficiently small, we obtain

k∇ϕ^{ε,δ,h}k^{2}_{L}2(Ω^{ε}_{f,T})6C h^{2}δ^{−2(}^{κ}^{`}^{+}^{κ}^{g}^{)} (4.28)
and the bound (4.19) is proved.

The proof of the bound (4.20) can be done by arguments similar to ones used in the proof of (4.19). This completes the proof of Lemma 4.5.

Now we are in a position to estimate the modulus of continuity of the function
Θe^{ε}_{`,f}.

Lemma 4.6 (Modulus of continuity with respect to time ). Under our standing assumptions, for all h∈(0, T)andδ sufficiently small, we have

Z T

h

Z

Ω

eΘ^{ε}_{`,f}(x, t)−Θe^{ε}_{`,f}(x, t−h)

2dx dt6Ch^{σ} with σ:= min1
2, 1

2η , (4.29) whereη is defined in (4.21)andCis a constant which does not depend onεandh.

Proof. Let us insert the functionϕ^{ε,δ,h}=ϕ^{ε,δ,h}(x, t) in equation (2.25). We have
Z

ΩT

Φ^{ε}(x)∂Θ^{ε}_{`}

∂t ϕ^{ε,δ,h}dx dt

+ Z

Ω_{T}

K^{ε}(x)eρ^{ε}_{`}n

λ_{`}(S^{ε}) ∇P^{ε}−ρe^{ε}_{`}~g

+∇β(S^{ε})o

· ∇ϕ^{ε,δ,h}dx dt= 0.

(4.30)

Taking into account the definition of the functionϕ^{ε,δ,h}and condition (A1), for the
first term in the left-hand side of (4.30), we have

I_{1}^{ε}(ϕ^{ε,δ,h}) : =
Z

ΩT

Φ^{ε}(x)∂Θ^{ε}_{`}

∂t ϕ^{ε,δ,h}dx dt

= Z T

0

Z

Ω^{ε}_{f}

Φf

∂Θ^{ε}_{`,f}

∂t

hZ min(t+h,T)

max(t,h)

h∂^{−h}Θ^{ε,δ}_{`,f}(x, τ)dτi
dx dt

+ Z T

0

Z

Ω^{ε}_{m}

Φm

∂Θ^{ε}_{`,m}

∂t

hZ min(t+h,T)

max(t,h)

h∂^{−h}Θ^{ε,δ}_{`,m}(x, τ)dτi
dx dt.