\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{amssymb} \usepackage{euscript} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2016 (2016), No. 52, pp. 1--28.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2016 Texas State University.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2016/52\hfil Immiscible compressible two-phase flow] {Homogenization of immiscible compressible two-phase flow in double porosity media} \author[L. Ait Mahiout, B. Amaziane, A. Mokrane, L. Pankratov \hfil EJDE-2016/52\hfilneg] {Latifa Ait Mahiout, Brahim Amaziane, \\ Abdelhafid Mokrane, Leonid Pankratov} \address{Latifa Ait Mahiout \newline LEDP, \'Ecole Normale Sup\'erieure de Kouba, BP 92, 16050, Vieux Kouba, Alger, Algeria} \email{latifaaitmahiout@gmail.com} \address{Brahim Amaziane \newline Laboratoire de Math\'ematiques et de leurs Applications, CNRS-UMR 5142, Universit\'e de Pau, Av. de l'Universit\'e, 64000 Pau, France} \email{brahim.amaziane@univ-pau.fr} \address{Abdelhafid Mokrane \newline LEDP, \'Ecole Normale Sup\'erieure de Kouba, BP 92, 16050, Vieux Kouba, Alger, Algeria} \email{abdelhafid.mokrane@gmail.com} \address{Leonid Pankratov \newline Laboratoire de Math\'ematiques et de leurs Applications, CNRS-UMR 5142, Universit\'e de Pau, Av. de l'Universit\'e, 64000 Pau, France.\newline Laboratory of Fluid Dynamics and Seismics, Moscow Institute of Physics and Technology, 9 Institutskiy per., Dolgoprudny, Moscow Region, 141700, Russian Federation} \email{leonid.pankratov@univ-pau.fr} \thanks{Submitted February 2, 2016. Published February 18, 2016.} \subjclass[2010]{35B27, 35K65, 76S05, 76T10} \keywords{Compressible immiscible; double porous media; two-phase flow; \hfill\break\indent fractured media homogenization; two-scale convergence} \begin{abstract} A double porosity model of multidimensional immiscible compressible 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 written 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 parabolic equation for the pressure and a doubly nonlinear degenerate parabolic diffusion-convection equation for the saturation, subject to appropriate boundary 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 temporal term of the saturation equation. The aim of this paper is to extend the results of \cite{ap-m2as} to this more general case. With the help of a new compactness result 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. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{remark}[theorem]{Remark} \newtheorem{definition}[theorem]{Definition} \newtheorem{corollary}[theorem]{Corollary} \allowdisplaybreaks \section{Introduction}\label{int} The modeling of displacement process involving two immiscible fluids is of considerable 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 $CO_2$. Furthermore, fractured rock domains corresponding to the so-called Excavation Damaged Zone (EDZ) receives increasing attention in connection with the behaviour of geological isolation of radioactive waste after the drilling of the wells or shafts, see, e.g., \cite{Shaw}. 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 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. \cite{bar}. 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., \cite{bear,panf,van}. During recent decades mathematical analysis and numerical simulation of multiphase flows in porous media have been the subject of investigation of many researchers 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 mathematical analysis of such models. We refer, for instance, to the books \cite{ant-kaz-mon1990,GC-JJ,Chen-Huan-Ma,GG-MM-1996,HR,hor,Vazquez} and the references therein. The mathematical analysis and the homogenization 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 \cite{alt-Bene,ant-kaz-mon1990,Arbogast-92,CC-PM,GC-JJ,Chen-I,Chen-II,GG-MM-1996, Kro-Luck} 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 \cite{Siam,HOS-2012,hor}. We refer for instance to \cite{Bourgeat-Gip2004,Bourgeat-Gip2010, Bourgeat-Marusic2005,Bourgeat-Marusic2010,Bourgeat-Piat2010,Gipouloux-Smai,Gloria-Etal} for more information on the homogenization of incompressible, 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 \cite{adh}, and was then revisited in the mathematical literature by many other authors. Here we restrict ourself to the mathematical homogenization method as described in \cite{hor} 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 \cite{APR,adh,blm,Chen,choq,GWC-RES,yeh2} and the references therein. However, as reported in \cite{ap-m2as}, the situation is quite different for immiscible compressible 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 \cite{CG-MS1,CG-MS2,CG-MS3}. Namely, in \cite{CG-MS1} certain terms related to the compressibility are neglected, and in \cite{CG-MS2,CG-MS3} the mass densities are assumed not to depend on the physical pressure, but on Chavent's global pressure. In the articles \cite{caro,CG-MS4,khal-saad1,khal-saad2}, a more general immiscible compressible two-phase flow model in porous media is considered for fields with a single rock type and \cite{app-1} treated the case with several types of rocks. In \cite{Siam,app-2} 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 \cite{AJ,AJZ1} 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 \cite{AJZ2} and homogenization results in \cite{AJV}. 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 \cite{ap-m2as} 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 \cite{ap-m2as} is no longer valid. To obtain these results we elaborated a new approach based on the ideas from \cite{Brezis,Gallouet} 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 transport 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 $\varepsilon^2$ surrounded by a system of connected fissures, $\varepsilon$, is a small parameter which characterizes the periodicity of the blocks. There are two kinds of degeneracy 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~\ref{micromodel} 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 \ref{uni-est} is devoted to the presentation of some \emph{a priori} estimates for the solutions of the problem. They are essentially based on an energy equality. In Section~\ref{ex-comp-sf}, 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., \cite{Brezis,Gallouet}). Finally, we formulate the corresponding two--scale convergence results. In Section \ref{dil-oper} 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~\ref{main-res} and its proof is given in Section~\ref{proof-main}. The proof is based on the two-scale convergence and the dilation techniques. \section{Formulation of the problem} \label{micromodel} The outline of the section is as follows. First, in subsection~\ref{micromodel-simpl} 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 \ref{def-glp-simp}. Then in the last subsection~\ref{main-assump}, we give the definition of a weak solution to our system. \subsection{Microscopic model} \label{micromodel-simpl} We consider a reservoir $\Omega \subset \mathbb{R}^d$ ($d = 2, 3$) which is assumed to be a bounded, connected Lipschitz domain with a periodic microstructure. More precisely, we will scale this periodic structure by a parameter $\varepsilon$ which represents the ratio of the cell size to the whole region $\Omega$ and we assume that $0 < \varepsilon \ll 1$ is a small parameter tending to zero. Let $Y := (0,1)^d$ be a periodicity cell; we pave $\mathbb{R}^d$ with $Y$. We assume that $Y_m$ is an open set with piecewise smooth boundary $\partial Y_m$ such that $Y_m \Subset Y$ and we reproduce $Y_m$ by periodicity, obtaining a periodic open set ${\mathsf M}$ in $\mathbb{R}^d$. We denote by ${\mathsf F}$ the periodic set ${\mathsf F} := \mathbb{R}^d \setminus {\mathsf M}$, which is obtained from the set $Y_{f} := Y \setminus Y_m$. Thus $Y = Y_m \cup Y_{f} \cup \Gamma_{fm}$, where $\Gamma_{fm}$ denotes the interface between the two media. Finally, we denote by $\chi_{f}$ and $\chi_m$ the characteristic functions of the sets ${\mathsf F}$ and ${\mathsf M}$. Then $\chi_m(\frac{x}{\varepsilon})$ is the periodic function of period $\varepsilon Y$ which takes the value $1$ in the set ${\mathsf M}^\varepsilon$, union of the sets obtained from $\varepsilon Y_m$ by translations of vectors $\varepsilon \sum_{i=1}^n k_i \vec e_i$, where $k_i \in \mathbb{Z}$ and $\vec e_i$, $1 \leq i \leq d$, is the canonical basis of $\mathbb{R}^d$, and which takes the value $0$ in the set ${\mathsf F}^\varepsilon$, complementary in $\mathbb{R}^d$ of this union. In other words, $\chi_m(\frac{x}{\varepsilon})$ is the characteristic function of the set ${\mathsf M}^\varepsilon$, while $\chi_{f}(\frac{x}{\varepsilon})$ is the characteristic function of ${\mathsf F}^\varepsilon$. Now we can define the subdomains $\Omega^\varepsilon_r$ with $r = ``f''$ or $``m''$ corresponding to the porous medium with the index $``r''$. We set: $$ \Omega_m^\varepsilon := \{x \in \Omega: \chi_m^\varepsilon (x) = 1\} \quad \text{and} \quad \Omega_{f}^\varepsilon := \left\{ x\in\Omega: \chi_{f}^\varepsilon (x) = 1\right\}. $$ Then $\Omega = \Omega^\varepsilon_m \cup \Gamma^\varepsilon_{fm} \cup \Omega^\varepsilon_{f}$, where $\Gamma^\varepsilon_{fm} := \partial \Omega^\varepsilon_{f} \cap \partial \Omega^\varepsilon_m \cap \Omega$ and the subscript $m$ and $f$ refer to the matrix and fracture, respectively. For the sake of simplicity, we assume that $\Omega^\varepsilon_m \cap \partial \Omega = \emptyset$. We also set: \begin{equation} \label{omeg12++} \Omega_{T} := \Omega \times (0, T), \quad \Omega^\varepsilon_{r,T} := \Omega^\varepsilon_r \times (0,T), \quad \text{and} \quad \Gamma_{f,m}^{\varepsilon,T} := \Gamma_{f,m}^{\varepsilon} \times (0, T), \end{equation} where $T > 0$ is fixed. We consider an immiscible compressible two-phase flow system in a porous medium which fills the domain $\Omega$. We focus here on the general case where both phases are compressible, the phases being $\ell$ and $g$. Let $\Phi^{\varepsilon}(x)$ be the porosity of the reservoir $\Omega$; $K^{\varepsilon}(x)$ be the absolute permeability tensor of $\Omega$; $S_{\ell}^{\varepsilon} = S_{\ell}^{\varepsilon}(x, t)$, $S_g^{\varepsilon} = S_g^{\varepsilon}(x, t)$ be the phase saturations; $k_{r,\ell} = k_{r,\ell}(S_{\ell}^{\varepsilon})$, $k_{r,g} = k_{r,g}(S_g^{\varepsilon})$ be the relative permeabilities of the phases; $p_{\ell}^{\varepsilon} = p_{\ell}^{\varepsilon}(x,t)$, $p_g^{\varepsilon} = p_g^{\varepsilon}(x,t)$ be the phase pressures; $\rho_{\ell}$, $\rho_g$ be the phase densities and $P_{c}$ the capillary pressure. In what follows, for the sake of presentation simplicity we neglect the source terms, and we denote $S^{\varepsilon} = S_{\ell}^{\varepsilon}$. The model for the two-phase flow is described by (see, e.g., \cite{GC-JJ,Chen-Huan-Ma,HR}): \begin{equation} \label{eq:1.1.15} \begin{gathered} 0 \leqslant S^\varepsilon \leqslant 1 \quad \text{in } \Omega_{T};\\ \Phi^\varepsilon(x) \frac{\partial\, \Xi^\varepsilon_{\ell}}{\partial t} - \operatorname{div} \Big\{K^\varepsilon(x) \lambda_{\ell} (S^\varepsilon) \rho_{\ell}(p^\varepsilon_{\ell}) \big(\nabla p^\varepsilon_{\ell} - \rho_{\ell}(p^\varepsilon_{\ell}) \vec{g}\big)\Big\} = 0 \quad \text{in } \Omega_{T}; \\ \Phi^\varepsilon(x) \frac{\partial\, \Xi^\varepsilon_g}{\partial t} - \operatorname{div} \Big\{K^\varepsilon(x) \lambda_g (S^\varepsilon) \rho_g(p^\varepsilon_g)\big(\nabla p^\varepsilon_g - \rho_g(p^\varepsilon_g) \vec{g} \big)\Big\} = 0 \quad \text{in } \Omega_{T}; \\ P_{c}(S^\varepsilon) = p^\varepsilon_g - p^\varepsilon_{\ell} \quad \text{in } \Omega_{T}, \end{gathered} \end{equation} where $\lambda_g(S^{\varepsilon}) = \widetilde\lambda_g(1-S^{\varepsilon})$; $\Xi^\varepsilon_{\ell} := S^{\varepsilon}\rho_{\ell}(p_{\ell}^{\varepsilon})$ and $\Xi^\varepsilon_g := (1 - S^{\varepsilon})\rho_g(p_g^{\varepsilon})$; each function $\gamma^{\varepsilon} := S^\varepsilon, p_{\ell}$, $p_g$, $\Xi^\varepsilon_{\ell}$, and $\Xi^\varepsilon_g$ is defined as: \begin{equation} \label{gggg} \gamma^{\varepsilon}(x,t)=\chi_{f}^{\varepsilon}(x)\gamma_{f}^{\varepsilon}(x,t)+ \chi_m^{\varepsilon}(x)\gamma_m^{\varepsilon}(x,t). \end{equation} The velocities of the phases $\vec{q}_{\ell}^{\varepsilon}$, $\vec{q}_g^{\varepsilon}$ are defined by Darcy--Muskat's law: \begin{gather} \vec{q}_{\ell}^{\varepsilon} := - K^{\varepsilon}(x) \lambda_{\ell}(S_{\ell}^{\varepsilon}) \Big(\nabla p_{\ell}^{\varepsilon}-\rho_{\ell}(p_{\ell}^{\varepsilon}) \vec{g}\Big) \quad \text{with } \lambda_{\ell}(S_{\ell}^{\varepsilon}) := \frac{k_{r,\ell}}{\mu_{\ell}}(S_{\ell}^{\varepsilon}); \label{eq:1.7} \\ \overrightarrow{q_g^{\varepsilon}} := - K^{\varepsilon}(x) \widetilde\lambda_g (S_g^{\varepsilon})\Big(\nabla p_g^{\varepsilon} - \rho_g(p_g^{\varepsilon}) \vec{g}\Big) \quad \text{with } \widetilde\lambda_g(S_g^{\varepsilon}) := \frac{k_{r,g}}{\mu_g}(S_g^{\varepsilon}) \label{eq:1.1.7} \end{gather} with $\vec{g}$, $\mu_{\ell}$, $\mu_g$ being the gravity vector and the viscosities, respectively. Now we specify the boundary and initial conditions. We suppose that the boundary $\partial\Omega$ consists of two parts $\Gamma_1$ and $\Gamma_{2}$ such that $\Gamma_1\cap\Gamma_{2} = \emptyset$, $\partial\Omega = \overline{\Gamma_1}\cup\overline{\Gamma_{2}}$. The boundary conditions are given by \begin{gather} p_g^{\varepsilon}(x, t) = 0 = p_{\ell}^{\varepsilon}(x, t) \quad \text{on } \Gamma_1 \times (0,T); \label{eq:1.1.22} \\ \vec{q}_{\ell}^{\varepsilon} \cdot \overrightarrow{\nu} = \overrightarrow{q_g^{\varepsilon}} \cdot \overrightarrow{\nu} = 0 \quad\text{on}\,\,\, \Gamma_{2}\times(0,T). \label{eq:1.1.24} \end{gather} Finally, the initial conditions read \begin{equation} S^{\varepsilon}(x, 0) = S^{0}(x) \quad \text{and} \quad p_g^{\varepsilon}(x, 0) = p_g^{0}(x) \quad\text{in}\,\, \Omega. \label{eq:1.1.25} \end{equation} \subsection{A fractional flow formulation} \label{def-glp-simp} In the sequel, we use a formulation obtained after transformation using the concept of the global pressure introduced in \cite{ant-kaz-mon1990,GC-JJ}. The global pressure is defined as follows: \begin{equation} P^{\varepsilon} = p_{\ell}^{\varepsilon}-G_{\ell}(S^{\varepsilon})= p_g^{\varepsilon}-G_g(S^{\varepsilon}), \label{eq:1.1.28} \end{equation} where the functions $G_{\ell}(S^{\varepsilon})$, $G_g(S^{\varepsilon})$ are given by: \begin{gather} G_g(S^{\varepsilon}) := G_g(0) + \int_{0}^{S^{\varepsilon}} \frac{\lambda_{\ell}(s)}{\lambda(s)}P_{c}'(s)ds, \label{eq:1.1.31} \\ G_{\ell}(S^{\varepsilon})=G_g(S^{\varepsilon})-P_{c}(S^{\varepsilon}), \label{eq:1.1.32} \end{gather} with $\lambda(s) := \lambda_{\ell}(s) + \lambda_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: \begin{gather} \lambda_{\ell}(S^{\varepsilon})\nabla p_{\ell}^{\varepsilon} + \lambda_g(S^{\varepsilon})\nabla p_g^{\varepsilon} = \lambda(S^{\varepsilon})\nabla P^{\varepsilon}\,, \label{eq:1.1.27} \\ \nabla G_{\ell}(S^{\varepsilon}) = -\frac{\lambda_g(S^{\varepsilon})}{\lambda(S^{\varepsilon})}P_{c}' (S^{\varepsilon})\nabla S^{\varepsilon}. \label{eq:1.1.33} \end{gather} Notice that from \eqref{eq:1.1.31}, \eqref{eq:1.1.33} we obtain \begin{equation} \lambda_{\ell}(S^{\varepsilon})\nabla G_{\ell}(S^{\varepsilon}) = \nabla \beta(S^{\varepsilon}) \quad \text{and} \quad \lambda_g(S^{\varepsilon})\nabla G_g(S^{\varepsilon}) = - \nabla \beta(S^{\varepsilon})\,, \label{eq:1.20} \end{equation} where \begin{equation} \beta(S^{\varepsilon}) := \int_{0}^{S^{\varepsilon}} \alpha(u)\, du \quad \text{with } \alpha(s) := \frac{\lambda_g(s)\lambda_{\ell}(s)}{\lambda(s)} |P_{c}'(s)|. \label{eq:1.1.38} \end{equation} Furthermore, we have the important relation \begin{equation} \lambda_g(S^{\varepsilon}) |\nabla p_g^{\varepsilon}|^2 + \lambda_{\ell}(S^{\varepsilon})|\nabla p_{\ell}^{\varepsilon}|^2 = \lambda(S^{\varepsilon}) |\nabla P^{\varepsilon}|^2 + |\nabla b(S^{\varepsilon})|^2, \label{eq:1.22} \end{equation} where \begin{equation} b(s) := \int_{0}^{s} a(\xi)d\xi \quad \text{with } a(s) := \sqrt{\frac{\lambda_g(s)\lambda_{\ell}(s)}{\lambda(s)}} |P_{c}'(s)|. \label{eq:1.23} \end{equation} If we use the global pressure and the saturation as new unknown functions, then problem \eqref{eq:1.1.15} reads \begin{equation} \label{eq:1.1.34} \begin{gathered} 0 \leqslant S^\varepsilon \leqslant 1 \quad \text{in } \Omega_{T}; \\ \Phi^{\varepsilon}(x) \frac{\partial\Theta_{\ell}^{\varepsilon}}{\partial t} - \operatorname{div} \Big\{ K^{\varepsilon}(x)\widetilde \rho_{\ell}^{\varepsilon} \big[\lambda_{\ell}(S^{\varepsilon})\nabla P^{\varepsilon} + \nabla \beta(S^{\varepsilon}) - \lambda_{\ell}(S^{\varepsilon}) \widetilde\rho_{\ell}^{\varepsilon}\vec{g} \big] \Big\} = 0 \quad \text{in } \Omega_{T}; \\ \Phi^{\varepsilon}(x)\frac{\partial\Theta_g^{\varepsilon}}{\partial t} -\operatorname{div} \Big\{ K^{\varepsilon}(x)\widetilde\rho_g^{\varepsilon} \big[\lambda_g(S^{\varepsilon})\nabla P^{\varepsilon} - \nabla\beta(S^{\varepsilon}) - \lambda_g(S^{\varepsilon})\widetilde{\rho_g}^{\varepsilon}\vec{g} \big] \Big\} = 0 \quad \text{in } \Omega_{T}, \end{gathered} \end{equation} we introduced the notation \begin{gather} \widetilde\rho_{\ell}^{\varepsilon} := \rho_{\ell}(P^{\varepsilon} + G_{\ell}(S^{\varepsilon})) \quad \text{and} \quad \widetilde\rho_g^{\varepsilon} := \rho_g(P^{\varepsilon}+G_g(S^{\varepsilon})); \label{eq:1.1.37} \\ \Theta_{\ell}^{\varepsilon} = \Theta_{\ell}(S^{\varepsilon},P^{\epsilon}) := S^{\varepsilon}\widetilde\rho_{\ell}^{\varepsilon} \quad \text{and} \quad \Theta_g^{\varepsilon} = \Theta_g(S^{\varepsilon},P^{\varepsilon}) := (1-S^{\varepsilon}) \widetilde\rho_g^{\varepsilon}. \label{eq:1.27} \end{gather} The system \eqref{eq:1.1.34} is completed by the following boundary and initial conditions. \begin{gather} P^{\varepsilon} = 0 \quad \text{on } \Gamma_1\times(0,T);\label{eq:1.1.39} \\ Q_{\ell}^{\varepsilon} \cdot \overrightarrow{\nu} = Q_g^{\varepsilon} \cdot \overrightarrow{\nu} \quad \text{on } \Gamma_{2} \times (0,T) \label{eq:1.1.40} \end{gather} where $Q_{\ell}^{\varepsilon}$ and $Q_g^{\varepsilon}$ are defined by: \begin{gather*} Q_{\ell}^{\varepsilon} := -K^{\varepsilon}(x)\widetilde\rho_{\ell}^{\varepsilon} \big[ \lambda_{\ell}(S^{\varepsilon})\nabla P^{\varepsilon} + \nabla\beta(S^{\varepsilon}) - \lambda_{\ell}(S^{\varepsilon}) \widetilde\rho_{\ell}^{\varepsilon}\vec{g} \big], \\ Q_g^{\varepsilon} := -K^{\varepsilon}(x) \widetilde\rho_g^{\varepsilon} \big[\lambda_g(S^{\varepsilon})\nabla P^{\varepsilon} - \nabla\beta(S^{\varepsilon}) - \lambda_g(S^{\varepsilon}) \widetilde\rho_g^{\varepsilon} \vec{g}\big]. \end{gather*} Finally, the initial conditions read \begin{equation} S^{\varepsilon}(x,0) = S^{0}(x) \quad \text{and} \quad P^{\varepsilon}(x,0) = P^{0}(x) \quad \text{in} \quad \Omega. \label{eq:0.1.33} \end{equation} \subsection{A weak formulation of the problem} \label{main-assump} Let us begin this section by stating the following assumptions. \begin{itemize} \item[(A1)] The porosity $\Phi=\Phi(y)$ is a $Y$-periodic function defined by: $\Phi^{\varepsilon}(x)=\chi_{f}^{\varepsilon}(x)\Phi_{f} + \chi_m^{\varepsilon}(x)\Phi_m$ with $0 < \Phi_{f}, \Phi_m < 1$, where $\Phi_{f}$ and $\Phi_m$ are constant that do not depend on $\varepsilon$. \item[(A2)] The absolute permeability tensor $K^\varepsilon$ is given by: \[ K^{\varepsilon}(x) := K_{f} \chi_{f}^{\varepsilon}(x)\,\mathbb{I} + \varepsilon^2 K_m \chi_m^{\varepsilon}(x)\mathbb{I}, \] where $\mathbb{I}$ is the unit tensor and $K_{f}$, $K_m$ are positive constants that do not depend on $\varepsilon$. \item[(A3)] The density $\rho_{\rm k} = \rho_{\rm k}(p)$, (${\rm k} = \ell,g$) is a monotone $C^{1}$-function in $\mathbb{R}$ such that \begin{equation} \begin{gathered} \rho_{\rm k}(p) = \rho_{\rm min} \quad \text{for } p \leqslant p_{\rm min};\quad \rho_{\rm k}(p)=\rho_{\rm max} \quad \text{for } p \geqslant p_{\rm max}; \\ \rho_{\rm min} < \rho_{\rm k}(p) < \rho_{\rm max} \quad \text{for } p_{\rm min}

