Mathematical modeling of cyclic treatments of chronic myeloid leukemia

1.   Introduction.

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$.

In this paper, we focus our attention on the modeling of immiscible compressible two-phase flow through heterogeneous reservoirs in the framework of the geological disposal of radioactive waste. The long-term safety of the disposal of nuclear waste is an important issue in all countries with a significant nuclear program. One of the solutions envisaged for managing waste produced by nuclear industry is to dispose the radioactive waste in deep geological formations chosen for their ability to delay and to attenuate possible releases of radionuclides in the biosphere. Repositories for the disposal of high-level and long-lived radioactive waste generally rely on a multi-barrier system to isolate the waste from the biosphere. The multibarrier system typically comprises the natural geological barrier provided by the repository host rock and its surroundings and an engineered barrier system, i.e. engineered materials placed within a repository, including the waste form, waste canisters, buffer materials, backfill and seals, for more details see for instance [43]. An important task of the safety assessment process is the handling of heterogeneities of the geological formation.

In the frame of designing nuclear waste geological repositories, a problem of possible two-phase flow of water and gas, mainly hydrogen, appears, for more details see for instance [43]. Multiple recent studies have established that in such installrations important amounts of gases are expected to be produced in particular due to the corrosion of metallic components used in the repository design, see e.g. [27,42] and the references therein. The French Agency for the Management of Radioactive Waste (Andra) [11] is currently investigating the feasibility of deep geological disposal of radioactive waste in an argillaceous formation. A question related to the long-term performance of the repository concerns the impact of the hydrogen gas generated in the wastes on the pressure and saturation fields in the repository and the host rock.

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 [14,23,26,29,36,38,44] 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 been shown under various assumptions on physical data; see for instance [3,14,15,21,24,25,23,29,41] and the references therein. A recent review of the mathematical homogenization methods developed for incompressible immiscible two-phase flow in porous media and compressible miscible flow in porous media can be viewed in [4,37,38]. We refer for instance to [16,17,18,19,20,34,35] 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.

However, as reported in [9], 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 [30,31,32]. Namely, in [30] certain terms related to the compressibility are neglected, and in [31,32] the mass densities are assumed not to depend on the physical pressure, but on Chavent's global pressure. In the articles [22,33,39,40], a more general immiscible compressible two-phase flow model in porous media is considered for fields with a single rock type and [9] treated the case with several types of rocks. In [4,10] 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 new global pressure concept was introduced in [5,7] for modeling immiscible, compressible two-phase flow in porous media without any simplifying assumptions. The resulting equations are written in a fractional flow formulation and lead to a coupled system which consists of a nonlinear parabolic equation (the global pressure equation) and a nonlinear diffusion-convection one (the saturation equation). This new formulation is fully equivalent to the original phase equations formulation, i.e. where the phase pressures and the phase saturations are primary unknowns. For this model, an existence result is obtained in [8] and homogenization results in [6].

Let us note that all the aforementioned works are restricted to the case where the gas density is bounded from below and above, contrarily to the present work. This assumption is too restrictive for some realistic problems, such as gas migration through engineered and geological barriers for a deep repository for radioactive waste. In this case the gas obeys the ideal gas law, i.e. the equation of state is given by $\varrho_g(p) \overset{def}{\mathop =} \sigma p$ where $ \varrho_g$ is the gas density, $ p_g$ is the gas pressure and $\sigma$ is a given constant. Then a new degeneracy appears in the evolution term of the gas pressure equation. In this paper we extend our previous results obtained in [4] to the more complex case of an ideal gas 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. To obtain these results we elaborated a new approach based on the ideas from [4] and regularization.

The rest of the paper is organized as follows. In Section 2 we describe the physical model and formulate the corresponding mathematical problem. We also provide the assumptions on the data.

The goal of Section 3 is to prove the existence result for the corresponding system of equations. The proof is divided into a number of steps. In subsection 3.1 we consider an auxiliary $\delta$-problem where the gas density admits a positive lower bound. The existence result for the $\boldsymbol\delta$-problem is given in subsection 3.2. In subsection 3.3 we obtain a number of a priori estimates for a solution of $\delta$-problem. Then in subsection 3.4 we prove a compactness result adapted to our model. Finally, in subsection 3.5 we pass to the limit in $\delta$-problems and complete the proof of the existence.

Section 4 is devoted to the corresponding homogenization problem. In subsection 4.1 we introduce the model with a periodic microstructure. We assume that both porosity and absolute permeability tensor are periodic rapidly oscillating functions.

Then subsection 4.2 we formulate the homogenization result. This result is proved in subsection4.3. The proof is based on the two-scale convergence technique. Our analysis relies essentially on a compactness result [4] which is rather involved due to the degeneracy and the nonlinearity of the system.

The last section is followed by some concluding remarks.

2.   Formulation of the problem.