0$ in $[0, 1]$. \end{itemize} Moreover, $\lambda_\ell(s) \sim s^{\varkappa_\ell}$ as $s \to 0$ and $\lambda_g(s) \sim (1 - s)^{\varkappa_g}$ as $s \to 1$ $(\varkappa_\ell, \varkappa_g > 0)$. \item[(A6)] The function $\alpha$ given by \eqref{eq:1.1.38} is a continuous function in $[0,1].$ Moreover, $\alpha(0) = \alpha(1) = 0$ and $\alpha>0$ in $(0,1)$. \item[(A7)] The function $\beta^{-1}$, inverse of $\beta$ defined in \eqref{eq:1.1.38} is a H\"older function of order $\theta$ with $\theta\in(0,1)$ on the interval $[0,\beta(1)].$ Namely, there exists a positive constant $C_{\beta}$ such that for all $s_1,s_{2}\in[0,\beta(1)]$, $|\beta^{-1}(s_1)-\beta^{-1}(s_{2})| \leqslant C_{\beta}|s_1-s_{2}|^{\theta}$. \item[(A8)] The initial data for the global pressure and the saturation defined in \eqref{eq:0.1.33} are such that ${P}^{\bf 0} \in L^2(\Omega)$ and $0 \leqslant {S}^{\bf 0} \leqslant 1$. \end{itemize} 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 \cite{app-1} that dealt with the existence of a weak solution of the studied problem. Next we introduce the Sobolev space $$ H_{\Gamma_1}^{1} := \{ u\in H^{1}(\Omega): u=0 \text{ on } \Gamma_1\}, $$ wich is a Hilbert space when it is equipped with the norm \[ \| u\|_{H_{\Gamma_1(\Omega)}^{1}} = \|\nabla u\|_{(L^2(\Omega))^{d}}\,. \] \begin{definition} \label{def2.1} \rm We say that the pair of functions $\langle P^{\varepsilon},S^{\varepsilon}\rangle$ is a weak solution to problem \eqref{eq:1.1.34}--\eqref{eq:0.1.33} if \begin{itemize} \item[(i)] $0\leqslant S^{\varepsilon} \leqslant 1$ a.e in $\Omega_{T}$. \item[(ii)] $P^{\varepsilon} \in L^2(0,T;H_{\Gamma_1}^{1}(\Omega)).$ \item[(iii)] The boundary conditions \eqref{eq:1.1.39}--\eqref{eq:1.1.40} are satisfied. \item[(iv)] For any $\varphi_{\ell},\varphi_g\in C^{1}([0,T];H_{\Gamma_1}^{1}(\Omega))$ satisfying $\varphi_{\ell}(T) = \varphi_g(T) = 0$, we have \begin{gather} \begin{aligned} &- \int_{\Omega_{T}}\Phi^{\varepsilon}(x)\Theta_{\ell}^{\varepsilon} \frac{\partial\varphi_{\ell}^{\varepsilon}}{\partial t}\,\,dx\,dt + \int_{\Omega} \Phi^{\varepsilon}(x)\Theta_{\ell}^{0}\varphi_{\ell}^{0}\,dx \\ & + \int_{\Omega_{T}} K^{\varepsilon}(x) \widetilde\rho_{\ell}^{\varepsilon} \Big\{ \lambda_{\ell}(S^{\varepsilon})\big(\nabla P^{\varepsilon} - \widetilde\rho_{\ell}^{\varepsilon}\vec{g}\big) + \nabla\beta(S^{\varepsilon})\Big\} \cdot\nabla\varphi_{\ell}^{\varepsilon}\,\,dx\,dt = 0; \end{aligned} \label{eq:1.2.2} \\ \begin{aligned} &- \int_{\Omega_{T}} \Phi^{\varepsilon}(x) \Theta_g^{\varepsilon} \frac{\partial\varphi_g^{\varepsilon}}{\partial t}\,\,dx\,dt + \int_{\Omega} \Phi^{\varepsilon}(x) \Theta_g^{0} \varphi_g^{0}\, dx \\ & + \int_{\Omega_{T}} K^{\varepsilon}(x) \widetilde\rho_g^{\varepsilon} \Big\{ \lambda_g(S^{\varepsilon}) \big(\nabla P^{\varepsilon}-\widetilde\rho_g^{\varepsilon} \vec{g}\big) - \nabla\beta(S^{\varepsilon})\Big\} \cdot\nabla\varphi_{\ell}^{\varepsilon}\,\,dx\,dt = 0, \end{aligned}\label{eq:1.2.3} \end{gather} where $\widetilde\rho_g^{\varepsilon}$ and $\widetilde\rho_{\ell}^{\varepsilon}$ are defined in \eqref{eq:1.1.37}; $\varphi_{\ell}^{0} = \varphi_{\ell}(0,x)$, $\varphi_g^{0} = \varphi_g(0,x)$; $\Theta_{\ell}^{0} = S^{0}\rho_{\ell}(P^{0} + G_{\ell}(S^{0}))$ and $\Theta_g^{0} = (1-S^{0}) \rho_g (P^{0} + G_g(S^{0}))$. \end{itemize} \end{definition} According to \cite{app-1}, under conditions {(A1)--(A8)}, for each $\varepsilon>0$, problem \eqref{eq:1.2.2}--\eqref{eq:1.2.3} has at least one weak solution. In what follows $C, C_1,\dots$ denote generic constants that do not depend on $\varepsilon$. \section{A priori estimates} \label{uni-est} To obtain the needed uniform estimates for the solution of problem \eqref{eq:1.1.15} (or the equivalent problem \eqref{eq:1.1.34}), we follow the choice of the test functions as in \cite{ap-m2as}: \begin{equation} R_{\ell}(p_{\ell}^{\varepsilon}) = \int_{0}^{p_{\ell}^{\varepsilon}} \frac{d\xi}{\rho_{\ell}(\xi)} \quad \text{and} \quad R_g(p_g^{\varepsilon}) = \int_{0}^{p_g^{\varepsilon}}\frac{d\xi}{\rho_g(\xi)}. \label{eq:1.3.1} \end{equation} Then, as in \cite{ap-m2as}, the following results hold. \begin{lemma} Let $\langle p_g^{\varepsilon},p_{\ell}^{\varepsilon}\rangle$ be a solution to \eqref{eq:1.1.15}. Then we have the energy equality \begin{equation} \begin{aligned} &\frac{d}{dt} \int_{\Omega}\Phi^{\varepsilon}(x)\zeta^{\varepsilon}(x, t)\, dx + \int_{\Omega} K^{\varepsilon}(x) \Big\{ \lambda_{\ell}(S^{\varepsilon}) \nabla p_{\ell}^{\varepsilon} \cdot \big[\nabla p_{\ell}^{\varepsilon} - \rho_{\ell}(p_{\ell}^{\varepsilon})\vec{g}\big] \\ &+ \lambda_g(S^{\varepsilon}) \nabla p_g^{\varepsilon} \cdot \big[\nabla p_g^{\varepsilon} - \rho_g(p_g^{\varepsilon}) \vec{g} \big] \Big\}\, dx = 0 \end{aligned} \label{eq:1.3.2} \end{equation} in the sense of distributions. Here \[ \zeta^{\varepsilon} := S^{\varepsilon}\mathcal{R}_{\ell} (p_{\ell}^{\varepsilon}) + (1-S^{\varepsilon})\mathcal{R}_g (p_g^{\varepsilon}) + F(S^{\varepsilon}), \] where $\mathcal{R}_{\rm k}(p) := \rho_{\rm k}(p) R_{\rm k}(p) - p$, $({\rm k} = \ell, g)$ and $F(s) := \int_1^{s} P_{c}(u)\, du$. Moreover, $\zeta^{\varepsilon} \geqslant 0$ in $\Omega_{T}$. \end{lemma} \begin{lemma} \label{lem7} Let $\langle p_g^{\varepsilon},p_{\ell}^{\varepsilon}\rangle$ be a solution to \eqref{eq:1.1.15}. Then \begin{equation} \| \sqrt{K^{\varepsilon}(x)\lambda_{\ell}(S^{\varepsilon})}\nabla p_{\ell}^{\varepsilon}\|_{L^2(\Omega_{T})} + \| \sqrt{K^{\varepsilon}(x) \lambda_g(S^{\varepsilon})} \nabla p_g^{\varepsilon} \|_{L^2(\Omega_{T})}\leqslant C. \label{eq:1.3.12} \end{equation} \end{lemma} \begin{corollary} \label{cora-1} Let $\langle p_g^{\varepsilon},p_{\ell}^{\varepsilon}\rangle$ be a solution to \eqref{eq:1.1.15}. Then \begin{equation} \label{hhh8-cor1} \begin{aligned} &\| \sqrt{\lambda_{\ell} (S^\varepsilon_{f})} \nabla p^\varepsilon_{\ell,f}\|_{L^2(\Omega^\varepsilon_{f,T})} +\| \sqrt{\lambda_g (S^\varepsilon_{f})} \nabla p^\varepsilon_{g,f} \|_{L^2(\Omega^\varepsilon_{f,T})}\\ &+ \varepsilon \| \sqrt{\lambda_{\ell} (S^\varepsilon_m)} \nabla p^\varepsilon_{\ell,m} \|_{L^2(\Omega^\varepsilon_{m,T})} + \varepsilon \| \sqrt{\lambda_g (S^\varepsilon_m)} \nabla p^\varepsilon_{g,m} \|_{L^2(\Omega^\varepsilon_{m,T})} \leqslant C. \end{aligned} \end{equation} \end{corollary} Then we obtain the following uniform \emph{a priori} estimates for the functions $P^{\varepsilon}$ and $\beta(S^{\varepsilon})$. \begin{lemma} \label{lem8} Let the pair of functions $\langle P^{\varepsilon},S^{\varepsilon}\rangle$ be a solution to \eqref{eq:1.1.34}. Then \begin{equation}\label{beta-un} \begin{aligned} &\| \nabla \beta(S^\varepsilon_{f})\|_{L^2(\Omega^\varepsilon_{f,T})} +\| \nabla P^\varepsilon_{f} \|_{L^2(\Omega^\varepsilon_{f,T})} +\varepsilon\| \nabla \beta(S^\varepsilon_m)\|_{L^2(\Omega^\varepsilon_{m,T})}\\ &+\varepsilon \| \nabla P^\varepsilon_m \|_{L^2(\Omega^\varepsilon_{m,T})} \leqslant C. \end{aligned} \end{equation} Moreover, \begin{gather} \| P_{f}^{\varepsilon}\|_{L^2(\Omega^\varepsilon_{f,T})} + \| \beta(S^{\varepsilon})\|_{L^2(\Omega_{T})} \leqslant C, \label{eq:1.3.25} \\ \| P_m^{\varepsilon}\|_{L^2(\Omega^\varepsilon_{m,T})} \leqslant C. \label{eq:1.3.26} \end{gather} \end{lemma} Now we pass to the uniform bounds for the time derivatives of the functions $\Theta^\varepsilon_g$ and $S^\varepsilon$. In a standard way (see, e.g., \cite{Siam,ap-m2as}) we can prove the following lemma. \begin{lemma} \label{cor-ps-l5} Let the pair of functions $\langle P^{\varepsilon},S^{\varepsilon}\rangle$ be a solution to \eqref{eq:1.1.34}. Then for $r = f, m$, \begin{gather} \label{theta-tt-1-cor} \{\partial_t(\Phi_r \Theta^\varepsilon_{\ell,r}) \}_{\varepsilon>0} \text{ is uniformly bounded in } L^2(0,T;H^{-1}(\Omega^\varepsilon_r)); \\ \label{theta-tt-2-cor} \{\partial_t(\Phi_r \Theta^\varepsilon_{g,r}) \}_{\varepsilon>0} \text{ is uniformly bounded in } L^2(0,T;H^{-1}(\Omega^\varepsilon_r)). \end{gather} \end{lemma} \section{Convergence of $\{P^\varepsilon_{f}\}_{\varepsilon>0}$, $\{S^\varepsilon_{f}\}_{\varepsilon>0}$, $\{\Theta^\varepsilon_{\ell,f}\}_{\varepsilon>0}$, $\{\Theta^\varepsilon_{g,f}\}_{\varepsilon>0}$} \label{ex-comp-sf} In this section, we obtain compactness results that will be used in passing to the limit as $\varepsilon$ tends to zero in the weak formulation. The compactness result used in \cite{ap-m2as} is no longer valid. To obtain these results we elaborated a new approach based on the ideas from \cite{Brezis,Gallouet} 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~\ref{ss-sf-1} we extend the function $S^\varepsilon_{f}$ from the subdomain $\Omega^\varepsilon_{f}$ to the whole $\Omega$ and obtain uniform estimates for the extended function $\widetilde S^\varepsilon_{f}$. Then in Section \ref{ss-sf-2}, using the uniform estimates for the function $\widetilde P^\varepsilon_{f}$ which follow from Lemma~\ref{lem8}, the definition of the extension operator, and the corresponding bounds for $\widetilde S^\varepsilon_{f}$, we prove the compactness result for the families $\{\widetilde \Theta^\varepsilon_{\ell,f}\}_{\varepsilon>0}$ and $\{\widetilde \Theta^\varepsilon_{g,f}\}_{\varepsilon>0}$, where $\widetilde\Theta^\varepsilon_{\ell,f}$, $\widetilde\Theta^\varepsilon_{g,f}$ are extensions of the functions $\Theta^\varepsilon_{\ell,f}$, $\Theta^\varepsilon_{g,f}$ from the subdomain $\Omega^\varepsilon_{f}$ to the whole $\Omega$ which will be specified at the end of subsection~\ref{ss-sf-1}. Finally, in subsection~\ref{con-results} we formulate the two--scale convergence which will be used in the derivation of the homogenized system. \subsection{Extensions of the functions $S^\varepsilon_{f}$, $\Theta^\varepsilon_{\ell,f}$, $\Theta^\varepsilon_{g,f}$} \label{ss-sf-1} Consider the function $S^\varepsilon_{f}$. To extend $S^\varepsilon_{f}$, following the ideas of \cite{blm}, we use the monotone function $\beta$ defined in \eqref{eq:1.1.38}. Let us introduce the function \begin{equation} \label{ue-23} \beta^\varepsilon_{f}(x, t) := \beta(S^\varepsilon_{f}) = \int_0^{S^\varepsilon_{f}} \alpha(u)\, du. \end{equation} Then it follows from condition (A6) that $0 \leqslant \beta^\varepsilon_{f} \leqslant \max_{s\in[0, 1]}\alpha(s)$ a.e. in $\Omega^\varepsilon_{{f},T}$. Furthermore, from \eqref{beta-un} we have \begin{equation} \label{ue-25} \| \nabla \beta^\varepsilon_{f} \|_{L^2(\Omega^\varepsilon_{{f},T})} \leqslant C. \end{equation} Let $\Pi^\varepsilon : H^1(\Omega^\varepsilon_{f}) \to H^1(\Omega)$ be the standard extension operator cf. \cite{acdp}. Then we have \begin{gather*} 0 \leqslant \widetilde\beta^\varepsilon_{f} := \Pi^\varepsilon \beta^\varepsilon_{f} \leqslant \max_{s\in[0, 1]} \alpha(s) \quad \text{a.e. in } \Omega_T\,, \| \nabla \widetilde \beta^\varepsilon_{f} \|_{L^2(\Omega_{T})} \leqslant C. \end{gather*} Now we can extend the function $S^\varepsilon_{f}$ from the subdomain $\Omega^\varepsilon_{f}$ to the whole $\Omega$. We denote this extension by $\widetilde S^\varepsilon_{f}$ and define it as: $$ \widetilde S^\varepsilon_{f} := (\beta)^{-1}(\widetilde \beta^\varepsilon_{f}). $$ This implies that \begin{gather} \label{ue-31} \int_{\Omega_T} |\nabla \beta(\,\widetilde S^\varepsilon_{f}\,)|^2 \,dx\, dt \leqslant C\,, 0 \leqslant \widetilde S^\varepsilon_{f} \leqslant 1 \quad \text{a.e. in } \Omega_T. \end{gather} Finally consider the sequences $\{\Theta^\varepsilon_{\ell,f}\}_{\varepsilon>0}$ and $\{\Theta^\varepsilon_{g,f}\}_{\varepsilon>0}$. We recall that $\Theta^\varepsilon_{\ell,f} := \rho_{\ell}\,\big(P^\varepsilon_{f} + G_{\ell}(S^\varepsilon_{f})\big) S^\varepsilon_{f}$ and $\Theta^\varepsilon_{g,f} := \rho_g \big(P^\varepsilon_{f} + G_g(S^\varepsilon_{f})\big) (1 - S^\varepsilon_{f})$. Then we define the extension of the function $\Theta^\varepsilon_{f}$ to the whole $\Omega$ by \begin{equation} \widetilde\Theta_{\ell,f}^{\varepsilon} := \rho_{\ell} (\widetilde P_{f}^{\varepsilon} + G_{\ell}(\widetilde S_{f}^{\varepsilon})) \widetilde S_{f}^{\varepsilon} \quad \text{and} \quad \widetilde \Theta_{g,f}^{\varepsilon} := \rho_g(\widetilde P_{f}^{\varepsilon} + G_g(\widetilde S_{f}^{\varepsilon})) (1-\widetilde S_{f}^{\varepsilon}), \label{eq:1.4.2} \end{equation} where $\widetilde P^\varepsilon_{f} := \Pi^\varepsilon P^\varepsilon_{f}$ is the extension of the function $P^\varepsilon_{f}$. \subsection{Compactness of the sequences $\{\widetilde\Theta_{\ell,f}^{\varepsilon}\}_{\varepsilon>0}$, $\{\widetilde\Theta_{g,f}^{\varepsilon}\}_{\varepsilon>0}$} \label{ss-sf-2} The following convergence result is valid. \begin{proposition} \label{prop9} Under our standing assumptions there exist the functions ${\mathsf L_1}, {\mathsf L_2}$ such that \begin{gather}\label{L-1} \widetilde\Theta_{\ell,f}^{\varepsilon} \to {\mathsf L_1} \text{ strongly in } L^2(\Omega_T), \\ \label{L-2} \widetilde\Theta_{g,f}^{\varepsilon} \to {\mathsf L_2} \text{ strongly in } L^2(\Omega_T). \end{gather} \end{proposition} \begin{proof} We will prove the convergence result for the sequence $\{\widetilde\Theta_{\ell,f}^{\varepsilon}\}_{\varepsilon>0}$, the corresponding result for the sequence $\{\widetilde\Theta_{g,f}^{\varepsilon}\}_{\varepsilon>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., \cite{Brezis,Gallouet}) in the space $L^{1}(\Omega_{T})$. To this end we have to obtain the moduli of continuity with respect to the space and temporal variables (see Lemmata \ref{lem10}, \ref{lem11} below). Finally, the uniform boundedness of the function $\Theta_{\ell,f}^{\varepsilon}$ imply the desired convergence result \eqref{L-1} in the space $L^2(\Omega_{T})$. \end{proof} We start with the following result which can be proved by arguments similar to those from Lemma~4.2 in \cite{ap-m2as}. \begin{lemma}[Modulus of continuity with respect to the space variable] \label{lem10} Let conditions {\rm (A1)--(A8)} be fulfilled. Then for $|\Delta x|$ sufficiently small, \begin{equation} \int_{0}^{T} \int_{\Omega} \big| \widetilde\Theta_{\ell,f}^{\varepsilon}(x + \Delta x, t) - \widetilde\Theta_{\ell,f}^{\varepsilon}(x, t) \big|^2\, dx\, dt \leqslant C\, |\Delta x|^{\theta}, \label{eq:1.5.1} \end{equation} where $\theta \in (0, 1)$ is defined in condition {(A7)}. \end{lemma} Now we turn to the derivation of the modulus of continuity with respect to the temporal variable. To do this, for any $\delta>0$, we introduce the functions $$ S^{\varepsilon,\delta}_r := \min\big\{1-\delta, \max(\delta, S^\varepsilon_r)\big\} \quad \text{with } r = f, m. $$ Let us estimate the norm of the function $S^{\varepsilon,\delta}_r$ ($r = f, m$) in $L^2(0,T;H^1(\Omega^\varepsilon_{r,T}))$. \begin{lemma} \label{lem-cdelta} Under our standing assumptions, \begin{equation} \label{nouv-idee-2} \| S^{\varepsilon,\delta}_f \|_{L^2(0,T;H^1(\Omega^\varepsilon_{f,T}))} + \varepsilon \| S^{\varepsilon,\delta}_m \|_{L^2(0,T;H^1(\Omega^\varepsilon_{m,T}))} \leqslant C\,\delta^{-\eta}, \end{equation} where $\eta := (\varkappa_\ell+\varkappa_g)$ and $C$ is a constant that does not depend on ${\varepsilon,\delta}$. \end{lemma} \begin{proof} We consider the function $S^{\varepsilon,\delta}_f$, the estimate for the function $S^{\varepsilon,\delta}_m$ can be obtained in a similar way. It is evident that $S^{\varepsilon,\delta}_f$ (as well as $S^{\varepsilon}_f$) belongs to the space $L^2(\Omega^\varepsilon_{f,T})$. Now we are going to estimate $\nabla S^{\varepsilon,\delta}_f$. From Lemma \ref{lem8} we know that $$ \| \nabla \beta(S^\varepsilon_{f})\|_{L^2(\Omega^\varepsilon_{f,T})} \leqslant C. $$ Then it is clear that \begin{equation} \label{nouv-idee-4} \| \nabla \beta(S^{\varepsilon,\delta}_{f})\|_{L^2(\Omega^\varepsilon_{f,T})} \leqslant \| \nabla \beta(S^\varepsilon_{f})\|_{L^2(\Omega^\varepsilon_{f,T})} \leqslant C, \end{equation} where $C$ is a constant that does not depend on $\varepsilon, \delta$. Moreover, condition (A5) implies the inequalities: \begin{equation}\label{nouv-idee-5} \lambda_\ell(S^{\varepsilon,\delta}_{f}) \geqslant C_1\delta^{\varkappa_\ell} \quad \text{and} \quad \lambda_g(S^{\varepsilon,\delta}_{f}) \geqslant C_2\delta^{\varkappa_g}. \end{equation} Then taking into account the definition of the functions $\beta$ and $\alpha$ (see \eqref{eq:1.1.38}) and conditions (A4), (A5), from \eqref{nouv-idee-5}, we have that \begin{equation} \label{nouv-idee-6} \alpha(S^{\varepsilon,\delta}_{f}) \geqslant \mathbb{C}_1 \delta^{(\varkappa_\ell+\varkappa_g)} \quad {\rm with} \,\, \mathbb{C}_1 := \frac{C_1 C_2}{2} \min_{s\in[0,1]} |P'_c(s)|. \end{equation} Then from \eqref{nouv-idee-4}, \eqref{nouv-idee-6} \begin{equation}\label{nouv-idee-7} \begin{aligned} C &\geqslant \| \nabla \beta(S^{\varepsilon,\delta}_{f})\|^2_{L^2(\Omega^\varepsilon_{f,T})} = \int_0^T \int_{\Omega^\varepsilon_{f}} \alpha^2(S^{\varepsilon,\delta}_{f}) |\nabla S^{\varepsilon,\delta}_{f}|^2\, dx\, dt \\ &\geqslant \mathbb{C}^2_1 \delta^{2(\varkappa_\ell+\varkappa_g)} \int_0^T \int_{\Omega^\varepsilon_{f}} |\nabla S^{\varepsilon,\delta}_{f}|^2\,dx\, dt \\ &= \mathbb{C}^2_1 \delta^{2(\varkappa_\ell+\varkappa_g)} \| \nabla S^{\varepsilon,\delta}_{f}\|^2_{L^2(\Omega^\varepsilon_{f,T})}. \end{aligned} \end{equation} The inequality \eqref{nouv-idee-6} implies that \begin{equation} \label{nouv-idee-7b} \| \nabla S^{\varepsilon,\delta}_{f}\|^2_{L^2(\Omega^\varepsilon_{f,T})} \leqslant C\, \delta^{-2(\varkappa_\ell+\varkappa_g)} \end{equation} and inequality \eqref{nouv-idee-2} is proved. This completes the proof of Lemma \ref{lem-cdelta}. \end{proof} In what follows we use a technical lemma which can be proved using the Fubini theorem. \begin{lemma} \label{lem12} For $h$ sufficiently small, $0 < h < \frac{T}{2}$ and for integrable functions ${\EuScript G}_1(t)$, ${\EuScript G}_{2}(t)$ it holds: $$ \int_{0}^{T} {\EuScript G}_1(t) \Big(\int_{\max(t,h)}^{\min(t+h,T)} {\EuScript G}_{2}(\tau)\, d\tau\Big)\, dt = \int_{h}^{T} {\EuScript G}_{2}(t)\Big(\int_{t-h}^{t} {\EuScript G}_1(\tau)\, d\tau \Big)\, dt. $$ \end{lemma} Now, for $\varepsilon > 0$ and $0 < h < \frac{T}{2}$, let us introduce the function \begin{equation} \begin{gathered} \varphi^{{\varepsilon,\delta}, h}(x, t) := \int_{\max(t, h)}^{\min(t+h,T)} h \partial^{-h} \Theta_{\ell}^{{\varepsilon,\delta}}(x, \tau)\, d\tau, \\ \text{with}\quad \partial^{-h} v (t) := \frac{v(t) - v(t - h)}{h}, \end{gathered} \label{eq:1.5.2} \end{equation} where \begin{equation} \label{nouv-idee-10} \Theta_{\ell}^{{\varepsilon,\delta}}(x, t) := \rho_{\ell}(P^{\varepsilon} + G_{\ell}(S^{{\varepsilon,\delta}}))\, S^{{\varepsilon,\delta}}. \end{equation} The properties of $\varphi^{{\varepsilon,\delta},h}$ are described by the next lemma. \begin{lemma} \label{lem13} Let $\varepsilon>0$, $0 < \delta <1$, and let $h > 0$ be small enough. There exist a constant $C$ which does not depend on ${\varepsilon,\delta}$, and $h$ such that for the sequence of functions defined by \eqref{eq:1.5.2} it holds \begin{gather} \varphi^{{\varepsilon,\delta}, h} \in L^2(0,T; H_{\Gamma_1}^{1}(\Omega)); \label{eq:1.5.3-1} \\ \varphi^{{\varepsilon,\delta},h}(x, T) = 0; \label{eq:1.5.3} \\ \| \varphi^{{\varepsilon,\delta},h} \|_{L^2(\Omega_{T})}\leqslant C\, h; \label{eq:1.5.4} \\ \| \nabla\varphi^{{\varepsilon,\delta},h}\|_{L^2(\Omega^{\varepsilon}_{f,T})} \leqslant C\, h\, \delta^{-\eta}; \label{eq:1.5.5} \\ \varepsilon \| \nabla\varphi^{{\varepsilon,\delta},h}\|_{L^2(\Omega^{\varepsilon}_{m,T})} \leqslant C\, h\, \delta^{-\eta}. \label{eq:1.5.5+1} \end{gather} Here \begin{equation} \label{last-point-eta} \eta := (\varkappa_\ell+\varkappa_g). \end{equation} \end{lemma} \begin{proof} The regularity property \eqref{eq:1.5.3-1} follows immediately from Lemma \ref{lem8} and Lemma \ref{lem-cdelta}. Moreover, taking into account that $S^\varepsilon = 1$ and $P^\varepsilon = {\rm const}$ on $\Gamma_1 \times(0, T)$ (see Section \ref{micromodel}), according to the definition of the function $S^{\varepsilon,\delta}$, we have that $\varphi^{{\varepsilon,\delta},h} = 0$ on $\Gamma_1 \times(0, T)$. Result \eqref{eq:1.5.3} follows directly from the definition of the function $\varphi^{{\varepsilon,\delta},h}$. In fact, $$ \varphi^{{\varepsilon,\delta},h}(x, T) = \int_{\max(T,h)}^{\min(T+h,T)} h \partial^{-h}\Theta_{\ell}^{{\varepsilon,\delta}}(x, \tau)\, d\tau = \int_{T}^{T} h\, \partial^{-h}\Theta_{\ell}^{{\varepsilon,\delta}}(x, \tau)\, d\tau = 0. $$ Bound \eqref{eq:1.5.4} also follows immediately from the definition of $\varphi^{{\varepsilon,\delta},h}$ since $\min(t+h, T) - \max(t, h) \leqslant h$ and the function $\Theta_{\ell}^{{\varepsilon,\delta}}$ is uniformly bounded in $L^\infty(\Omega_T)$. For bound \eqref{eq:1.5.5}, we have \begin{equation} \begin{aligned} &\|\nabla\varphi^{{\varepsilon,\delta},h}\|_{L^2(\Omega^{\varepsilon}_{f,T})}^2\\ &\leqslant \int_{\Omega^{\varepsilon}_{f,T}}\, \Big[\int_{\max(t,h)}^{\min(t+h,T)} \big|\nabla \Theta_{\ell,f}^{{\varepsilon,\delta}}(x, \tau) - \nabla \Theta_{\ell,f}^{{\varepsilon,\delta}}(x,\tau-h)\big|\, d\tau \Big]^2\,dx\,dt := {\EuScript J}^{\varepsilon,\delta}. \end{aligned} \label{eq:1.5.6} \end{equation} Let us estimate the right-hand side of \eqref{eq:1.5.6}. Since $[\min(t+h, T) - \max(t, h)] \leqslant h$, from Cauchy's inequality, for a.e. $(x, t) \in \Omega^{\varepsilon}_{f,T}$ we obtain \begin{align*} &\int_{\max(t,h)}^{\min(t+h,T)} \big|\nabla\Theta_{\ell,f}^{{\varepsilon,\delta}}(x, \tau) - \nabla\Theta_{\ell,f}^{{\varepsilon,\delta}}(x,\tau-h)\big| d\tau \\ &\leqslant {h}^{1/2} \Big[\int_{\max(t,h)}^{\min(t+h,T)} \big|\nabla\Theta_{\ell,f}^{{\varepsilon,\delta}}(x, \tau) - \nabla\Theta_{\ell,f}^{{\varepsilon,\delta}}(x,\tau-h)\big|^2 d\tau \Big]^{1/2}. \end{align*} Therefore, from this inequality we obtain \begin{equation} {\EuScript J}^{\varepsilon,\delta} \leqslant C\, h\, \int^{T}_0 \int_{\Omega^{\varepsilon}_{f}} \Big[ \int_{\max(t,h)}^{\min(t+h,T)}\big|\nabla\Theta_{\ell,f}^{{\varepsilon,\delta}}(x,\tau) - \nabla\Theta_{\ell,f}^{{\varepsilon,\delta}}(x,\tau-h)\big|^2\, d\tau\Big]\,dx\,dt. \label{eq:1.5.7} \end{equation} Now we apply Lemma \ref{lem12} with ${\EuScript G}_1(t) := 1$ and ${\EuScript G}_{2}(t) := |\nabla\Theta_{\ell,f}^{{\varepsilon,\delta}} (x,\tau) - \nabla\Theta_{\ell,f}^{{\varepsilon,\delta}}(x,\tau-h)|^2$ in the right--hand side of \eqref{eq:1.5.7}. We have: \begin{align*} &\int^{T}_0 \Big[ \int_{\max(t, h)}^{\min(t+h,T)}\big|\nabla\Theta_{\ell,f}^{{\varepsilon,\delta}} (x, \tau) -\nabla\Theta_{\ell,f}^{{\varepsilon,\delta}}(x,\tau-h)\big|^2\, d\tau\Big] \,dt \\ &=\int_{h}^{T} \big|\nabla\Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t) - \nabla\Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t-h)\big|^2 \, \Big[\int_{t-h}^{t}1\, d\tau\Big]\, dt \\ &= h \int_{h}^{T} \big|\nabla\Theta_{\ell,f}^{{\varepsilon,\delta}} (x, t) -\nabla\Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t-h)\big|^2 \, dt. \end{align*} Then, by \eqref{eq:1.5.7}, we deduce that \begin{equation} \label{point-7} \begin{aligned} &{\EuScript J}^{\varepsilon,\delta}\\ &\leqslant C\, h^2\, \int_{\Omega^{\varepsilon}_{f}} \int^{T}_h \big|\nabla\Theta_{\ell,f}^{{\varepsilon,\delta}} (x, t) -\nabla\Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t-h)\big|^2 \, dt\, dx \\ &\leqslant C h^2 \Big[\int_{\Omega^{\varepsilon}_{f}} \int^{T}_h \big|\nabla\Theta_{\ell,f}^{{\varepsilon,\delta}} (x, t)\big|^2 \, dt\, dx + \int_{\Omega^{\varepsilon}_{f}} \int^{T}_h \big|\nabla \Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t-h)\big|^2 \, dt\, dx \Big] \\ &\leqslant C\, h^2 \| \nabla\Theta_{\ell,f}^{{\varepsilon,\delta}} \|_{L^2(\Omega^{\varepsilon}_{f,T})}^2. \end{aligned} \end{equation} It remains to estimate the right-hand side of \eqref{point-7}. To this end we rewrite $\nabla\Theta_{\ell,f}^{{\varepsilon,\delta}}$ as follows: \begin{equation} \label{nouv-idee-12} \nabla\Theta_{\ell,f}^{{\varepsilon,\delta}} =\rho_{\ell}\,\big(P^{\varepsilon} + G_{\ell}(S^{{\varepsilon,\delta}})\big)\, \nabla S^{{\varepsilon,\delta}} + \rho_\ell' \, S^{{\varepsilon,\delta}} \nabla P^{\varepsilon} + \rho_\ell' \, S^{{\varepsilon,\delta}} \nabla G_{\ell}(S^{{\varepsilon,\delta}}). \end{equation} Then, from \eqref{nouv-idee-12}, conditions {(A3)}, {(A5)}, and the definition of the function $G_\ell$, we have \begin{equation} \label{nouv-idee-13} \begin{aligned} &\| \nabla\Theta_{\ell,f}^{{\varepsilon,\delta}} \|_{L^2(\Omega^{\varepsilon}_{f,T})}^2 \\ &\leqslant \rho^2_{\rm max} \| \nabla S^{{\varepsilon,\delta}}\|_{L^2(\Omega^{\varepsilon}_{f,T})}^2 + C \Big[\| \nabla P^\varepsilon \|_{L^2(\Omega^{\varepsilon}_{f,T})}^2 + \| P'_c \nabla S^{{\varepsilon,\delta}} \|_{L^2(\Omega^{\varepsilon}_{f,T})}^2 \Big]. \end{aligned} \end{equation} To estimate the right-hand side of \eqref{nouv-idee-13}, we make use of condition (A4), bound \eqref{nouv-idee-2}, and Lemma \ref{lem-cdelta}. Taking into account that $\delta$ is sufficiently small, we obtain \begin{equation} \label{nouv-idee-14} \| \nabla\Theta_{\ell,f}^{{\varepsilon,\delta}} \|_{L^2(\Omega^{\varepsilon}_{f,T})}^2 \leqslant C\, \delta^{-2(\varkappa_\ell+\varkappa_g)}. \end{equation} Now, from \eqref{eq:1.5.6}, \eqref{point-7}, and \eqref{nouv-idee-14}, for $\delta$ sufficiently small, we obtain \begin{equation} \label{nouv-idee-16} \|\nabla\varphi^{{\varepsilon,\delta},h}\|_{L^2(\Omega^{\varepsilon}_{f,T})}^2 \leqslant C\, h^2\, \delta^{-2(\varkappa_\ell+\varkappa_g)} \end{equation} and the bound \eqref{eq:1.5.5} is proved. The proof of the bound \eqref{eq:1.5.5+1} can be done by arguments similar to ones used in the proof of \eqref{eq:1.5.5}. This completes the proof of Lemma \ref{lem13}. \end{proof} Now we are in a position to estimate the modulus of continuity of the function $\widetilde\Theta_{\ell,f}^{\varepsilon}$. \begin{lemma}[Modulus of continuity with respect to time ] \label{lem11} Under our standing \\ assumptions, for all $h \in (0,T)$ and $\delta$ sufficiently small, we have \begin{equation} \label{point-000} \int_{h}^{T} \int_{\Omega}\, \big| \widetilde \Theta_{\ell,f}^{\varepsilon}(x, t) - \widetilde \Theta_{\ell,f}^{\varepsilon}(x, t-h)\big|^2\,dx\,dt \leqslant C h^\sigma \quad \text{with } \sigma := \min\big\{\frac{1}{2},\frac{1}{2\eta}\big\}, \end{equation} where $\eta$ is defined in \eqref{last-point-eta} and $C$ is a constant which does not depend on $\varepsilon$ and $h$. \end{lemma} \begin{proof} Let us insert the function $\varphi^{{\varepsilon,\delta},h} = \varphi^{{\varepsilon,\delta},h}(x, t)$ in equation \eqref{eq:1.2.2}. We have \begin{equation} \begin{aligned} &\int_{\Omega_{T}}\Phi^{\varepsilon}(x)\frac{\partial \Theta_{\ell}^{\varepsilon}}{\partial t} \varphi^{{\varepsilon,\delta},h}\,dx\,dt\\ &+ \int_{\Omega_{T}} K^{\varepsilon}(x) \widetilde \rho_{\ell}^{\varepsilon} \Big\{ \lambda_{\ell}(S^{\varepsilon})\big(\nabla P^{\varepsilon} - \widetilde \rho_{\ell}^{\varepsilon}\vec{g}\big) + \nabla\beta(S^{\varepsilon})\Big\} \cdot \nabla \varphi^{{\varepsilon,\delta},h}\,dx\,dt = 0. \end{aligned}\label{eq:1.2.2-lem} \end{equation} Taking into account the definition of the function $\varphi^{{\varepsilon,\delta},h}$ and condition (A1), for the first term in the left-hand side of \eqref{eq:1.2.2-lem}, we have \begin{align*} I_1^{\varepsilon}(\varphi^{{\varepsilon,\delta},h}) :&= \int_{\Omega_{T}}\Phi^{\varepsilon}(x)\frac{\partial \Theta_{\ell}^{\varepsilon}}{\partial t} \varphi^{{\varepsilon,\delta},h}\,dx\,dt\\ & = \int_{0}^{T} \int_{\Omega_{f}^{\varepsilon}} \Phi_{f} \frac{\partial\Theta_{\ell,f}^{\varepsilon}}{\partial t} \Big[ \int_{\max(t, h)}^{\min(t+h,T)} h\partial^{-h} \Theta_{\ell,f}^{{\varepsilon,\delta}}(x,\tau)\, d\tau\Big]\, dx\, dt \\ &\quad + \int_{0}^{T} \int_{\Omega_m^{\varepsilon}} \Phi_m \frac{\partial\Theta_{\ell,m}^{\varepsilon}}{\partial t} \Big[ \int_{\max(t,h)}^{\min(t+h,T)} h\partial^{-h}\Theta_{\ell,m}^{{\varepsilon,\delta}}(x,\tau) \, d\tau \Big]\,dx\, dt. \end{align*} Now Lemma \ref{lem12} implies \begin{equation} \label{new-goods-1} \begin{aligned} I_1^{\varepsilon}(\varphi^{{\varepsilon,\delta},h}) &= \int_{h}^{T} \int_{\Omega_{f}^{\varepsilon}} \Phi_{f}\, h\, \partial^{-h}\Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t) \Big[\int_{t-h}^{t} \frac{\partial\Theta_{\ell,f}^{\varepsilon}}{\partial\tau}\ , d\tau \Big]\,dx\,dt \\ &\quad +\int_{h}^{T} \int_{\Omega_m^{\varepsilon}} \Phi_m\, h\, \partial^{-h}\Theta_{\ell,m}^{{\varepsilon,\delta}}(x, t) \Big[\int_{t-h}^{t}\frac{\partial\Theta_{\ell,m}^{\varepsilon}}{\partial\tau}\, d\tau \Big]\,dx\,dt \\ &= \Phi_{f} \int_{h}^{T} \int_{\Omega_{f}^{\varepsilon}} \big\{ \Theta_{\ell,f}^{\varepsilon}(x, t) - \Theta_{\ell,f}^{\varepsilon}(x, t-h)\big\} \\ &\quad\times \big\{ \Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t) - \Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t-h)\big\} \,dx\,dt \\ &\quad + \Phi_m \int_{h}^{T} \int_{\Omega_m^{\varepsilon}} \big\{\Theta_{\ell,m}^{\varepsilon}(x, t) -\Theta_{\ell,m}^{\varepsilon}(x, t-h)\big\}\\ &\quad\times \big\{\Theta_{\ell,m}^{{\varepsilon,\delta}}(x, t) -\Theta_{\ell,m}^{{\varepsilon,\delta}}(x, t-h)\big\} \,dx\,dt \\ :&= {\EuScript J}^{\varepsilon,\delta}_f + {\EuScript J}^{\varepsilon,\delta}_m. \end{aligned} \end{equation} Considering the integrand in the first term of \eqref{new-goods-1}, we have \begin{equation} \label{new-goods-2} \begin{aligned} &\big[ \Theta_{\ell,f}^{\varepsilon}(x, t) - \Theta_{\ell,f}^{\varepsilon}(x, t-h)\big] \, \big[ \Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t) - \Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t-h)\big]\\ & =\big| \Theta_{\ell,f}^{\varepsilon}(x, t) - \Theta_{\ell,f}^{\varepsilon}(x, t-h)\big|^2 \\ &\quad + \big[ \Theta_{\ell,f}^{\varepsilon}(x, t) - \Theta_{\ell,f}^{\varepsilon}(x, t-h)\big]\, \big[\Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t) - \Theta_{\ell,f}^{\varepsilon}(x, t) \big] \\ &-\big[ \Theta_{\ell,f}^{\varepsilon}(x, t) - \Theta_{\ell,f}^{\varepsilon}(x, t-h)\big] \, \big[\Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t-h) - \Theta_{\ell,f}^{\varepsilon}(x, t-h) \big] \\ :&=\big| \Theta_{\ell,f}^{\varepsilon}(x, t) - \Theta_{\ell,f}^{\varepsilon}(x, t-h)\big|^2 + \mathbb{T}^{\,{\varepsilon,\delta},h}_1(x, t) - \mathbb{T}^{\,{\varepsilon,\delta},h}_2(x, t). \end{aligned} \end{equation} Then the first term in the right--hand side of \eqref{new-goods-1} can be rewritten as \begin{equation} \label{new-goods-3} \begin{aligned} {\EuScript J}^{\varepsilon,\delta}_f &= \Phi_{f} \int_{h}^{T} \int_{\Omega_{f}^{\varepsilon}} \big| \Theta_{\ell,f}^{\varepsilon}(x, t) - \Theta_{\ell,f}^{\varepsilon}(x, t-h)\big|^2 \,dx\,dt \\ &\quad + \Phi_{f} \int_{h}^{T} \int_{\Omega_{f}^{\varepsilon}} \mathbb{T}^{\,{\varepsilon,\delta},h}_1(x, t) \,dx\,dt -\Phi_{f} \int_{h}^{T} \int_{\Omega_{f}^{\varepsilon}} \mathbb{T}^{\,{\varepsilon,\delta},h}_2(x, t) \,dx\,dt. \end{aligned} \end{equation} It is easy to see that $$ \big|\Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t) - \Theta_{\ell,f}^{\varepsilon}(x, t) \big| \leqslant C\,\delta \quad \text{and} \quad \big|\Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t-h) - \Theta_{\ell,f}^{\varepsilon}(x, t-h) \big| \leqslant C\,\delta, $$ where $C$ is a constant that does not depend on ${\varepsilon,\delta},h$. Let us consider, for example, the first bound, the second one can be obtained by similar arguments. We have \begin{equation} \label{new-goods-30} \begin{aligned} &\big|\Theta_{\ell,f}^{\varepsilon}(x, t) - \Theta_{\ell,f}^{{\varepsilon,\delta}}(x, t) \big| \\ &\leqslant \big|\rho_{\ell}(P^{\varepsilon} + G_{\ell}(S^{\varepsilon}_f)) S^{\varepsilon}_f - \rho_{\ell}(P^{\varepsilon} + G_{\ell}(S^{\varepsilon}_f)) S^{{\varepsilon,\delta}}_f\big|\\ &\quad +\big|\rho_{\ell}(P^{\varepsilon} + G_{\ell}(S^{\varepsilon}_f)) S^{{\varepsilon,\delta}}_f - \rho_{\ell}(P^{\varepsilon} + G_{\ell}(S^{{\varepsilon,\delta}}_f)) S^{{\varepsilon,\delta}}_f\big| \\ &\leqslant \rho_{\rm max}\big|S^{\varepsilon}_f - S^{{\varepsilon,\delta}}_f \big| + \big|\rho_{\ell}(P^{\varepsilon} + G_{\ell}(S^{\varepsilon}_f)) - \rho_{\ell}(P^{\varepsilon} + G_{\ell}(S^{{\varepsilon,\delta}}_f)) \big| \\ &\leqslant \rho_{\rm max}\delta + \rho'_{\ell} \big| G_{\ell}(S^{\varepsilon}_f) - G_{\ell}(S^{{\varepsilon,\delta}}_f) \big| \\ &\leqslant \rho_{\rm max}\delta + \rho'_{\ell}\, \big|P_c(S^{\varepsilon}_f) - P_c(S^{{\varepsilon,\delta}}_f) \big| \\ &\leqslant \rho_{\rm max}\delta + \rho'_{\ell}\,P'_c \big|S^{\varepsilon}_f - S^{{\varepsilon,\delta}}_f \big| \leqslant C\, \delta. \end{aligned} \end{equation} Since the function $\Theta_{\ell,f}^{\varepsilon}$ is uniformly bounded in $\Omega_T$, from \eqref{new-goods-3} and \eqref{new-goods-30} we have \begin{equation} \label{new-goods-4} {\EuScript J}^{\varepsilon,\delta}_f = \Phi_{f} \int_{h}^{T} \int_{\Omega_{f}^{\varepsilon}} \big| \Theta_{\ell,f}^{\varepsilon}(x, t) - \Theta_{\ell,f}^{\varepsilon}(x, t-h)\big|^2 \,dx\,dt + {\mathsf j}^{\varepsilon,\delta}_f \end{equation} \ with $\big|{\mathsf j}^{\varepsilon,\delta}_f \big| \leqslant C\delta$, where $C$ is a constant that does not depend on ${\varepsilon,\delta}$, $h$. In a similar way, one obtains \begin{equation} \label{new-goods-5} {\EuScript J}^{\varepsilon,\delta}_m = \Phi_m \int_{h}^{T} \int_{\Omega_m^{\varepsilon}} \big| \Theta_{\ell,m}^{\varepsilon}(x, t) - \Theta_{\ell,m}^{\varepsilon}(x, t-h)\big|^2 \,dx\,dt + {\mathsf j}^{\varepsilon,\delta}_m \end{equation} with $\big|{\mathsf j}^{\varepsilon,\delta}_m \big| \leqslant C\delta$, where $C$ is a constant that does not depend on ${\varepsilon,\delta}$, $h$. Then we obtain \begin{equation} \label{point-8} I_1^{\varepsilon}(\varphi^{{\varepsilon,\delta},h}) \geqslant \Phi_{f} \int_{h}^{T} \int_{\Omega_{f}^{\varepsilon}}\big| \Theta_{\ell,f}^{\varepsilon}(x, t) - \Theta_{\ell,f}^{\varepsilon}(x, t-h)\big|^2\,dx\,dt + {\mathsf j}^{\varepsilon,\delta}_f + {\mathsf j}^{\varepsilon,\delta}_m. \end{equation} Next, we estimate the right-hand side of \eqref{eq:1.2.2-lem}. To this end we use the notation: $$ \int_{\Omega_{T}} K^{\varepsilon}(x) \widetilde \rho_{\ell}^{\varepsilon} \Big\{ \lambda_{\ell}(S^{\varepsilon})\big(\nabla P^{\varepsilon} - \widetilde \rho_{\ell}^{\varepsilon}\vec{g}\big) + \nabla\beta(S^{\varepsilon})\Big\} \cdot \nabla \varphi^{{\varepsilon,\delta},h}\,dx\,dt =: I_{f}^{\varepsilon}(\varphi^{{\varepsilon,\delta},h}) + I_m^{\varepsilon}(\varphi^{{\varepsilon,\delta},h}). $$ We have \begin{equation} \label{eq:1.5.8} \begin{aligned} &\big| I_{f}^{\varepsilon}(\varphi^{{\varepsilon,\delta},h})\big| + \big|I_m^{\varepsilon}(\varphi^{{\varepsilon,\delta},h})\big| \\ &\leqslant \| K_{f}\widetilde\rho_{\ell}^{\varepsilon} \big\{\lambda_{\ell}(S_{f}^{\varepsilon})\nabla P_{f}^{\varepsilon} + \nabla\beta(S_{f}^{\varepsilon}) - \widetilde\rho_{\ell}^{\varepsilon} \lambda_{\ell}(S_{f}^{\varepsilon}) \vec{g}\big\} \|_{L^2(\Omega^{\varepsilon}_{f,T})} \| \nabla\varphi^{{\varepsilon,\delta},h} \|_{L^2(\Omega^{\varepsilon}_{f,T})} \\ &\quad + \varepsilon \| K_m \widetilde\rho_{\ell}^{\varepsilon} \big\{\lambda_{\ell}(S_m^{\varepsilon})\nabla P_m^{\varepsilon} + \nabla\beta(S_m^{\varepsilon}) - \widetilde\rho_{\ell}^{\varepsilon} \lambda_{\ell}(S_m^{\varepsilon}) \vec{g}\big\} \|_{L^2(\Omega^{\varepsilon}_{m,T})}\\ &\quad\times \varepsilon\| \nabla \varphi^{{\varepsilon,\delta},h} \|_{L^2(\Omega^{\varepsilon}_{m,T})}. \end{aligned} \end{equation} Then for $I_{f}^{\varepsilon}(\varphi^{{\varepsilon,\delta},h})$, $I_m^{\varepsilon}(\varphi^{{\varepsilon,\delta},h})$ from (A3), (A5), the \emph{a priori} estimates \eqref{beta-un}, \eqref{eq:1.5.5}, and \eqref{eq:1.5.5+1} we obtain \begin{equation} \label{eq:1.5.9} \big| I_{f}^{\varepsilon}(\varphi^{{\varepsilon,\delta},h})\big| + \big| I_m^{\varepsilon}(\varphi^{{\varepsilon,\delta},h})\big| \leqslant C_1\, h\, \delta^{-\eta}. \end{equation} From \eqref{eq:1.2.2-lem}, \eqref{point-8}, and \eqref{eq:1.5.9}, we obtain the bound \begin{equation} \label{new-goods-6} \int_{h}^{T} \int_{\Omega_{f}^{\varepsilon}}\big| \Theta_{\ell,f}^{\varepsilon}(x, t) - \Theta_{\ell,f}^{\varepsilon}(x, t-h)\big|^2\,dx\,dt \leqslant C_1 h \delta^{-\eta} + C_2 \delta. \end{equation} Now we set: $\delta = \delta(h) = h^{1/2\eta}$. Then from \eqref{new-goods-6} we have \begin{equation} \label{nouv-idee-18} \int_{h}^{T} \int_{\Omega^\varepsilon_f}\, \big| \Theta_{\ell,f}^{\varepsilon}(x, t) - \Theta_{\ell,f}^{\varepsilon}(x, t-h)\big|^2\,dx\,dt \leqslant C\, h^\sigma \quad \text{with } \sigma := \min\{\frac{1}{2}, \frac{1}{2\eta}\}. \end{equation} Finally, taking into account the fact that the extension is by reflection, from \eqref{nouv-idee-18}, we obtain the desired estimate \eqref{point-000} and Lemma \ref{lem11} is proved. \end{proof} By the Kolmogorov-Riesz-Fr\'echet compactness criterion, the compactness of the sequences $\{\widetilde\Theta_{\ell,f}^{\varepsilon}\}_{\varepsilon>0}$, $\{\widetilde\Theta_{g,f}^{\varepsilon}\}_{\varepsilon>0}$ is a consequence of Lemmata \ref{lem10}, \ref{lem11} (evidently, with the corresponding change of the index $"\ell"$ by $"g"$ for the sequence $\{\widetilde\Theta_{g,f}^{\varepsilon}\}_{\varepsilon>0}$). It assures the following convergence results: \begin{gather} \label{L-100} \widetilde\Theta_{\ell,f}^{\varepsilon} \to {\mathsf L_1} \text{ strongly in } L^1(\Omega_T), \\ \label{L-200} \widetilde\Theta_{g,f}^{\varepsilon} \to {\mathsf L_2} \text{ strongly in } L^1(\Omega_T). \end{gather} Now taking into account the $L^\infty$-boundedness of the functions $\widetilde\Theta_{\ell,f}^{\varepsilon}$, $\widetilde\Theta_{g,f}^{\varepsilon}$ we, finally, arrive to the desired convergence results \eqref{L-1}, \eqref{L-2}. This completes the proof of Proposition \ref{prop9}. \begin{corollary} \label{cor16} There exist the functions $P, S \in L^2(\Omega_{T})$ such that (up to a subsequence), we have \begin{gather} \widetilde P_{f}^{\varepsilon} \to P \quad \text{a.e in } \Omega_{T}; \label{eq:1.5.13} \\ \widetilde S_{f}^{\varepsilon} \to S \quad \text{a.e in } \Omega_{T}. \label{eq:1.5.14} \end{gather} \end{corollary} \begin{proof} By Proposition \ref{prop9} we conclude that there exist the functions ${\mathsf L_1}, {\mathsf L_2}$ such that \begin{gather*} \widetilde\Theta_{\ell,f}^{\varepsilon} \to {\mathsf L_1} \text{ strongly in $L^2(\Omega_T)$ and a.e. in } \Omega_T; \\ \widetilde \Theta_{g,f}^{\varepsilon} \to {\mathsf L_2} \text{ strongly in $L^2(\Omega_T)$ and a.e. in } \Omega_T. \end{gather*} Now, we will use the fact that the function ${\EuScript D}$ given by $$ {\EuScript D}(\widetilde S_{f}^{\varepsilon}, \widetilde P_{f}^{\varepsilon}) = (\widetilde\Theta_{\ell,f}^{\varepsilon}(\widetilde S_{f}^{\varepsilon}, \widetilde P_{f}^{\varepsilon}), \widetilde\Theta_{g,f}^{\varepsilon}(\widetilde S_{f}^{\varepsilon}, \widetilde P_{f}^{\varepsilon})) $$ is a diffeomorphism from $[0,1] \times \mathbb{R}$ to ${\EuScript D}([0,1]\times\mathbb{R})$ (for more details see \cite{AJZ1}) so it has continuous inverse. Therefore, almost everywhere in $\Omega_{T}$ convergence of $\widetilde\Theta_{r,f}^{\varepsilon}(\widetilde S_{f}^{\varepsilon}, \widetilde P_{f}^{\varepsilon})$ ($r\in \{ \ell,g\} $) implies \eqref{eq:1.5.13} and \eqref{eq:1.5.14}. Corollary \ref{cor16} is proved. \subsection{Two--scale convergence results} \label{con-results} In this section, taking into account the compactness results from the previous section, we formulate the convergence results for the sequences $\{\widetilde P^\varepsilon_{f}\}_{\varepsilon>0}$, $\{\widetilde S^\varepsilon_{f}\}_{\varepsilon>0}$, $\{\widetilde \Theta_{r,f}^{\varepsilon} \}_{\varepsilon>0}$ with $r\in \{ \ell,g\} $. In this paper the homogenization process for the problem is rigorously obtained by using the two-scale approach, see, e.g., \cite{Allaire}. For reader's convenience, let us recall the definition of the two-scale convergence. \begin{definition} \label{two} \rm A sequence of functions $\{v^\varepsilon\}_{\varepsilon>0} \subset L^2(\Omega_T)$ two-scale converges to $v \in L^2(\Omega_T\times Y)$ if $\| v^\varepsilon\|_{L^2(\Omega_T)} \leqslant C$, and for any test function $\varphi \in C^\infty(\Omega_T; C_{\rm per}(Y))$ the following relation holds $$ \lim_{\varepsilon\to 0} \int_{\Omega_T} v^\varepsilon(x, t)\, \varphi (x, \frac{x}{\varepsilon}, t)\, dx\, dt = \int_{\Omega_T \times Y} v(x, y, t)\, \varphi(x, y, t) \, dy\, dx\, dt. $$ \end{definition} This convergence is denoted by $v^\varepsilon(x, t) \stackrel {2s} \rightharpoonup v(x, y, t)$. Now we summarize the convergence results for the sequences $\{\widetilde P^\varepsilon_{f}\}_{\varepsilon>0}$, $\{\widetilde S^\varepsilon_{f}\}_{\varepsilon>0}$, $\{\widetilde{\Theta_{r,f}^{\varepsilon}} \}_{\varepsilon>0}$ ($r\in \{ \ell,g\}$). They are given by the following lemma. \begin{lemma} \label{thm17} There exist a function $S$ such that $0 \leqslant S \leqslant 1$ a.e. in $\Omega_T$, and functions $P \in L^2(0, T; H^1(\Omega))$, $P_1, \beta_1 \in L^2(\Omega_T; H^1_{\rm per}(Y))$ such that up to a subsequence, \begin{gather} \widetilde S_{f}^{\varepsilon} \to S \text{ strongly in $ L^2(\Omega_T)$ and a.e. in } \Omega_T; \label{eq:1.6.3} \\ \beta(\widetilde S_{f}^{\varepsilon}) \to \beta(S) \text{ strongly in } L^2(\Omega_T); \label{eq:1.6.4} \\ \widetilde P_{f}^{\varepsilon} \to P \text{ strongly in $L^2(\Omega_T)$ and a.e. in } \Omega_T; \label{eq:1.5.13-1} \\ \widetilde P_{f}^{\varepsilon} \rightharpoonup P \text{ weakly in } L^2(0,T;H^{1}(\Omega)); \label{eq:1.6.1} \\ \widetilde\Theta_{\ell,f}^{\varepsilon} \rightharpoonup \rho_{\ell} (P + {\mathsf G}_{\ell}\,(S))S \text{ strongly in } L^2(\Omega_T); \label{eq:1.6.2} \\ \widetilde\Theta_{g,f}^{\varepsilon} \rightharpoonup \rho_g(P + {\mathsf G}_g(S))(1 - S) \text{ strongly in } L^2(\Omega_T); \label{eq:1.6.2-1} \\ \nabla\widetilde P_{f}^{\varepsilon}(x,t)\; \underrightarrow{2s}\; \nabla P(x,t)+\nabla_{y}P_1(x,t,y); \label{eq:1.6.5} \\ \nabla\beta(\widetilde S_{f}^{\varepsilon})(x,t) \; \underrightarrow{2s}\; \nabla\beta(S)(x,t) + \nabla_{y}\beta_1(x,t,y). \label{eq:1.6.6} \end{gather} \end{lemma} The proof of the above lemma is based on the a priori estimates obtained in Section \ref{uni-est}, Proposition \ref{prop9}, Corollary \ref{cor16}, and two-scale convergence arguments similar to those in \cite{Allaire} (see also Lemma 4.8 from \cite{Siam} for similar arguments). \section{Dilation operator and convergence results} \label{dil-oper} It is known that by the nonlinearities and the strong coupling of the problem, the two--scale convergence does not provide an explicit form for the source terms appearing in the homogenized model, see for instance \cite{blm,choq,yeh2}. To overcome this difficulty the authors make use of the {dilation operator}. Here we refer to \cite{adh,blm,choq,yeh2} for the definition and main properties of the dilation operator. Let us also notice that the notion of the dilation operator is closely related to the notion of the {unfolding operator}. We refer here, e.g., to \cite{ddg} for the definition and the properties of this operator. The outline of the section is as follows. First, in Section \ref{def-base-dil} we introducethe definition of the dilation operator and describe its main properties. Then in Section \ref{dil-func-prop} we obtain the equations for the dilated saturation and the global pressure functions and the corresponding uniform estimates. Finally, in Section \ref{dil-func-conv} we consider the convergence results for the dilated functions. \subsection{Definition of the dilation operator and its main properties} \label{def-base-dil} \begin{definition} \label{def-dilop} \rm For a given $\varepsilon > 0$, we define a dilation operator $D^\varepsilon$ mapping measurable functions defined in $\Omega^\varepsilon_{m,T}$ to measurable functions defined in $\Omega_T \times Y_m$ by \begin{equation} \label{do-1} (D^\varepsilon \psi )(x, y, t) := \begin{cases} \psi( c^\varepsilon(x) + \varepsilon y, t), & \text{if } c^\varepsilon(x) + \varepsilon y \in \Omega^\varepsilon_m; \\ 0, & \text{otherwise}, \end{cases} \end{equation} where $c^\varepsilon(x) := \varepsilon k$ if $x \in \varepsilon (Y_m + k)$. Here $k \in \mathbb{Z}^d$ denotes the lattice translation point of the $\varepsilon$-cell domain which contains $x$. \end{definition} The basic properties of the dilation operator are given by the following lemma (see \cite{adh,yeh2}). \begin{lemma} \label{lemma-dilop1} Let $v, w \in L^2(0, T; H^1(\Omega^\varepsilon_m))$. Then we have \begin{gather} \label{dilll} \nabla_y D^\varepsilon v = \varepsilon D^\varepsilon (\nabla_x v) \quad \text{a.e. in } \Omega_T \times Y_m; \\ \| D^\varepsilon v \|_{L^2(\Omega_T \times Y_m)} = \| v \|_{L^2(\Omega_T)}; \nonumber \\ \| \nabla_y D^\varepsilon v \|_{L^2(\Omega_T \times Y_m)} = \varepsilon \| D^\varepsilon \nabla_x\, v \|_{L^2(\Omega^\varepsilon_{m,T})}; \nonumber \\ (D^\varepsilon v, D^\varepsilon w )_{L^2(\Omega_T \times Y_m)} = (v, w)_{L^2(\Omega^\varepsilon_{m,T})}; \nonumber \\ (D^\varepsilon v, w )_{L^2(\Omega_T \times Y_m)} = (v, D^\varepsilon w )_{L^2(\Omega_T \times Y_m)}. \nonumber \end{gather} \end{lemma} The following lemma gives the link between the two-scale and the weak convergence (see, e.g., \cite{blm}). \begin{lemma} \label{lemma-dilop} Let $\{v^\varepsilon\}_{\varepsilon>0}$ be a uniformly bounded sequence in $L^2(\Omega^\varepsilon_{m,T})$ that satisfies the conditions: \begin{itemize} \item[(i)] $D^\varepsilon v^\varepsilon \rightharpoonup v^0$ weakly in $L^2(\Omega_T; L^2_{\rm per}(Y_m))$; \item[(ii)] $\mathbf{1}^\varepsilon_m(x) v^\varepsilon \stackrel {2s} \rightharpoonup v^\star \in L^2(\Omega_T; L^2_{\rm per}(Y_m))$. \end{itemize} Then $v^0 = v^\star$ a.e. in $\Omega_T \times Y_m$. \end{lemma} Also we have the following result (see, e.g., \cite{choq,yeh2}). \begin{lemma} \label{lemma-dilop2} If $v^\varepsilon \in L^2(\Omega^\varepsilon_{m,T})$ and $\mathbf{1}^\varepsilon_m(x) v^\varepsilon \stackrel {2s} \to v \in L^2(\Omega_T; L^2_{\rm per}(Y_m))$ then $D^\varepsilon v^\varepsilon$ converges to $v$ strongly in $L^q(\Omega_T \times Y_m)$, where $"\stackrel {2s} \to"$ denotes the strong two--scale convergence. If $v \in L^2(\Omega_T)$ is considered as an element of $L^2(\Omega_T \times Y_m)$ constant in $y$, then $D^\varepsilon v$ converges strongly to $v$ in $L^2(\Omega_T \times Y_m)$. \end{lemma} \subsection{Dilated functions $D^\varepsilon S^\varepsilon_m, D^\varepsilon P^\varepsilon_m$ and their properties} \label{dil-func-prop} Let us introduce the following notation for the dilated functions defined on the matrix part of the reservoir $\Omega$: \begin{equation} p_m^{\varepsilon} := D^{\varepsilon}P_m^{\varepsilon}, \quad s_m^{\varepsilon} := D^{\varepsilon}S_m^{\varepsilon}, \quad \vartheta_m^{\varepsilon} := D^{\varepsilon}\beta_m^{\varepsilon} = \beta(s_m^{\varepsilon}), \label{eq:1.9.1} \end{equation} and \begin{equation} \begin{gathered} \theta_{\ell,m}^{\varepsilon} := D^{\varepsilon}\Theta_{\ell,m}^{\varepsilon} = \rho_{\ell}(p_m^{\varepsilon} + G_{\ell}\,(s_m^{\varepsilon}))\,s_m^{\varepsilon} , \\ \theta_{g,m}^{\varepsilon} := D^{\varepsilon}\Theta_{g,m}^{\varepsilon} = \rho_g(p_m^{\varepsilon} + G_g\,(s_m^{\varepsilon}))\,(1 - s_m^{\varepsilon}). \end{gathered} \label{eq:1.9.1+1} \end{equation} The goal of this section is to derive the equations for the dilated functions $s_m^{\varepsilon}, p_m^{\varepsilon}$ and using Lemma \ref{lemma-dilop1} to obtain the corresponding uniform estimates. First, we establish the following result. \begin{lemma}\label{lem:25} For any $\varepsilon > 0$ and for a.e. $x \in \Omega$, the dilated functions $s_m^{\varepsilon}, p_m^{\varepsilon}$ defined by \eqref{eq:1.9.1} are the solutions of the system \begin{equation} \begin{gathered} 0 \leqslant s_m^{\varepsilon} \leqslant 1; \\ \begin{aligned} &\Phi_m \frac{\partial\theta_{\ell,m}^{\varepsilon}}{\partial t} - \operatorname{div}_y\Big\{ K_m \rho_{\ell}\,(s_m^{\varepsilon},p_m^{\varepsilon}) \big[\lambda_{\ell}(s_m^{\varepsilon})\nabla p_m^{\varepsilon} \\ &+ \nabla\vartheta_m^{\varepsilon} - \varepsilon \lambda_{\ell}(s_m^{\varepsilon})\rho_{\ell}\,(s_m^{\varepsilon},p_m^{\varepsilon}) \vec{g}\big] \Big\} = 0; \end{aligned} \\ \begin{aligned} &\Phi_m\frac{\partial\theta_{g,m}^{\varepsilon}}{\partial t} - \operatorname{div}_y\Big\{ K_m\, \rho_g(s^{\varepsilon},p_m^{\varepsilon}) \big[\lambda_g(s_m^{\varepsilon})\nabla p_m^{\varepsilon} \\ &- \nabla \vartheta_m^{\varepsilon} - \varepsilon \lambda_g(s_m^{\varepsilon}) \rho_g(s_m^{\varepsilon},p_m^{\varepsilon})\vec{g}\big]\Big\} = 0 \end{aligned}\\ \end{gathered} \label{eq:1.9.2} \end{equation} in $L^2(0,T;H^{-1}(Y_m))$. Here $\rho_{\rm k}\,(s_m^{\varepsilon},p_m^{\varepsilon}) := \rho_{\rm k}(p_m^{\varepsilon} + G_{\rm k}\,(s_m^{\varepsilon}))$ (${\rm k} = \ell, g$). \end{lemma} The proof of the above can be done by the arguments similar to those used in the proof of \cite[Lemma 6.7]{yeh2}. We complete system \eqref{eq:1.9.2} with the boundary and initial conditions. The boundary conditions are \begin{equation} p_m^{\varepsilon}(x, y, t) = P_{f}^{\varepsilon}(\varepsilon y + c^{\varepsilon}(x), t), \quad \vartheta_m^{\varepsilon}(x, y, t) = \beta_{f}^{\varepsilon}(\varepsilon y + c^{\varepsilon}(x), t) \label{eq:1.9.5} \end{equation} in $H^{1/2}(\Gamma_{fm})$ for $(x, t) \in \Omega^{\varepsilon}_{m,T}$, where $\beta_{f}^{\varepsilon} := \beta(S^\varepsilon_{f})$. The initial conditions are \begin{equation} p_m^{\varepsilon}(x, y, 0) = D^{\varepsilon} P^{0}_m(x, y), \quad s_m^{\varepsilon}(x , y, 0) = D^{\varepsilon} S^{0}_m(x, y) \quad \text{in } \Omega\times Y_m. \label{eq:1.9.6} \end{equation} For $x$ outside the matrix part by the definition \eqref{do-1}, we set \begin{equation} p_m^{\varepsilon}(x, y, t) = s_m^{\varepsilon}(x, y, t) = 0\quad \text{in } H^{1/2}(\Gamma_{fm}) \quad\text{for } x \in \Omega \setminus \Omega_m^{\varepsilon},\; t \in (0, T). \label{eq:1.9.7} \end{equation} \begin{remark} \label{rem26} \rm There is no confusion that in \eqref{eq:1.9.6} we use again the notation $D^{\varepsilon}$, for the dilatation operator defined for $L^2(\Omega)$-functions which maps them to the functions in $L^2(\Omega\times Y_m)$ by the formula \begin{equation} (D^{\varepsilon}\varphi)(x, y) := \varphi(\varepsilon y+c^{\varepsilon}(x)). \label{eq:1.9.8} \end{equation} \end{remark} \begin{remark} \label{rem26+1} \rm We note that problem \eqref{eq:1.9.2} with non-homogeneous Dirichlet boundary conditions \eqref{eq:1.9.5} and initial conditions \eqref{eq:1.9.6} corresponds to a family of problems in $Y_m\times(0, T)$ parameterized by the slow variable $x$ and depending on it through the boundary data \eqref{eq:1.9.5}. \end{remark} The following lemma contains the uniform \emph{a priori} estimates for the dilated solutions $s_m^{\varepsilon}, p_m^{\varepsilon}$ of \eqref{eq:1.9.2}. \begin{lemma} \label{lem:27} Let $(p_m^{\varepsilon}, s^\varepsilon_m)$ be a solution to \eqref{eq:1.9.2}-\eqref{eq:1.9.6}. There exists a constant $C$ that does not depend on $\varepsilon$ and such that \begin{gather} \label{dilo-8} 0 \leqslant s^\varepsilon_m \leqslant 1 \quad \text{a.e. in } \Omega_T \times Y_m; \\ \| p_m^{\varepsilon} \|_{L^2(\Omega_{T}; H_{\rm per}^{1}(Y_m))} + \| \vartheta_m^{\varepsilon} \|_{L^2(\Omega_{T}; H_{\rm per}^{1}(Y_m))} \leqslant C; \label{eq:1.9.9} \\ \label{dilo-11} \| \partial_t(\Phi_m\, \theta^\varepsilon_{\ell,m}) \|_{L^2(\Omega_T; H^{-1}_{\rm per}(Y_m))} + \| \partial_t(\Phi_m\,\theta^\varepsilon_{g,m}) \|_{L^2(\Omega_T; H^{-1}_{\rm per}(Y_m))} \leqslant C; \end{gather} where $\vartheta_m^{\varepsilon} := D^{\varepsilon}\beta_m^{\varepsilon} = \beta(s_m^{\varepsilon})$. \end{lemma} \begin{proof} Statement \eqref{dilo-8} is evident. The uniform estimate \eqref{eq:1.9.9} easy follow from the uniform bounds \eqref{beta-un}, \eqref{eq:1.3.26}, and Lemma \ref{lemma-dilop1}. The bound \eqref{dilo-11} follows from Lemmata \ref{lemma-dilop1}, \ref{cor-ps-l5}. The proof is complete. \end{proof} \subsection{Convergence results for the dilated functions $p_m^{\varepsilon}, s^\varepsilon_m, \theta_{\ell,m}^{\varepsilon}, \theta_{g,m}^{\varepsilon}$} \label{dil-func-conv} In this section we establish convergence results which will be used below to obtain the homogenized system. From Lemmata \ref{lemma-dilop}, \ref{lem:27} we get the following convergence results. \begin{lemma} \label{cor:28} Let $(p_m^{\varepsilon}, s^\varepsilon_m)$ be a solution to \eqref{eq:1.9.2}--\eqref{eq:1.9.6}. Then (up to a subsequence), \begin{gather} \label{weak-dil-10} \chi_m^{\varepsilon}\, S^\varepsilon_m \stackrel {2s} \rightharpoonup s \in L^2(\Omega_T; L^2_{\rm per}(Y_m)), \quad s^\varepsilon_m \rightharpoonup s \text{ weakly in } L^2(\Omega_T \times Y_m); \\ \chi_m^{\varepsilon}\, P^\varepsilon_m \stackrel {2s} \rightharpoonup p \in L^2(\Omega_T; L^2_{\rm per}(Y_m)), \quad p^\varepsilon_m \rightharpoonup p \text{ weakly in } L^2(\Omega_T; H^1(Y_m)); \label{eq:1.6.7} \\ \varepsilon \nabla_x P^\varepsilon_m \stackrel {2s} \rightharpoonup \nabla_y p \in L^2(\Omega_T; L^2_{\rm per}(Y_m)); \label{eq:1.6.10} \\ \begin{gathered} \chi_m^{\varepsilon} \beta(S^\varepsilon_m) \stackrel {2s} \rightharpoonup \vartheta^\star; \quad \beta(s^\varepsilon_m) \rightharpoonup \vartheta^\star \text{ weakly in } L^2(\Omega_T; H^1(Y_m)); \\ \varepsilon \nabla_x \beta(S_m^{\varepsilon}) \stackrel {2s} \rightharpoonup \nabla_{y}\vartheta^\star, \end{gathered}\label{eq:1.6.8} \end{gather} where $\chi_m^\varepsilon = \chi_m^\varepsilon(x)$ is the characteristic function of the subdomain $\Omega_m^\varepsilon$. \end{lemma} \begin{remark} \label{rmk3} \rm The limit $\vartheta^\star$ will be identified below. More precisely, it will be shown that $\vartheta^\star = \beta(s)$. \end{remark} It is clear that the convergence results of Lemma \ref{cor:28} above are not sufficient for derivation of the equations for the limit functions $s, p$. In order to overcome this difficulty, we introduce the restrictions of the dilated functions $s^\varepsilon_m$, $p^\varepsilon_m$ which are defined below. The key feature of the newly introduced functions is that they possess more compactness properties than the dilated functions introduced in \eqref{eq:1.9.1}, \eqref{eq:1.9.1+1}. To this end, following \cite{choq}, we fix $x = x_0 \in \Omega$ and define the restrictions of $s^\varepsilon_m$, $p^\varepsilon_m$ to the $\varepsilon$--cell containing the point $x_0$. These functions are defined in the domain $Y_m \times (0, T)$ and are constants in the slow variable $x$. In order to obtain the uniform estimates for the restricted functions (they are similar to the corresponding estimates for $P^\varepsilon_{f}$, $S^\varepsilon_{f}$ from Section \ref{uni-est}) we make use of the estimates \eqref{dilo-8}--\eqref{eq:1.9.9}. Then we proceed as follows. First, we define the restrictions of $s^\varepsilon_m$, $p^\varepsilon_m$ to a cell $K^\varepsilon_{x_0}$ which is a cube containing $x_0 \in \Omega$ and the length $\varepsilon$. We denote by $k(x_0, \varepsilon) \in \mathbb{Z}^d$ such that $K^\varepsilon_{x_0} := \varepsilon \big(Y + k(x_0, \varepsilon)\big)$. Due to the definition of the dilation operator, the functions $s^\varepsilon_m$, $p^\varepsilon_m$ are constant in $x$ in $K^\varepsilon_{x_0}$. The restricted functions are given by: \begin{equation} \label{new-smx} s^\varepsilon_{m,x_0}(y, t) := \begin{cases} s^\varepsilon_m, &\text{if } x \in K^\varepsilon_{x_0}; \\ 0, & \text{otherwise}; \end{cases} \quad p^\varepsilon_{m,x_0}(y, t) := \begin{cases} p^\varepsilon_m, &\text{if } x \in K^\varepsilon_{x_0}; \\ 0, &\text{otherwise} . \end{cases} \end{equation} For any $\varepsilon > 0$, the pair of functions $(s^\varepsilon_{m,x_0}, p^\varepsilon_{m,x_0})$ is a solution to \eqref{eq:1.9.2}--\eqref{eq:1.9.6} in $Y_m \times (0,T)$. Now we obtain the uniform, in $\varepsilon$, for the functions $s^\varepsilon_{m,x_0}, p^\varepsilon_{m,x_0}$ defined in \eqref{new-smx}, for a.e. $x_0 \in \Omega$. We have. \begin{lemma} \label{lem:29} For a.e. $x_{0} \in \Omega$, there is a constant $C = C(x_{0})$ which does not depend on $\varepsilon$ and such that \begin{gather} \| p_{m,x_{0}}^{\varepsilon}\|_{L^2(0,T;H_{\rm per}^{1}(Y_m))} + \|\vartheta_{m,x_{0}}^{\varepsilon}\|_{L^2(0,T;H_{\rm per}^{1}(Y_m))} \leqslant C(x_{0}), \label{eq:1.9.16} \\ \label{dilo-11-x0} \| \partial_t(\Phi_m\, \theta^\varepsilon_{\ell,m,x_{0}}) \|_{L^2(0,T; H^{-1}_{\rm per}(Y_m))} +\| \partial_t(\Phi_m\,\theta^\varepsilon_{g,m,x_{0}}) \|_{L^2(0,T; H^{-1}_{\rm per}(Y_m))} \leqslant C; \end{gather} where \begin{gather*} \theta_{\ell,m,x_{0}}^{\varepsilon} := \rho_{\ell}(p_{m,x_{0}}^{\varepsilon} + G_{\ell}\,(s_{m,x_{0}}^{\varepsilon}))\,s_{m,x_{0}}^{\varepsilon}, \\ \theta_{g,m,x_{0}}^{\varepsilon} := \rho_g(p_{m,x_{0}}^{\varepsilon} + G_g\,(s_{m,x_{0}}^{\varepsilon}))\,(1 - s_{m,x_{0}}^{\varepsilon}). \end{gather*} \end{lemma} Next we establish a compactness for the families $\{\theta_{\ell,m,x_{0}}^{\varepsilon}\}_{\varepsilon>0}$ and $\{ \theta_{g,m,x_{0}}^{\varepsilon}\}_{\varepsilon>0}$. We obtain this compactness result by applying \cite[Lemma 4.2]{Siam}. \begin{proposition} \label{lem:31} Under our standing assumptions, for a.e. $x_{0}\in\Omega$, the families $\{\theta_{w,m,x_{0}}^{\varepsilon}\}_{\varepsilon>0}$ and $\{\theta_{g,m,x_{0}}^{\varepsilon}\}_{\varepsilon>0}$ are compact sets in $L^2(Y_m\times(0,T))$. \end{proposition} \section{Formulation of the main result} \label{main-res} We study the asymptotic behavior of the solution to problem \eqref{eq:1.1.15} (through the equivalent problem \eqref{eq:1.1.34}) as $\varepsilon \to 0$. In particular, we are going to show that the effective model is \begin{equation}\label{hom-0} \begin{gathered} 0 \leqslant S \leqslant 1 \quad \text{in } \Omega_T; \\ \Phi^\star \frac{\partial \Xi^\star_{\ell}}{\partial t} - \operatorname{div}_x \Big\{\mathbb{K}^\star\, \rho_{\ell}\,\lambda_{\ell}(S) \big[ \nabla P_{\ell} - \rho_{\ell}\vec{g}\big] \Big\} = {\EuScript Q}_{\ell} \quad \text{in } \Omega_T; \\ \Phi^\star \frac{\partial \Xi^\star_g}{\partial t} - \operatorname{div}_x \Big\{\mathbb{K}^\star\, \rho_g\,\lambda_g(S) \big[ \nabla P_g - \rho_g\vec{g}\big] \Big\} = {\EuScript Q}_g \quad \text{in } \Omega_T;\\ P_c(S) = P_g - P_{\ell} \quad \text{in } \Omega_T. \end{gathered} \end{equation} Here we use the following notation: \noindent$\bullet$ $S$, $P_{\ell}$, $P_g$ denote the homogenized water saturation, water pressure, and gas pressure, respectively. \noindent$\bullet$ $\Phi^\star$ denotes the effective porosity and is given by \begin{equation} \label{hom-1} \Phi^\star := \Phi_f\frac{|Y_f|}{|Y_m|}, \end{equation} where $|Y_r|$ is the measure of the set $Y_r$, $r = f, m$. \noindent$\bullet$ $\mathbb{K}^\star$ is the homogenized tensor with the entries $\mathbb{K}^\star_{ij}$ defined by: \begin{equation} \label{hom-2} \mathbb{K}^{\star}_{ij} := \frac{K_f}{|Y_m|}\, \int_{Y_\mathsf{f}}\, [\nabla_y \xi_i + \vec e_i ] [\nabla_y \xi_j + \vec e_j ]\, dy, \end{equation} where $\xi_j$ is a solution of the auxiliary problem \begin{equation} \label{hom-3} \begin{gathered} - \Delta_y \xi_j = 0 \quad \text{in } Y_\mathsf{f}; \\ \nabla_y \xi_j \cdot \vec \nu = - \vec e_j \cdot \vec \nu \quad \text{on } \Gamma_{fm};\\ y \mapsto \xi_j(y)\quad Y\text{-periodic}. \end{gathered} \end{equation} \noindent$\bullet$ The functions $\Xi^\star_{\ell}, \Xi^\star_g$ are \begin{equation} \label{hom-4} \Xi^\star_{\ell} := S \rho_{\ell}\,(P_{\ell}) \quad \text{and} \quad \Xi^\star_g := (1 - S) \rho_g\,(P_g). \end{equation} For almost all point $x \in \Omega$ a matrix block $Y_\mathsf{m} \subset \mathbb{R}^d$ is suspended topologically. Equations for flow in a matrix block are given by \begin{equation}\label{hom-5} \begin{gathered} 0 \leqslant s \leqslant 1 \quad \text{in } \Omega_T\times Y_m; \\ \Phi_m \frac{\partial \theta^\star_\ell}{\partial t} - \operatorname{div}_y \Big\{K_m\, \rho_{\ell} (p_{\ell}) \lambda_{\ell}\,(s) \nabla_y p_{\ell} \Big\} = 0 \quad \text{in } \Omega_T\times Y_m; \\ \Phi_m \frac{\partial \theta^\star_g}{\partial t} -\operatorname{div}_y \Big\{K_m \rho_g(p_g) \lambda_g(s) \nabla_y p_g \Big\} = 0 \quad \text{in } \Omega_T\times Y_m;\\ P_c(s) = p_g - p_{\ell} \quad \text{in } \Omega_T\times Y_m. \end{gathered} \end{equation} Here we use the following notation: \noindent$\bullet$ $s$, $p_{\ell}$, $p_g$ denote the water saturation, the water and gas pressures in the block $Y_m$, respectively. \noindent$\bullet$ The functions $\theta^\star_\ell, \theta^\star_g$ are defined by \begin{equation} \label{hom-6} \theta^\star_\ell := s \rho_{\ell}(p_{\ell}) \quad \text{and} \quad \theta^\star_g := (1 - s) \rho_g(p_g). \end{equation} For any $x \in \Omega$ and $t > 0$, the matrix-fracture sources are given by \begin{equation}\label{hom-7} {\EuScript Q}_{\ell} := - \frac{\Phi_m}{|Y_m|}\, \int_{Y_m} \frac{\partial \theta^\star_\ell}{\partial t} (x, y, t) \,dy \quad \text{and} \quad {\EuScript Q}_g := - \frac{\Phi_m}{|Y_m|}\, \int_{Y_m} \frac{\partial \theta^\star_g}{\partial t}(x, y, t) \,dy. \end{equation} The boundary conditions for the effective system \eqref{hom-0} are given by \begin{equation} \label{hom-8} \begin{gathered} P_g(x, t) = P_{\ell}(x, t) = 0 \quad \text{on } \Gamma_1 \times (0,T); \\ \mathbb{K}^\star\rho_g(P_g) \lambda_g(S) \big(\nabla P_g - \rho_g(P_g) \vec{g} \big) \cdot \vec \nu = 0 \quad \text{on } \Gamma_{2} \times (0,T);\\ \mathbb{K}^\star \rho_{\ell}(P_{\ell}) \lambda_{\ell}\,(S) \big(\nabla P_\ell - \rho_{\ell}(P_{\ell}) \vec{g} \big)\cdot \vec \nu = 0 \quad \text{on } \Gamma_{2} \times (0,T). \end{gathered} \end{equation} The interface conditions for the system \eqref{hom-5} are given by: \begin{equation}\label{hom-9} \begin{gathered} P_{\ell}(x, t) = p_{\ell}(x, y, t) \quad \text{on } \Omega_T \times \Gamma_{fm}; \\ P_g(x, t) = p_g(x, y, t) \quad \text{on } \Omega_T \times \Gamma_{fm}. \end{gathered} \end{equation} Finally, the initial conditions read \begin{gather} \label{hom-10} S(x, 0) = S^0(x) \quad \text{and} \quad P_g(x, 0) = p_g^{0}(x) \quad \text{in } \Omega, \\ \label{hom-11} s(x, y, 0) = S^0(x) \quad \text{and} \quad p_g(x, y, 0) = p_g^{0}(x) \quad \text{in } \Omega \times Y_m. \end{gather} The main result of the paper reads as follows. \begin{theorem} \label{t-hom-main} Let assumptions {\rm (A1)--(A8)} be fulfilled. Then the solution of initial problem \eqref{eq:1.1.15} converges (up to a subsequence) in the two-scale sense to a weak solution of homogenized problem \eqref{hom-0}. \end{theorem} \section{Proof of Theorem \ref{t-hom-main}} \label{proof-main} The proof of the homogenization result will be done in several steps. Using the convergence and compactness results from Section \ref{ex-comp-sf} we pass to the limit in equations \eqref{eq:1.2.2}, \eqref{eq:1.2.3}. This is done in Section \ref{passage1}. The passage to the limit in the matrix blocks makes use of the dilation operator defined in Section \ref{dil-oper} above. In the passage to the limit in \eqref{eq:1.2.2}, \eqref{eq:1.2.3} we see that the homogenized equations contain, first of all, the corrector functions $P_1, \beta_1$ defined in Lemma \ref{thm17} as well as non--explicit limits ${\mathsf L}_{\ell,m}, {\mathsf L}_{g,m}$ (see Section \ref{passage1}). The corrector functions are then identified in Section \ref{ident-p-beta}. Having obtained the homogenized equations in terms of the global pressure and saturation, we can reformulate these equations in terms of new unknown functions that it is naturally to call the {homogenized phase pressures}. They are defined as follows: $P_{\ell} := P + G_{\ell}(S)$ and $P_g := P + G_g(S)$. However, in these equations the limits ${\mathsf L}_{\ell,m}, {\mathsf L}_{g,m}$ are still non identified. To obtain the explicit form of these limits we make use of the ideas from \cite{choq}. This completes the proof of the main result of the paper. \subsection{Passage to the limit in the equation \eqref{eq:1.2.2}, \eqref{eq:1.2.3}} \label{passage1} We set \begin{equation}\label{paslim-1} \begin{aligned} \varphi_{\ell}(x, \frac{x}{\varepsilon}, t ) :&= \varphi(x, t) + \varepsilon \zeta(x, \frac{x}{\varepsilon}, t )\\ &= \varphi(x, t) + \varepsilon \zeta_1(x, t)\, \zeta_2( \frac{x}{\varepsilon})\\ &=\varphi(x, t) + \varepsilon \zeta^\varepsilon(x, t), \end{aligned} \end{equation} where $\varphi \in {\EuScript D}(\Omega_T), \zeta_1 \in {\EuScript D}(\Omega_T), \zeta_2 \in C^\infty_{\rm per}(Y)$, and plug the function $\varphi_{\ell}$ in \eqref{eq:1.2.2}. This yields \begin{equation} \label{eq:1.7.1} \begin{aligned} &-\Phi_{f}\,\int_{\Omega_{T}} \chi_{f}^{\varepsilon}(x) \widetilde\Theta_{\ell,f}^{\varepsilon} \Big[ \frac{\partial \varphi}{\partial t} + \varepsilon \frac{\partial \zeta^\varepsilon}{\partial t} \Big] \, dx\, dt \\ &+ K_{f}\, \int_{\Omega_{T}} \chi_{f}^{\varepsilon}(x) \widehat \rho^{\varepsilon}_{\ell,f} \Big\{\lambda_{\ell}(\widetilde S^\varepsilon_{f}) \Big[\nabla \widetilde P^\varepsilon_{f} - \widehat \rho^{\varepsilon}_{\ell,f}\, \vec{g} \Big] + \nabla \beta(\widetilde S^\varepsilon_f) \Big\} \\ &\times [\nabla \varphi + \varepsilon \nabla_x \zeta^\varepsilon + \nabla_y \zeta^\varepsilon ]\, dx\, dt\\ &-\Phi_m\, \int_{\Omega_{T}} \chi_{f}^{\varepsilon}(x)\, \Theta^\varepsilon_{\ell,m} \Big[ \frac{\partial \varphi}{\partial t} + \varepsilon \frac{\partial \zeta^\varepsilon}{\partial t} \Big]\, dx\, dt\\ &+ \varepsilon^2\,K_m\, \int_{\Omega_{T}} \chi_{f}^{\varepsilon}(x)\,\widehat \rho^{\varepsilon}_{\ell,m} \Big\{\lambda_{\ell}(S^\varepsilon_m) (\nabla P^\varepsilon_m - \widehat \rho^{\varepsilon}_{\ell,m} \vec{g} ) + \nabla \beta(S^\varepsilon_m) \Big\}\\ &\times [\nabla \varphi + \varepsilon \nabla_x \zeta^\varepsilon + \nabla_y \zeta^\varepsilon ]\, dx\, dt = 0, \end{aligned} \end{equation} where \begin{equation}\label{rho-ell} \begin{gathered} \widehat \rho^{\varepsilon}_{\ell,f} := \rho_{\ell}(\widetilde S_{f}^{\varepsilon}, \widetilde P_{f}^{\varepsilon}) := \rho_{\ell} (\widetilde P_{f}^{\varepsilon} + G_{\ell}(\widetilde S_{f}^{\varepsilon})), \\ \widehat \rho^{\varepsilon}_{\ell,m} := \rho_{\ell}\,(S_m^{\varepsilon}, P_m^{\varepsilon}) := \rho_{\ell} ( P_m^{\varepsilon} + G_{\ell}(S_m^{\varepsilon})). \end{gathered} \end{equation} Taking into account the convergence results from Lemmata \ref{thm17}, \ref{cor:28}, we pass to the two-scale limit in \eqref{eq:1.7.1} as $\varepsilon \to 0$ and get the first homogenized equation \begin{equation} \label{eq:1.7.2} \begin{aligned} &-|Y_{f}| \Phi_{f} \int_{\Omega_{T}} \Theta^\star_{\ell} \frac{\partial\varphi}{\partial t}\,dx\, dt \\ &+ K_{f} \int_{\Omega_{T}\times Y_{f}} \rho_{\ell}\Big\{ \lambda_{\ell}(S) \big[\nabla P + \nabla_{y} P_1 - \rho_{\ell}\vec{g} \big] + \big[\nabla\beta + \nabla_{y}\beta_1\big] \Big\} \cdot \Upsilon_\varphi\,dy\,dx\,dt\\ & = \Phi_m\, \int_{\Omega_{T}\times Y_m} {\mathsf L}_{\ell,m}(x, y, t)\frac{\partial\varphi}{\partial t}(x, t)\,dy\,dx\,dt, \end{aligned} \end{equation} where \begin{equation} \label{rho-theta-h-ell} \rho_{\ell}:= \rho_{\ell}\,\big(P + G_{\ell}(S)\big) \quad \text{and} \quad \Theta^\star_{\ell} := S\, \rho_{\ell}; \end{equation} the function $\Upsilon_\varphi$ is $\Upsilon_\varphi := [\nabla \varphi + \zeta_1 \nabla_y \zeta_2]$, finally, the function ${\mathsf L}_{\ell,m} = {\mathsf L}_{\ell,m}(x, y, t)$ is defined as a two-scale limit of the sequence $\{\Theta^\varepsilon_{\ell,m}\}_{\varepsilon>0}$ as $\varepsilon \to 0$. Namely, \begin{equation} \label{unknown-2} \Theta^\varepsilon_{\ell,m} \stackrel {2s} \rightharpoonup {\mathsf L}_{\ell,m} \in L^2(\Omega_T; L^2_{\rm per}(Y_m)). \end{equation} In a similar way we pass to the two--scale limit in equation \eqref{eq:1.2.3} and get the second homogenized equation \begin{equation} \label{eq:1.7.3} \begin{aligned} &-|Y_{f}| \Phi_{f} \int_{\Omega_{T}} \Theta^\star_g \frac{\partial\varphi}{\partial t}\,dx\, dt \\ &+ K_{f} \int_{\Omega_{T}\times Y_{f}} \rho_g\Big\{ \lambda_g(S) \big[\nabla P + \nabla_{y} P_1 - \rho_g\vec{g} \big] - \big[\nabla\beta + \nabla_{y}\beta_1\big] \Big\} \cdot \Upsilon_\varphi\,dy\,dx\,dt \\ & = \Phi_m \int_{\Omega_{T}\times Y_m} {\mathsf L}_{g,m}(x, y, t) \frac{\partial\varphi}{\partial t}(x, t)\,dy\,dx\,dt, \end{aligned} \end{equation} where \begin{equation} \label{rho-theta-h-g} \rho_g:= \rho_g\,\big(P + G_g(S)\big) \quad \text{and} \quad \Theta^\star_g := (1 - S)\, \rho_g; \end{equation} the function ${\mathsf L}_{g,m} = {\mathsf L}_{g,m}(x, y, t)$ is defined as a two-scale limit of the sequence $\{\Theta^\varepsilon_{g,m}\}_{\varepsilon>0}$ as $\varepsilon \to 0$. Namely, \begin{equation} \label{unknown-3} \Theta^\varepsilon_{g,m} \stackrel {2s} \rightharpoonup {\mathsf L}_{g,m} \in L^2(\Omega_T; L^2_{\rm per}(Y_m)). \end{equation} \subsection{Identification of the corrector functions $P_1$, $\beta_1$ and homogenized equations} \label{ident-p-beta} \noindent\textbf{Step 1.} Identification of $P_1$. Consider the equations \eqref{eq:1.7.2}, \eqref{eq:1.7.3}. We set $\varphi \equiv 0$. We get: Taking into account that $\rho_{\ell}, \rho_g$ are strictly positive and do not depend on $y$, we obtain \begin{gather} \label{hernya-3} \int_{Y_{f}} \Big\{ \lambda_{\ell}(S) \big[\nabla P + \nabla_{y} P_1 - \rho_{\ell}\vec{g} \big] + \big[\nabla\beta + \nabla_{y}\beta_1\big] \Big\} \cdot \nabla \zeta_2(y)\,dy = 0, \\ \label{hernya-4} \int_{Y_{f}} \Big\{ \lambda_g(S) \big[\nabla P + \nabla_{y} P_1 - \rho_g\vec{g} \big] - \big[\nabla\beta + \nabla_{y}\beta_1\big] \Big\} \cdot \nabla \zeta_2(y)\,dy = 0. \end{gather} Now adding \eqref{hernya-3} and \eqref{hernya-4}, we obtain \begin{equation} \label{hernya-5} \int_{Y_{f}} \Big\{ \lambda(S)\big[\nabla P + \nabla_{y} P_1 \big] - \big[\lambda_{\ell}(S)\, \rho_{\ell} + \lambda_g(S)\, \rho_g \big]\, \vec{g} \Big\} \cdot \nabla_y \zeta_2(y)\,dy = 0. \end{equation} Taking into account condition (A4), from \eqref{hernya-5}, we obtain \begin{equation} \label{hernya-6} \int_{Y_{f}} \Big\{\nabla P - {\mathsf X}(S, P)\, \vec{g} + \nabla_{y} P_1 \Big\} \cdot \nabla_y \zeta_2(y)\,dy = 0, \end{equation} where \begin{equation} \label{hernya-7} {\mathsf X}(S, P) := \frac{\lambda_{\ell}(S)\, \rho_{\ell} + \lambda_g(S) \rho_g}{\lambda(S)}. \end{equation} Now we proceed in a standard way (see, e.g., \cite{hor}). Let $\xi_j$ be a solution of \begin{equation} \label{1newpoint-4} \begin{gathered} - \Delta_y \xi_j = 0 \quad \text{in } Y_{f}; \\ \nabla_y \xi_j \cdot \vec \nu_y = - \vec e_j \cdot \vec \nu_y \quad \text{on } \Gamma_{fm}\\ y \mapsto \xi_j(y)\quad Y\text{-periodic}. \end{gathered} \end{equation} Then the function $P_1$ can be represented as \begin{equation} \label{1bvp-2} P_1(x, y, t) = \sum^d_{j=1} \xi_j(y) \Big(\frac{\partial P}{\partial x_j}(x, t) - {\mathsf X}(S, P) g_j\Big). \end{equation} \noindent\textbf{Step 2.} Identification of $\beta_1$. Now we turn to the identification of the function $\beta_1$. From \eqref{hernya-3} and \eqref{hernya-6}, we obtain $$ \int_{Y_{f}} \Big\{ \lambda_{\ell}(S)\, {\mathsf X}(S, P)\, \vec{g} - \lambda_{\ell}(S)\, \rho_{\ell} \vec{g} + \big[\nabla\beta + \nabla_{y}\beta_1\big] \Big\} \cdot \nabla_y \zeta_2(y)\,dy = 0. $$ From this equation we have \begin{equation} \label{hernya-8} \int_{Y_{f}} \big[\nabla\beta + \nabla_{y}\beta_1\big] \cdot \nabla_y \zeta_2(y)\,dy = -\int_{Y_{f}} \Big\{ \lambda_{\ell}(S)\, \big[ {\mathsf X}(S, P) - \rho_{\ell} \big]\, \vec{g} \Big\} \cdot \nabla_y \zeta_2(y)\,dy. \end{equation} Now, from \eqref{hernya-8} we get the following equation for the function $\beta_1$: \begin{equation} \label{hernya-10} \int_{Y_{f}} \big[\nabla\beta - {\mathsf Z}(S, P)\, \vec{g} + \nabla_{y}\beta_1\big] \cdot \nabla_y \zeta_2(y)\,dy = 0, \end{equation} where \begin{equation} \label{hernya-9} {\mathsf Z}(S, P) := \frac{\lambda_{\ell}(S) \lambda_g(S)}{\lambda(S)} \big[\rho_{\ell}- \rho_g\big]. \end{equation} Then the function $\beta_1$ can be represented as \begin{equation} \label{1bvp-2-new} \beta_1(x, y, t) = \sum^d_{j=1} \xi_j(y) \Big(\frac{\partial\, \beta(S)}{\partial x_j}(x, t) - {\mathsf Z}(S, P)\, g_j\Big). \end{equation} \noindent\textbf{Step 3.} Homogenized equation for the $l$ phase. Now we are in position to obtain the first homogenized equation. Choosing $\zeta_2 = 0$, from \eqref{eq:1.7.2}, we obtain \begin{equation} \label{hernya-11} \begin{aligned} &\Phi_{f} |Y_{f}| \frac{\partial \Theta^\star_{\ell}}{\partial t} - \operatorname{div}_x \Big\{K_{f} \lambda_{\ell}(S) \rho_{\ell} \int_{Y_{f}} \big[\nabla P + \nabla_{y} P_1 - {\mathsf X}(S, P) \vec{g} \big]\, dy\\ &+ K_{f} \int_{Y_{f}} \lambda_{\ell}(S) \rho_{\ell}\big[{\mathsf X}(S, P) - \rho_{\ell}\big]\vec{g}\, dy \Big\} \\ &- \operatorname{div}_x \Big\{ K_{f} \int_{Y_{f}} \rho_\ell \big[\nabla\beta + \nabla_{y}\beta_1 - {\mathsf Z}(S, P)\, \vec{g}\big] \, dy + K_{f} \int_{Y_{f}} \rho_\ell {\mathsf Z}(S, P)\, \vec{g} \, dy \Big\} \\ & = - \Phi_m\, \int_{Y_m} \frac{\partial {\mathsf L}_{\ell,m}}{\partial t}(x, y, t)\,dy. \end{aligned} \end{equation} Then taking into account the definitions of the functions ${\mathsf X}(S, P)$ and ${\mathsf Z}(S, P)$, from \eqref{hernya-11}, we have \begin{equation} \begin{aligned} &\Phi_{f} |Y_{f}| \frac{\partial \Theta^\star_{\ell}}{\partial t} - \operatorname{div}_x \Big\{K_{f} \lambda_{\ell}(S) \rho_{\ell}\int_{Y_{f}} \big[\nabla P + \nabla_{y} P_1 - {\mathsf X}(S, P) \vec{g} \big]\, dy \Big\} \\ &- \operatorname{div}_x \Big\{ K_{f} \int_{Y_{f}} \rho_\ell \big[\nabla\beta + \nabla_{y}\beta_1 - {\mathsf Z}(S, P)\, \vec{g}\big] \, dy \Big\}\\ &= - \Phi_m\, \int_{Y_m} \frac{\partial {\mathsf L}_{\ell,m}}{\partial t}(x, y, t) \,dy. \end{aligned} \label{hernya-12} \end{equation} Now we use the representations for the corrector functions $P_1$ and $\beta_1$, to obtain \begin{equation} \begin{aligned} &\Phi^\star \frac{\partial \Theta^\star_{\ell}}{\partial t} - \operatorname{div}_x \Big\{\mathbb{K}^\star\,\lambda_{\ell}(S) \rho_\ell \big[\nabla P - {\mathsf X}(S, P)\, \vec{g}\big]\\ &+ \mathbb{K}^\star\,\rho_\ell \big[\nabla\beta(S) - {\mathsf Z}(S, P) \vec{g} \big] \Big\}\\ &= - \frac{\Phi_m}{|Y_m|}\, \int_{Y_m} \frac{\partial {\mathsf L}_{\ell,m}}{\partial t}(x, y, t)\,dy, \end{aligned} \label{hernya-13} \end{equation} where the effective porosity $\Phi^\star$ and $\mathbb{K}^\star$, the effective permeability tensor, are given by $$ \Phi^\star := \Phi_{f}\frac{|Y_{f}|}{|Y_m|} \quad \text{and} \quad \mathbb{K}^{\star}_{ij} := \frac{K_{f}}{|Y_m|}\, \int_{Y_\mathsf{f}}\, [\nabla_y \xi_i + \vec e_i ]\,\cdot [\nabla_y \xi_j + \vec e_j ]\, dy \quad (1 \leqslant i, j \leqslant d). $$ Since $\lambda_{\ell}(S) {\mathsf X}(S, P) + {\mathsf Z}(S, P) = \lambda_{\ell}(S) \rho_\ell$, from \eqref{hernya-13} we obtain \begin{equation} \begin{aligned} &\Phi^\star \frac{\partial \Theta^\star_{\ell}}{\partial t} - \operatorname{div}_x \Big\{\mathbb{K}^\star\,\rho_{\ell} \big[\lambda_{\ell}(S)\,\nabla P + \nabla\beta(S) - \lambda_{\ell}(S) \rho_{\ell} \vec{g}\big]\Big\} \\ &= - \frac{\Phi_m}{|Y_m|}\, \int_{Y_m} \frac{\partial {\mathsf L}_{\ell,m}}{\partial t}(x, y, t)\,dy. \end{aligned} \label{1newpoint-7:2ok} \end{equation} \noindent\textbf{Step 4.} Homogenized equation for the $g$ phase. In a similar way, choosing $\zeta_2 = 0$, from the equation \eqref{eq:1.7.3}, we derive the second homogenized equation \begin{equation} \begin{aligned} &\Phi^\star \frac{\partial \Theta^\star_g}{\partial t} - \operatorname{div}_x \Big\{\mathbb{K}^\star\,\rho_g \big[\lambda_g(S)\,\nabla P - \nabla\beta(S) - \lambda_g(S)\, \rho_g \vec{g}\big] \Big\} \\ &= - \frac{\Phi_m}{|Y_m|}\, \int_{Y_m} \frac{\partial {\mathsf L}_{g,m}}{\partial t}(x, y, t)\,dy. \end{aligned} \label{hernya-18ok} \end{equation} Let us introduce now the functions that is naturally to call the homogenized phase pressures. Namely, we define \begin{equation} \label{phase-pres} P_{\ell} := P + G_{\ell}(S) \quad \text{and} \quad P_g := P + G_g(S). \end{equation} Then from \eqref{1newpoint-7:2ok}--\eqref{phase-pres} we get the equations \begin{gather} \Phi^\star \frac{\partial \Xi^\star_{\ell}}{\partial t} - \operatorname{div}_x \Big\{\mathbb{K}^\star \rho_{\ell} \lambda_{\ell}(S) \big[ \nabla P_{\ell} - \rho_{\ell} \vec{g}\big] \Big\} = - \frac{\Phi_m}{|Y_m|} \int_{Y_m} \frac{\partial {\mathsf L}_{\ell,m}}{\partial t}(x, y, t)\,dy; \label{newpoint-9} \\ \Phi^\star \frac{\partial \Xi^\star_g}{\partial t} - \operatorname{div}_x \Big\{\mathbb{K}^\star \rho_g \lambda_g(S) \big[ \nabla P_g - \rho_g \vec{g}\big] \Big\} = - \frac{\Phi_m}{|Y_m|} \int_{Y_m} \frac{\partial {\mathsf L}_{g,m}}{\partial t}(x, y, t)\,dy, \label{newpoint-10} \end{gather} where $$ \Xi^\star_{\ell} := S \rho_{\ell}\,(P_{\ell}) \quad \text{and} \quad \Xi^\star_g := (1 - S) \rho_g\,(P_g). $$ At this stage, the homogenization process is finished by obtaining equations \eqref{newpoint-9}, \eqref{newpoint-10}. However, if in addition we want to express the functions ${\mathsf L}_{\ell,m}, {\mathsf L}_{g,m}$ in an explicit form, we need additional assumptions. Here, we assume that the cell problem~\eqref{hom-5} has a unique solution, then following the ideas of \cite{blm,choq}, we obtain as in \cite{ap-m2as} these limits explicitly. This completes the proof of Theorem \ref{t-hom-main}. \end{proof} \begin{remark} \rm We conclude this section with a remark about the obtained homogenized model. The scaling is such that, in the final homogenized equations, the less permeable part of the matrix contributes as a nonlinear memory term which we can represent explicitly in the case where the cell problem has a unique solution. A uniqueness proof is not ready to the moment, it is out of the scope of this paper, it remains an open problem even for the incompressible model, but its aspects and their applications must be considered for further achievements for this homogenization approach. \end{remark} \subsection*{Acknowledgments} This work was done in the framework of the International Research Group of CNRS GDRI LEM2I (Laboratoire Euro-Maghr\'ebin de Math\'emati\-ques et de leurs Interactions). Most of the work on this paper was done when L.~Ait Mahiout, A.~Mokrane and L.~Pankratov were visiting the Applied Mathematics Laboratory of the University of Pau \& CNRS. They are grateful for the invitations and the hospitality. The work of L. Pankratov was partially supported by the Russian Academic Excellence Project 5top100 and by the Russian Scientific Fund [grant number 15-11-00015]. \begin{thebibliography}{10} \bibitem{acdp} E. Acerbi, V. Chiad\`o Piat, G. Dal Maso, D. Percival; \emph{An extension theorem from connected sets, and homogenization in general periodic domains}, J. Nonlinear Analysis 18 (1992), 481--496. \bibitem{Allaire} G. Allaire; \emph{Homogenization and two-scale convergence}, SIAM J. Math. Anal. 23 (1992), 1482--1518. \bibitem{alt-Bene} H. W. Alt, E. di Benedetto; \emph{Nonsteady flow of water and oil through inhomogeneous porous media}, Ann. Scuola Norm. Sup. Pisa Cl. Sci. 12 (1985), 335--392. \bibitem{Siam} B. Amaziane, S. Antontsev, L. Pankratov, A. Piatnitski; \emph{Homogenization of immicible compressible two-phase flow in porous media: application to gas migration in a nulear repository}, SIAM J. Multiscale Model. Simul. 8 (5) (2010), 2023--2047. \bibitem{AJ} B. Amaziane, M. Jurak; \emph{A new formulation of immiscible compressible two-phase flow in porous media}, C. R. Mecanique 336 (2008), 600--605. \bibitem{AJV} B. Amaziane, M. Jurak, A. Vrba\v ski; \emph{Homogenization results for a coupled system modeling immiscible compressible two-phase flow in porous media by the concept of global pressure}, Appl. Anal. 92 (7) (2013), 1417--1433. \bibitem{AJZ1} B. Amaziane, M. Jurak, A. \v Zgalji\'c-Keko; \emph{Modeling and numerical simulations of immiscible compressible two-phase flow in porous media by the concept of global pressure}, Transp. Porous Media 84 (1) (2010), 133--152. \bibitem{AJZ2} B. Amaziane, M. Jurak, A. \v Zgalji\'c-Keko; \emph{An existence result for a coupled system modeling a fully equivalent global pressure formulation for immiscible compressible two-phase flow in porous media}, J. Differential Equations 250 (2011), 1685--1718. \bibitem{ap-m2as} B. Amaziane, L. Pankratov; \emph{Homogenization of a model for water-gas flow through double-porosity media}, Math. Meth. Appl. Sci. 39 (2016), 425--451. \bibitem{app-1} B. Amaziane, L. Pankratov, A. Piatnitski; \emph{The existence of weak solutions to immiscible compressible two-phase flow in porous media: the case of fields with different rock-types}, Discrete Continuous and Dynamical Systems, Ser. B 15 (2013), 1217--1251. \bibitem{app-2} B. Amaziane, L. Pankratov, A. Piatnitski; \emph{Homogenization of immiscible compressible two-phase flow in highly heterogeneous porous media with discontinuous capillary pressures}, Math. Models Methods Appl. Sci. 24 (7) (2014), 1421--1451. \bibitem{APR} B. Amaziane, L. Pankratov, V. Rybalko; \emph{On the homogenization of some double-porosity models with periodic thin structures}, Appl. Anal. 88 (10-11) (2009), 1469--1492. \bibitem{ant-kaz-mon1990} S.~N. Antontsev, A.~V. Kazhikhov, V.~N. Monakhov; \emph{Boundary value problems in mechanics of nonhomogeneous fluids}. North-Holland Publishing Co., Amsterdam, 1990. \bibitem{Arbogast-92} T. J. Arbogast; \emph{The existence of weak solutions to single porosity and simple dual-porosity models of two-phase incompressible flow}, Nonlinear Anal. 19 (1992), 1009--1031. \bibitem{adh} T. Arbogast, J. Douglas, U. Hornung; \emph{Derivation of the double porosity model of single phase flow via homogenization theory}, SIAM J. Appl. Math. Anal. 21 (1990), 823--836. \bibitem{bar} G. I. Barenblatt, Y. Zheltov, I. N. Kochina; \emph{Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks}, Journal of Applied Mathematics and Mechanics 24 (1960), 1286--1303. \bibitem{bear} J. Bear, C.F. Tsang, G. de Marsily; \emph{Flow and contaminant transport in fractured rock}. Academic Press Inc, London, 1993. \bibitem{Bourgeat-Gip2004} A. Bourgeat, O. Gipouloux, E. Maru\v si\'c-Paloka; \emph{Mathematical modeling of an under-ground waste disposal site by upscaling}, Math. Methods Appl. Sci. 27 (2004), 381--403. \bibitem{Bourgeat-Gip2010} A. Bourgeat, O. Gipouloux, F. Smai; \emph{Scaling up of source terms with random behavior for modelling transport migration of contaminants in aquifers}, Nonlinear Anal. Real World Appl. 11 (2010), 4513--4523. \bibitem{Bourgeat-Marusic2005} A. Bourgeat, E. Maru\v si\'c-Paloka; \emph{A homogenized model of an underground waste repository including a disturbed zone}, SIAM J. Multiscale Model. Simul. 3 (2005), 918--939. \bibitem{Bourgeat-Marusic2010} A. Bourgeat, E. Maru\v si\'c-Paloka, A. Piatnitski; \emph{Scaling up of an underground nuclear waste repository including a possibly damaged zone}, Asymptot. Anal. 67 (2010), 147--165. \bibitem{Bourgeat-Piat2010} A. Bourgeat, A. Piatnitski; \emph{Averaging of a singular random source term in a diffusion convection equation}, SIAM J. Math. Anal. 42 (2010), 2626--2651. \bibitem{blm} A. Bourgeat, S. Luckhaus, A. Mikelic; \emph{Convergence of the homogenization process for a double-porosity model of immiscible two-phase flow}, SIAM J. Appl. Math. Anal. 27 (1996), 1520--1543. \bibitem{Brezis} H. Brezis; \emph{Functional analysis, Sobolev spaces and partial differential equations}. Masson, Paris, 1983. \bibitem{CC-PM} C. Canc\`es, P. Michel; \emph{An existence result for multidimensional immiscible two-phase flows with discontinuous capillary pressure field}, SIAM J. Math. Anal. 44 (2012), 966--992. \bibitem{caro} F. Caro, B. Saad, M. Saad; \emph{Study of degenerate parabolic system modelling the hydrogen displacement in a nuclear waste repository}, Discrete and Continuous Dynamical Systems, Ser. S 7 (2014), 191--205. \bibitem{GC-JJ} G. Chavent, J. Jaffr\'e; \emph{Mathematical models and finite elements for reservoir simulation}. North-Holland, Amsterdam, 1986. \bibitem{Chen-I} Z. Chen; \emph{Degenerate two-phase incompressible flow. I. Existence, uniqueness and regularity of a weak solution}, J. Differential Equations 171 (2001), 203--232. \bibitem{Chen-II} Z. Chen; \emph{Degenerate two-phase incompressible flow. II. Regularity, stability and stabilization}, J. Differential Equations 186 (2002), 345--376. \bibitem{Chen} Z. Chen; \emph{Homogenization and simulation for compositional flow in naturally fractured reservoirs}, J. Math. Anal. Appl. 326 (1) (2007), 12--32. \bibitem{Chen-Huan-Ma} Z. Chen, G. Huan, Y. Ma; \emph{Computational methods for multiphase flows in porous media}, SIAM, Philadelphia, 2006. \bibitem{choq} C. Choquet; \emph{Derivation of the double porosity model of a compressible miscible displacement in naturally fractured reservoirs}, Appl. Anal. 83 (2004) 477--500. \bibitem{ddg} D. Cioranescu, A. Damlamian, G. Griso; \emph{Periodic unfolding and homogenization}, C. R. Acad. Sci. Paris, Ser. I 335 (2002), 99--104. \bibitem{GWC-RES}G. W. Clark, R. E. Showalter; \emph{Two-scale convergence of a model for flow in a partially fissured medium}, Electron. J. Diff. Eqns. 2 (1999), 1--20. \bibitem{Gallouet} R. Eymard, T. Gallouet, R. Herbin; \emph{Finite volume methods. Handbook of Numerical Analysis}, P.G. Ciarlet, J.L. Lions eds., 7 (2000), 713--1020. \bibitem{GG-MM-1996} G. Gagneux, M. Madaune-Tort; \emph{Analyse math\'ematique de mod\`eles non-lin\'eaires de l'ing\'e\-nierie p\'etroli\`ere}. Springer-Verlag, Berlin, 1996. \bibitem{CG-MS1} C. Galusinski, M. Saad; \emph{On a degenerate parabolic system for compressible, immiscible, two-phase flows in porous media}, Adv. Differential Equations 9 (2004), 1235--1278. \bibitem{CG-MS2} C. Galusinski, M. Saad; \emph{Water-gas flow in porous media}, Discrete Contin. Dyn. Syst., Ser. B, {\bf 9} (2008), 281-308. \bibitem{CG-MS3} C. Galusinski, M. Saad; \emph{Two compressible immiscible fluids in porous media}, J. Differential Equations 244 (2008) 1741--1783. \bibitem{CG-MS4} C. Galusinski, M. Saad; \emph{Weak solutions for immiscible compressible multifluid flows in porous media}, C. R. Acad. Sci. Paris, S\'er. I 347 (2009), 249--254. \bibitem{Gipouloux-Smai} O. Gipouloux, F. Smai; \emph{Scaling up of an underground waste disposal model with random source terms}, Internat. J. Multiscale Comput. Engin. 6 (2008), 309--325. \bibitem{Gloria-Etal} A. Gloria, T. Goudon, S. Krell; \emph{Numerical homogenization of a nonlinearly coupled elliptic-parabolic system}, reduced basis method, and application to nuclear waste storage, Math. Models Methods Appl. Sci. 23 (2013), 2523--2560. \bibitem{HR} R. Helmig; \emph{Multiphase flow and transport processes in the subsurface}. Springer, Berlin, 1997. \bibitem{HOS-2012} P. Henning, M. Ohlberger, B. Schweizer; \emph{Homogenization of the degenerate two-phase flow equations}, Math. Models Methods Appl. Sci. 23 (2013), 2323--2352. \bibitem{hor} U. Hornung; \emph{Homogenization and porous media}. Springer-Verlag, New York, 1997. \bibitem{khal-saad1} Z. Khalil, M. Saad; \emph{Solutions to a model for compressible immiscible two phase flow in porous media}, Electronic Journal of Differential Equations 122 (2010), 1--33. \bibitem{khal-saad2} Z. Khalil, M. Saad; \emph{On a fully nonlinear degenerate parabolic system modeling immiscible gas-water displacement in porous media}, Nonlinear Analysis: Real World Applications 12 (2011) 1591--1615. \bibitem{Kro-Luck} D. Kroener, S. Luckhaus; \emph{Flow of oil and water in a porous medium}, J. Differential Equations 55 (1984), 276--288. \bibitem{panf}M. Panfilov; \emph{Macroscale models of flow through highly heterogeneous porous media}. Kluwer Academic Publishers, London, 2000. \bibitem{Shaw} R. P. Shaw; \emph{Gas generation and migration in deep geological radioactive waste repositories}. Geological Society, 2015. \bibitem{van} T. D. Van Golf-Racht; \emph{Fundamentals of fractured reservoir engineering}. Elsevier Scientific Pulishing Company, Amsterdam, 1982. \bibitem{Vazquez} J. L. V\'azquez; \emph{The porous medium equation. Mathematical theory}. Oxford University Press, Oxford, 2007. \bibitem{yeh2} L. M. Yeh; \emph{Homogenization of two-phase flow in fractured media}, Math. Methods Appl. Sci. 16 (2006), 1627--1651. \end{thebibliography} \end{document}