We consider an immiscible compressible two-phase flow process in a porous reservoir $\Omega \subset \mathbb{R}^d$ ($d = 1, 2, 3$) which is a bounded Lipschitz domain. The time interval of interest is $(0,T)$. We focus here on the particular case of water and gas phases, but the consideration below is also valid for a general wetting phase and a non-wetting phase. Let $\Phi = \Phi(x)$ be the porosity of $\Omega$; $K = K(x)$ be the absolute permeability tensor of $\Omega$; $\varrho_w$, $\varrho_g$ are the densities of water and gas, respectively; $S_{w} = S_{w}(x, t)$, $S_{g} = S_{g}(x, t)$ are the saturations of water and gas in $\Omega\times (0,T)$; $k_{r,w} = k_{r,w}(S_{w})$, $k_{r,g} = k_{r,g}(S_{g})$ are the relative permeabilities of water and gas; $p_{w} = p_{w}(x,t)$, $p_{g} = p_{g}(x,t)$ are the pressures of water and gas in $\Omega\times (0,T)$.

In what follows, for the sake of presentation simplicity we neglect the source terms. Then the conservation of mass of each phase is described by (see, e.g., [23,26,36]):

where $\Omega_T \overset{def}{\mathop =} \Omega \times (0, T)$ with $T > 0$, the velocities of water and gas $\vec q_{w}$, $\vec q_{g}$ are defined by Darcy-Muskat's law:

Here $\vec g$, $\mu_w, \mu_g$ are the gravity vector and the viscosities of the water and gas, respectively.

From now on we assume that the density of the water is constant, which for the sake of simplicity will be taken equal to one, i.e. $ {\varrho_w(p_{w})} = {\rm Const}=1, $ and the gas density $\varrho_g$ obeys the ideal gas law and is given by the following function:

Here $\sigma$, $\varrho_{\rm max}$, and $p_{\rm max}$ are positive constants. Note that $\varrho_g$ is a continuous monotone and non-negative function.

The model is completed as follows. By the definition of saturations, one has

We set:

Then the curvature of the contact surface between the two fluids links the jump of pressure of two phases to the saturation by the capillary pressure law:

where $P^\prime_{c}(s)$ denotes the derivative of the function $P_{c}(s)$.

Now due to (6) and the assumption on the water density, we rewrite the system (1) as follows:

where $\lambda_g(S) := \lambda_g(1 - S)$ and

The system (8) have to be completed by appropriate boundary and initial conditions.

Boundary conditions.: We suppose that the boundary $\partial \Omega$ consists of two parts $\Gamma_{\rm inj}$ and $\Gamma_{\rm imp}$ such that $\Gamma_{\rm inj} \cap \Gamma_{\rm imp} = \emptyset$, $\partial \Omega = \overline\Gamma_{\rm inj} \cup \overline\Gamma_{\rm imp}$. The boundary conditions are given by:

where the velocities $\vec q_{w}, \vec q_{g}$ are defined in (2), (3).

Initial conditions.: The initial conditions read:

Notice that from (10) and (7) it follows that $S=1$ on $\Gamma_{\rm inj} \times (0,T)$. The initial condition for $S$ is uniquely defined by the equation

Then according to (9) the initial condition for $\Theta$ reads

Remark 1. It is important to underline that in the earlier works (see, e.g., [4,9,10,30,31,32,33]) it was assumed that the gas density admits a strictly positive lower bound:

2.1.   A fractional flow formulation.

In the sequel, we use a formulation obtained after transformation using the concept of the so called global pressure. In the case of incompressible two-phase flow this concept was introduced for the first time in [12,13]. Following [14,23], see also [26], we first recall the definition of the global pressure. It plays a crucial role, in particular, for compactness results. The idea of introducing the global pressure is as follows. We want to replace the water-gas flow by a flow of a fictive fluid obeying the Darcy law with a non-degenerating coefficient. Namely, we are looking for a pressure ${\mathsf P}$ and the coefficient $\gamma(S)$ such that $\gamma(S) > 0$ for all $S \in [0, 1]$, and

Then the global pressure, ${\mathsf P}$, is defined by:

the functions ${\mathsf G}_{w}(s)$ and ${\mathsf G}_{g}(s)$ will be introduced later on, in (19), (20). Now it is easy to see that

where

We set:

Then $\gamma(S) = \lambda(S)$. By construction, $\lambda(S) > 0$ for all $S \in [0, 1]$ (see the condition (A.5) below). Thus the relation (15) is established. Now we specify the functions ${\mathsf G}_{w}$, ${\mathsf G}_{g}$. We define ${\mathsf G}_{g}$ as follows:

The functions ${\mathsf G}_{w}$ are then defined by

Notice that from (19), (20) we get:

where

Now we link the capillary pressure and the mobilities. In a standard way (see, e.g., [40] or [9] for more details) we obtain the following identity:

where

Ii is also convenient to introduce the following function $\beta$:

where the function $\alpha$ is defined in (22). Notice that by the definition of the global pressure, (24), (25), and by the boundedness of $\lambda_{w}, \lambda_{g}$ (see the condition (A.5) below) the following relations holds:

In order to complete this section, let us calculate the value of the global pressure function ${\mathsf P}$ on $\Gamma_{\rm inj}$. In what follows (see condition (A.4) below) we assume that the capillary pressure function $P_{c}$ satisfies the condition $P_{c}(1) = 0$. Then from (10)$_1$ we have that the saturation $S$ equals one on $\Gamma_{\rm inj}$. Now the definition of the global pressure (16) implies that the function ${\mathsf P}$ on $\Gamma_{\rm inj}$ is a constant which we denote by ${\mathsf P}^1$.

2.2.   Main assumptions.

The main assumptions on the data are as follows:

(A.1) The porosity $\Phi \in L^{\infty}(\Omega)$, and there are positive constants $\phi_-, \phi^+$ such that $0 < \phi_- < \phi^+$ and

(A.2) The tensor $K \in (L^{\infty}(\Omega))^{d\times d}$ and there exist constants $K_-, K^+$ such that $0 < K_- < K^+$ and

(A.3) The function $\varrho_{g} = \varrho_g(p)$ is given by (4).

(A.4) The capillary pressure function $P_{c}(s) \in C^1([0, 1]; \mathbb{R}^+)$. Moreover, $P_{c}^\prime(s) < 0$ in $[0, 1]$ and $P_{c}(1) = 0$.

(A.5) The functions $\lambda_{w}, \lambda_{g}$ belong to the space $C([0, 1]; \mathbb{R}^+)$ and satisfy the following properties:

(ⅰ) $0 \leqslant \lambda_{w}, \lambda_{g} \leqslant 1$ in $[0, 1]$;

(ⅱ) $\lambda_{w}(0) = 0$ and $\lambda_{g}(1) = 0$;

(ⅲ) there is a positive constant $L_0$ such that $\lambda(s) = \lambda_{w}(s) + \lambda_{g}(s) \geqslant L_0 > 0$ in $[0, 1]$.

(A.6) The function $\alpha$ defined in (22) is an element of $C^1([0, 1]; \mathbb{R}^+)$. Moreover, $\alpha (0) = \alpha (1) = 0$, and $\alpha > 0$ in $(0, 1)$.

(A.7) The function $\beta^{-1}$, inverse of $\beta$ defined in (25) is a Hölder continuous function of order $\theta$ with $\theta \in (0, 1)$ on the interval $[0, \beta(1)]$. There exists a positive constant $C_\beta$ such that for all $s_1, s_2 \in [0, \beta(1)]$ the following inequality holds:

(A.8) The initial data for the pressures are such that $p^\bf{0}_{g}, p^\bf{0}_{w} \in L^2(\Omega)$.

(A.9) The initial data for the saturation is such that $S^\bf{0} \in L^\infty(\Omega)$ and $0 \leqslant S^\bf{0} \leqslant 1$ a.e. in $\Omega$.

The assumptions (A.1)-(A.9) are classical for two-phase flow in porous media.

3.   Existence result.

In order to define a weak solution of the above problem, we introduce the following Sobolev space:

The space $H^1_{\Gamma_{\rm inj}}(\Omega)$ is a Hilbert space. The norm in this space is given by

Theorem 3.1. Let assumptions (A.1)-(A.9) be fulfilled. Then there exist functions $\langle p_{g}, p_{w} , S\rangle$ such that:

(Ⅰ) 

where the function $\Theta$ is defined in (9); $S=1$ on $\Gamma_{\rm inj}$.

(Ⅱ) the maximum principle holds:

(Ⅲ) For any $\varphi_w, \varphi_g \in C^1([0, T]; H^1(\Omega))$ satisfying $\varphi_w = \varphi_g = 0$ on $\Gamma_{\rm inj} \times (0,T)$ and $\varphi_w(x,T) = \varphi_g(x,T) = 0$, we have:

with $\Theta$ defined in (9), and $P_{c}(S) = p_{g} - p_{w}$.

(Ⅳ) The initial conditions are satisfied in a weak sense as follows:

Furthermore, we have

and

with $S^\bf{0}$ and $\Theta^\bf{0}$ defined in (12) and (13), respectively.

The proof of Theorem 3.1 is divided into a several steps. It is based on a auxiliary existence result for the system obtained by approximation of the initial degenerate gas density $\varrho_g$ by a family of functions $\{\varrho^\delta_g\}_{\delta>0}$ that admit a positive lower bound. For such kind of system the desired existence result is proved in [9,30,31,32,33]. This result is formulated in subsection 3.1. Using the weak formulation of the regularized problem and the uniform in $\delta$ estimates for its solution, we, finally, prove Theorem 3.1.

3.1.   Auxiliary $\delta$-problem.

In this subsection we approximate the function $\varrho_g$ by a family of functions $\{\varrho^\delta_g\}_{\delta>0}$ that admit a positive lower bound. For each $\delta > 0$ we set:

Here $\sigma, \varrho_{\rm max}, p_{\rm max}$ are positive constants.

In addition to (8), consider the following family of problems:

where

System (40) have to be completed with the corresponding boundary and initial conditions.

Boundary conditions.: The boundary conditions read

where the velocities $\vec q^{\,\delta}_{w}, \vec q^{\,\delta}_{g}$ are given by:

Initial conditions.: The initial conditions read:

The remaining part of the Section is organized as follows. First, in subection 3.2 we recall the existence result for the system (40). Then we obtain the uniform in $\delta$ estimates for the solution of $\delta$-problem (40). In subection 3.4 we formulate the compactness and convergence results which we use in the proof of Theorem 3.1.

3.2.   An existence result for the $\boldsymbol\delta$-problem.

The goal of this subsection is to recall the existence result for the $\delta$-problem (40). First, we reformulate the condition (A.3) from subsection 2.2 in order to adapt it to our $\delta$-problem. For this problem the condition (A.3) becomes:

(A.3$_{\boldsymbol\delta}$) The function $\varrho^\delta_{g} = \varrho^\delta_g(p)$ is given by (39).

Now we are in position to formulate the existence result to $\delta$-problem (40). It reads.

Theorem 3.2. (see [9,33]) Let assumptions (A.1)-(A.2), (A.3$_{\boldsymbol\delta}$), (A.4)-(A.9) be fulfilled. Then, for each $\delta > 0$, there exist $\langle p^\delta_{g}, p^\delta_{w} , S^\delta\rangle$ such that:

(Ⅰ) 

where the function $\Theta^\delta$ is given by (41).

(Ⅱ) the maximum principle holds:

(Ⅲ) For any $\varphi_w, \varphi_g \in C^1([0, T]; H^1(\Omega))$ satisfying $\varphi_w = \varphi_g = 0$ on $\Gamma_{\rm inj} \times (0,T)$ and $\varphi_w(x,T) = \varphi_g(x,T) = 0$, the following integral identity holds:

Here $\Theta^\delta(x, 0) = \varrho^\delta_g(p_{g}^\bf{0})(1 - S^\bf{0})$ with the function $S^\bf{0}$ defined in condition

(A.9), and $P_{c}(S^\delta) = p^\delta_{g} - p^\delta_{w}$.

(Ⅳ) The initial conditions are satisfied in a weak sense as follows:

Furthermore, we have

and

3.3.   A priori estimates for a solution of $\delta$-problem (40).

We start this subsection by obtaining the energy equality for $\delta$-problem (40). The following result holds.

Lemma 3.3 (Energy equality for $\delta$-problem Let $\langle p^\delta_{w}, p^\delta_{g}, S^\delta\rangle$ be a solution to (40). Then

in the sense of distributions. Here

where

Notice that in the previous works (see, e.g., [4,9,10,30,31,32,33]), the function ${\mathsf R}_{g}$ was defined by ${\mathsf R}_{g}(p_g) \overset{def}{\mathop =} \int\limits_{0}^{p_g} \frac{d\xi}{\varrho_g(\xi)}$. However, in our case, with such a definition we do not have uniform in $\delta$ lower bound for the function ${\mathcal{R}}^\delta$. Thus, we have to modify the definition of ${\mathsf R}_{g}$, subtracting an appropriate constant $C = C(\delta)$. The properties of the functions ${\mathsf R}^\delta_{g}, {\mathcal{R}}^\delta$, and ${\mathcal{E}}^\delta$ are given in:

Lemma 3.4. Let ${\mathsf R}^\delta_{g}, {\mathcal{R}}^\delta$, and ${\mathcal{E}}^\delta$ be the functions defined by (55), (56). Then

(ⅰ) The function ${\mathcal{R}}^\delta$ is negative and bounded from below, that is

(ⅱ) The function ${\mathcal{E}}^\delta$ is bounded from below. Namely,

Proof of Lemma 3.4. Using the definition of the gas density $\varrho^\delta_g$ given by (39), it is easy to calculate that

Consider now the function ${\mathcal{R}}^\delta$. Due to (39) and (59), we have:

The last formula immediately implies (57). Now (58) follows easily from (57) and the estimate:

This completes the proof of Lemma 3.4.

In order to formulate a priori estimates for the solution to $\delta$-problem (40), we remark first that the global pressure ${\mathsf P}^\delta$ for the problem under consideration can be introduced in a way similar to one used in subsection 2.1 above. Then the desired a priori estimates for the solution of the $\delta$-problem can be easily derived from Lemmata 3.3, 3.4, and the equality (see subsection 2.1 for more details):

where the function $\mathfrak{b}(s)$ is defined in (24).

The following result holds.

Lemma 3.5. Let $\langle p^\delta_{w}, p^\delta_{g}, S^\delta \rangle$ be a solution to (40), the global pressure ${\mathsf P}^\delta$ is defined in (16), and the function $\beta(s)$ is defined in (25). Then

Here $C$ does not depend on $\delta$.

Proof of Lemma 3.5. Integrating (54) over the interval $(0, T)$, we get:

Let us estimate now the right-hand side of (66) from above. Due to the definition of the function ${\mathcal{E}}^\delta$, (55), and the initial conditions (44) we have that

where $S^\bf{0} = S^\bf{0}(x)$ is the initial condition of the saturation function (see condition (A.9) in Section 2.2). Now from condition (A.1) and the maximum principle (48), we easily obtain that

where $|\Omega|$ stands for the measure of the domain $\Omega$. Therefore,

where $\mathbb{C}_0$ is a constant which only depends on $\max\limits_{S \in [0, 1]} P_{c}(S)$, and the constant $\phi^+$. Combining condition (A.1), (66), (69), and the bound (61) yields:

Applying the Cauchy inequality, from (70), we deduce (63), and consequently (64).

The uniform estimates (65) can be obtained in the standard way from (40) with the help of (63). Lemma 3.5 is proved.

3.4.   Compactness and convergence results.

In this subsection we recall two compactness results that were obtained in [9].

Lemma 3.6 (Compactness lemma). Let $\Phi \in L^\infty(\Omega)$, and assume that there exist positive constants $\phi_1, \phi_2$ such that $0 < \phi_1 \leqslant \Phi(x) \leqslant \phi_2 < 1$ a.e. in $\Omega$. Assume moreover that a family $\{v^\delta \}_{\delta>0} \subset L^2(\Omega_T)$ satisfies the following properties:

1. the functions $v^\delta$ satisfy the ineqiality $0 \leqslant v^\delta \leqslant C$;

2. there exists a function $\varpi$ such that $\varpi(\xi) \rightarrow 0$ as $\xi\to 0$, and the following inequality holds true:

3. the estimate holds $\Vert \partial_t(\Phi\, v^\delta) \Vert_{L^2(0,T;H^{-1}(\Omega))} \leqslant C$.

Then the family $\{v^\delta \}_{\delta>0}$ is a compact set in $L^2(\Omega_T)$.

This result is a particular case of Lemma 4.2 proved in [4]. We apply the statement of Lemma 3.6 in order to prove the compactness of the sequences $\{\Theta^\delta\}_{\delta>0}$, $\{S^\delta\}_{\delta>0}$. As in [9] we obtain.

Proposition 1. Let $\{\Theta^\delta\}_{\delta>0} \subset L^2(\Omega_{T})$ be defined by (41). Then $\{\Theta^\delta\}_{\delta>0}$ is a compact set in the space $L^q(\Omega_{T})$ for all $q \in [1, +\infty)$.

Proposition 2. Let $\{S^\delta\}_{\delta>0} \subset L^2(\Omega_{T})$. Then, for all $q \in [1, +\infty)$, $\{S^\delta\}_{\delta>0}$ is a compact set in the space $L^q(\Omega_{T})$ for all $q \in [1, +\infty)$.

Now from Lemma 3.5 and Propositions 1, 2 we have.

Lemma 3.7. Up to a subsequence,

where $\Theta \overset{def}{\mathop =} \varrho_{g}({\mathsf P} + {\mathsf G}_{g}(S))\left(1 - S\right) \overset{def}{\mathop =} \varrho_g(p_{g})\,(1 - S)$;

for any $\psi \in C([0,1])$ such that $\psi(0) = 0$.

Proof of Lemma 3.7. The convergence (72) follows immediately from Proposition 2 and the limit function $S$ evidently satisfies (73). The relations (74) and (75) are the consequence of (64) from Lemma 3.5. The convergence (76) follows from (72)-(74), the inequality $0 \leqslant S^\delta \leqslant 1$, and the fact that $\varrho^\delta_g$ is monotone. In order to justify (76) we first observe that, due to the definitions of the functions $\varrho^\delta$ and $\varrho$,

Then for any $v \in L^\infty(\Omega_T)$, we have:

Denoting $\bar\Theta$ the limit of $\Theta^\delta$ and passing to the limit, as $\delta\to0$, in the last inequality, we obtain

Choosing $v = {\mathsf P} + \varkappa v_1$ and sending $\varkappa$ to zero yields

for any $v_1\in L^2(\Omega_T)$. This implies (76).

Finally, the convergence (77) can be proved by arguments similar to those from Lemma 4.2 in [40]. Lemma 3.7 is proved.

3.5.   Proof of Theorem 3.1.

We begin this subsection by studying the regularity properties of solution to (8).

3.5.1.   Regularity properties of a solution to system (8).

Taking into account the lower semi-continuity of the norm, by Lemma 3.5, we obtain:

Now we set:

We also recall the relation (23):

Then, taking into account (78), (80), and the last relation we obtain that the functions $p_{w}, p_{g}$ defined in (81) are such that

Thus properties (30)-(31) are established. The maximum principle (33) follows immediately from (48) and (72). Finally, the interpretation of the initial conditions can be done as in [40] (see also [9]).

3.5.2.   Passage to the limit in equations (49), (50).

Consider the equation (49), with $\varphi_w \in C^1([0, T]; H^1(\Omega))$ and is such that $\varphi_w = 0$ on $\Gamma_{\rm inj} \times (0,T)$ and $\varphi_w(x,T) = 0$.

Taking into account (72), one easily gets:

We then recall that

Then the third term on the left-hand side of (49) takes the form:

Now taking into account the convergence (72), (74), and (75), we obtain that

Returning now to the water pressure function $p_{w}$, we finally get:

Considering (72), one can check that the fourth term of (49) satisfies the relation

Thus, the saturation equation (34) is derived.

We turn to (50) with $\varphi_g \in C^1([0, T]; H^1(\Omega))$, $\varphi_g = 0$ on $\Gamma_{\rm inj} \times (0,T)$, and $\varphi_g(x,T) = 0$.

Taking into account (76), one easily gets:

where $\Theta \overset{def}{\mathop =} \varrho_g(p_{g})(1-S)$ (see (9) above).

Considering the definition of the functions $\varrho_g$ and $\varrho^\delta_g$ (see (4) and (39), respectively) we have

In order to pass to the limit in the third term of (50) we recall that (see relations (27))

Then

Now taking into account the convergence results (72), (74), and (75) we obtain that

Returning now to the gas pressure function $p_{g}$, we finally get:

Finally, in view if (77),

Thus the gas pressure equation (35) is obtained. Theorem 3.1 is proved.

4.   Homogenization result.

In this Section we consider the problem describing a reservoir with a periodic microstructure. Then in the model considered in the previous sections one has rapidly oscillating porosity function and absolute permeability tensor. Our goal is to prove the homogenization result for this model. The convergence of the homogenization process is justified by the technique of two-scale convergence [2].

4.1.   Formulation of the microscopic problem.

In this section, we present the mathematical model describing water-gas flow in a periodically heterogeneous porous medium. As above we suppose that the gas density vanishes as the gas pressure is zero. For simplicity, we assume no source/sink terms.

We consider a bounded Lipschitz domain $\Omega \subset \mathbb{R}^d$ ($d =1, 2, 3$) with a periodic microstructure. The microscopic length scale that represents the ratio between the cell size to the size of the whole region $\Omega$, is denoted by ${\varepsilon }$. We assume that $0<{\varepsilon } \ll 1$ is a small parameter tending to zero. We denote by $Y \overset{def}{\mathop =} (0, 1)^d$ the periodic cell. Let $\Phi^{\varepsilon }(x) = \Phi(x/{\varepsilon })$ be the porosity of $\Omega$; $K^{\varepsilon }(x) = K(x/{\varepsilon })$ be the absolute permeability tensor of $\Omega$; $S^{\varepsilon } \overset{def}{\mathop =} S^{\varepsilon }_{w}(x, t)$, is the saturations of water in $\Omega\times (0,T)$; $p^{\varepsilon }_{w} = p^{\varepsilon }_{w}(x,t)$, $p^{\varepsilon }_{g} = p^{\varepsilon }_{g}(x,t)$ are the pressures of water and gas in $\Omega\times (0,T)$, respectively;

System (8), in the case of a periodic porous medium, takes the form

where

System (92) has to be equipped with appropriate boundary and initial conditions.

Boundary conditions.: We suppose that the boundary $\partial \Omega$ consists of two parts $\Gamma_{\rm inj}$ and $\Gamma_{\rm imp}$ such that $\Gamma_{\rm inj} \cap \Gamma_{\rm imp} = \emptyset$, $\partial \Omega = \overline\Gamma_{\rm inj} \cup \overline\Gamma_{\rm imp}$. The boundary conditions are given by:

where the velocities $\vec q^{\,{\varepsilon }}_{w}, \vec q^{\,{\varepsilon }}_{g}$ are defined as follows:

Initial conditions.: The initial conditions read:

Let us formulate the main assumptions on the data. First, we replace conditions (A.1), (A.2) from Section 2.2 with the following assumptions:

(A.1$_{\boldsymbol{\varepsilon }}$) The function $\Phi = \Phi(y)$ is $Y$-periodic, $\Phi \in L^{\infty}(Y)$, and there are positive constants $\phi_1, \phi_2$ such that $0 < \phi_1 \leqslant \Phi(y) \leqslant \phi_2 < 1$ a.e. in $Y$.

(A.2$_{\boldsymbol{\varepsilon }}$) The tensor $K = K(y)$ is $Y$-periodic, it belongs to $(L^{\infty}(Y))^{d\times d}$. Moreover, there exist positive constants $K_-, K_+$ such that

We also suppose that conditions (A.3)-(A.9) from subsection 2.2 hold true.

We now provide a weak formulation of problem (92)-(96).

Definition 4.1. For each ${\varepsilon } > 0$ we say that $\langle p^{\varepsilon }_{w}, p^{\varepsilon }_{g}, S^{\varepsilon }\rangle$ is a weak solution of problem (92), (94), (96) if (30)-(38) are fulfilled for functions $ p^{\varepsilon }_{w}, p^{\varepsilon }_{g}, S^{\varepsilon }$ instead of $p_{w}, p_{g}, S$.

Let us recall that for any $\varphi_w, \varphi_g \in C^1([0, T]; H^1(\Omega))$ satisfying $\varphi_w = \varphi_g = 0$ on $\Gamma_{\rm inj} \times (0,T)$ and $\varphi_w(x,T) = \varphi_g(x,T) = 0$, we have:

and

where the function $\Theta^{\varepsilon }$ is defined in (93).

Notational convention. In what follows $C, C_1,..$ denote generic constants that do not depend on ${\varepsilon }$.

4.2.   Statement of the homogenization result.

We study the asymptotic behavior of the solution to problem (92), (94), (96) as ${\varepsilon } \to 0$. In particular, we are going to show that the effective model reads:

where $S$, $P_{w}$, $P_g$ denote the homogenized water saturation, water pressure, and gas pressure, respectively. $\langle\Phi\rangle$ denotes the mean value of the function $\Phi$ over the cell $Y$. $\mathbb{K}^\star$ is the homogenized tensor with the entries $\mathbb{K}^\star_{ij}$ defined by:

where the function $\xi_j$ is a $Y$-periodic solution of the following local problem:

with $\vec e_j$ being the $j$-th coordinate vector.

The function $\Theta^{\star}= \Theta^{\star}(x, t)$ is given by: $ \Theta^{\star} \overset{def}{\mathop =} (1 - S)\, \varrho_{g}\,(P_{g}). $

Remark 2. The homogenized system (100) generalizes the result obtained earlier in [4] in two ways. First, this system allows the gas density to degenerate. In addition, this system is written in terms of the homogenized phase pressures $P_{w}, P_{g}$ and not in terms of the homogenized global pressure and water saturation as was done in [4] (see (3.1)).

System (100) has to be completed with the following boundary and initial conditions.

Boundary conditions.: The boundary conditions are given by:

where the velocities $\vec q^{\,{\varepsilon }}_{w}, \vec q^{\,{\varepsilon }}_{g}$ are defined as follows:

Initial conditions.: The initial conditions read:

4.3.   Proof of Theorem 4.3.

The rigorous justification of the homogenization process relies on the two-scale convergence approach, see, e.g., [2]. For the reader's convenience, we recall the definition of the two-scale convergence.

Definition 4.2. 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 $\Vert v^{\varepsilon }\Vert_{L^2(\Omega_T)} \leqslant C$, and for any test function $\varphi \in C^\infty(\Omega_T; C_\#(Y))$ the following relation holds:

This convergence is denoted by $v^{\varepsilon }(x, t) \stackrel {2s} \rightharpoonup v(x, y, t)$.

The homogenization result reads.

Theorem 4.3. Let assumptions (A.1$_{\boldsymbol{\varepsilon }}$), (A.2$_{\boldsymbol{\varepsilon }}$), (A.3)-(A.9) be fulfilled. Then a solution of problem (92) converges (up to a subsequence) to a weak solution of the homogenized problem (100).

The proof is divided into a number of steps.

4.3.1. A priori estimates for solutions to problem (92)

In this section we derive the a priori estimates for problem (92). For any ${\varepsilon } >0$, we consider the following ${\varepsilon ,\delta }$-problem:

where the family of functions $\{\varrho^\delta_g\}_{\delta>0}$ is defined in (39) and

The $({\varepsilon ,\delta })$-problem is completed by the boundary and initial conditions (94) and (96), respectively.

The energy equality for problem (106) can be obtained as in Section 3.3.

Lemma 4.4 (Energy equality for $({\varepsilon ,\delta })$-problem. Let $\langle p^{\varepsilon ,\delta }_{w}, p^{\varepsilon ,\delta }_{g}, S^{\varepsilon ,\delta } \rangle$ be a solution to (106). Then

in the sense of distributions. Here

where the functions $\digamma(s)$ and ${\mathcal{R}}^\delta(p)$ are defined by (56).

Then following the lines of Section 3.3 one can prove the following statement which is similar to that of Lemma 3.5.

Lemma 4.5. Let $\langle p^{\varepsilon ,\delta }_{w}, p^{\varepsilon ,\delta }_{g} \rangle$ be a solution to (106), the global pressure ${\mathsf P}^{\varepsilon ,\delta }$ be defined in (16), and the function $\beta(s)$ be defined in (25). Then

Here $C$ does not depend on ${\varepsilon ,\delta }$.

Now, as in Section 3.5.1, we conclude that

where $C$ is a constant that does not depend on ${\varepsilon ,\delta }$.

The uniform estimates for the time derivatives of the functions $\Phi^{\varepsilon } \Theta^{\varepsilon }$ and $\Phi^{\varepsilon } S^{\varepsilon }$ can be derived from (92) using (115). These estimates read:

where $C$ is a constant that does not depend on ${\varepsilon }$.

4.3.2. Compactness and convergence results

First, we recall the following compactness result established in~[4].

Lemma 4.6 (Compactness lemma). Let $\Phi \in L^\infty(Y)$, and assume that there are positive constants $\phi_1, \phi_2$ such that $0 < \phi_1 \leqslant \Phi(y) \leqslant \phi_2 < 1$ a.e. in $Y$. Assume moreover that a family $\{v^{\varepsilon } \}_{{\varepsilon }>0} \subset L^2(\Omega_T)$ satisfies the following properties:

1. $v^{\varepsilon }\in L^\infty(\Omega_T)$, and $0 \leqslant v^{\varepsilon } \leqslant C$;

2. there exists a function $\varpi$ such that $\varpi(\xi) \rightarrow 0$ as $\xi\to 0$, and

3. $\left\Vert \partial_t(\Phi^{\varepsilon } v^{\varepsilon }) \right\Vert_{L^2(0,T;H^{-1}(\Omega))} \leqslant C$.

Then the family $\{v^{\varepsilon } \}_{{\varepsilon }>0}$ is a compact set in $L^2(\Omega_T)$.

Remark 3. In the formulation of the above compactness lemma the periodicity of $\Phi$ can be replaced with the assumption that $\Phi^{\varepsilon } \rightharpoonup 1$ weakly in $L^2(\Omega)$, as ${\varepsilon }\to0$.

Now we turn to the compactness result for the family $\{\Theta^{\varepsilon }\}_{{\varepsilon }>0}$.

Proposition 3. Under our standing assumptions, the set $\{\Theta^{\varepsilon }\}_{{\varepsilon }>0}$ is compact in the space $L^q(\Omega_{T})$ for all $q \in [1, +\infty)$.

A similar result holds for the set $\{S^{\varepsilon }\}_{{\varepsilon }>0}$.

Proposition 4. Under our standing assumptions, the set $\{S^{\varepsilon }\}_{{\varepsilon }>0}$ is compact in the space $L^q(\Omega_T)$ for all $q \in [1, +\infty)$.

Summarizing the above statements yields.

Lemma 4.7 There exist a function $S$ with $0 \leqslant S \leqslant 1$ a.e. in $\Omega_T$ and a function ${\mathsf P} \in L^2(0, T; H^1(\Omega))$ such that up to a subsequence:

The Proof of Lemma 4.7 relies on the arguments similar to those used in the proof of Lemma 4.8 in [4].

4.3.3. Passage to the limit in equations (98), (99)

In this subsection we apply the method of a cut-off function introduced in [10].

It is easy to justify the passage to the two-scale limit in the temporal terms using the convergence results (117) and (120) from Lemma 4.7 as it was done, for example, in [4]. Namely, let $\varphi_0 \in {\mathcal{D}}(\Omega_T)$. The first two terms in (98) become:

Now we pass to the limit on the right-hand side of (121). Taking into account (117), we have that

For any $\eta > 0$, we introduce the family of functions $\{S^{{\varepsilon },\eta}\}$ defined by:

These functions satisfy the estimate:

where $C(\eta) \to + \infty$ as $\eta \to 0$. Therefore,

Now, taking into account (117), (118), for a subsequence,

with ${\mathsf V}^\eta_w \in L^2(\Omega_T; H^1_\#(Y))$. We set:

with ${\mathsf Z}(s)$ being a smooth function equal to zero for $s\not\in(\eta, 1 - \eta)$; $\zeta(y)$ is smooth periodic, and $\varphi$ is a smooth function with a compact support in $\Omega_T$. Using $\varphi^{\varepsilon }_w$ as a test function in (98) and considering the global pressure definition, we get:

We pass to the two-scale limit in (125). Taking into account (117), (118), and (123), we obtain:

Therefore,

for all $(x,t)\in\Omega_T$ such that $S \in (\eta, 1 - \eta)$. Here $\xi \in \mathbb{R}^d$ is a vector with the components $\xi_j$ that are the solutions of the auxiliary problem (102).

Since $\eta$ is an arbitrary positive number, representation (127) is valid for all $(x,t)$ such that $S \in (0,1)$. In particular, ${\mathsf V}^\eta_\ell$ does not depend on $\eta$: ${\mathsf V}^\eta_w = {\mathsf V}_w$. This leads to the following equation:

Finally, with the help of our a priori estimates we deduce in a standard way that

where $\mathbb{I}$ is the unit matrix. This allows us, with the help of (122), to obtain the weak formulation of the homogenized saturation equation (100)$_2$.

The derivation of the weak formulation for the homogenized gas pressure equation can be done in a similar way. This completes the proof of Theorem 4.3.

5.   Concluding remarks.

We have presented new results for immiscible compressible two-phase flow in porous media. More precisely, we give a week formulation and an existence result for a degenerate system modeling water-gas flow through a porous medium. The water is assumed to be incompressible and the gas phase is supposed compressible and obeying the ideal gas law leading to a new degeneracy in the evolution term of the pressure equation. Furthermore, a homogenization result for the corresponding system is established in the case of a single rock-type model. The extension to a porous medium made of several types of rocks, i.e. the porosity, the absolute permeability, the capillary and relative permeabilities curves are different in each type of porous media, is straightforward. Let us also mention that this homogenization result has been used successfully in [1] to simulate numerically a benchmark test proposed in the framework of the European Project FORGE: Fate Of Repository Gases [28]. The study still needs to be improved in several areas such as the cases of unbounded capillary pressure and double porosity media. These more complicated cases appear in the applications. Further work on these important issues is in progress.

Acknowledgments

Most of the work on this paper was done when L. Pankratov and A. Piatnitski 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 has been partially supported by the RScF, project No. 15-11-00015. This work was partially supported by the Carnot Institute, ISIFoR project (Institute for the sustainable engineering of fossil resources). The supports are gratefully acknowledged.