Homogenized Model of Two-Phase Flow with Local Nonequilibrium in Double Porosity Media

We consider two-phase flow in a heterogeneous porous medium with highly permeable fractures and low permeable periodic blocks. The flow in the blocks is assumed to be in local capillary disequilibrium and described by Barenblatt’s relaxation relationships for the relative permeability and capillary pressure. It is shown that the homogenization of such equations leads to a new macroscopic model that includes two kinds of long-memory effects: the mass transfer between the blocks and fractures and the memory caused by the microscopic Barenblatt disequilibrium. We have obtained a general relationship for the double nonequilibrium capillary pressure which represents great interest for applications. Due to the nonlinear coupling and the nonlocality in time, the macroscopic model remains incompletely homogenized in general case. The completely homogenized model was obtained for two different regimes. The first case corresponds to a linearized flow in the blocks. In the second case, we assume a low contrast in the block-fracture permeability. Numerical results for the two-dimensional problem are presented for two test cases to demonstrate the effectiveness of the methodology.


Introduction
The classical theory of multiphase flow of liquids and gases through porous media is based on the hypothesis of the local thermodynamic equilibrium. This means that, in each representative elementary volume (REV) of a porous medium, the thermobaric, the capillary, and the chemical equilibrium between the fluid phases is established instantly. Until the phases form a small-scale mixture in porous space with the characteristic size of droplets (bubbles, ganglions, menisci, films, etc.) about the radius of pores, the thermobaric and chemical equilibrium between two phases is totally acceptable, which is not however the case of the capillary equilibrium.
The capillary equilibrium determines a special distribution of two phases in the pore network such that the resulting capillary force is nil. According to the classical theory, such an equilibrium phase distribution is reached when the wetting fluid prefers to occupy narrow pores, while the nonwetting fluid is confined in large pores. However if we displace oil by water, the displaced oil (nonwetting) ahead of the front occupies all the pores including narrow, while the injected water (wetting) behind the front moves in all the pores including the large ones. Such a fluid distribution does not correspond to the equilibrium state; then a relaxation process spontaneously starts in the medium and tends to redistribute the two fluids over the pore network. This redistribution consists of the rearrangement of phase clusters, implies the fluid replacement to long distances, which concerns the overall fluid, and, consequently, is nonlocal. The redistribution of phases represents the intrinsic flow of both phases caused by the capillary forces (applied to the curved interfaces between two phases) and slowed by the friction against the pore walls. Such a flow occurs within a pore network 2 Advances in Mathematical Physics having usually a low connectivity, which reduces the rate of phase redistribution. The lower the medium permeability, the higher the friction, and the lower the connectivity of the pore network. Consequently, the redistribution time in low permeable media is expected to be large, such that the hypothesis of the local capillary equilibrium can no longer be accepted.
The modeling of such redistribution/relaxation was introduced first by Barenblatt [1] with modifications in [2,3] and developed further by a number of authors. In [1], it is assumed that the relative permeability ( , ) and the capillary pressure ( ) depend not on the true instantaneous saturation , but on an effective saturation that corresponds to the equilibrium value that might be reached by the system after a phase redistribution in space: where is the redistribution time. Note that / is the rate of the phase redistribution; therefore ( / ) means the saturation variation from a nonequilibrium value up to the equilibrium one. Theoretical justifications of relationship (1) have been obtained by several authors using the approach of the nonequilibrium thermodynamics: first by Nikolaevskii et al. [4], published only in Russian, by Marle in [5], and later by Hassanizadeh and Gray [6]. The method consists in formulating the relationship for the entropy production as a bilinear form made up of the products between thermodynamic forces and fluxes and assuming Onsager's linear relationships between the forces and fluxes of the same tensor dimension.
Let us now consider a double porosity medium, which represents a connected network of highly permeable porous channels ("the fractures") and very low permeable inclusions ("the blocks"). As mentioned above, the capillary nonequilibrium should be retained in the blocks, which is the objective of the present paper.
The two-phase flow in double porosity media at local equilibrium is well studied; see, for example, [7][8][9]. It was shown that, in such systems, despite the local equilibrium, another nonequilibrium phenomenon appears on the macroscale, which is caused by the high delay in capillary redistribution of the fluids between blocks and fractures. This disequilibrium is described in terms of long-memory operators, which appear in the macroscopic model of flow. In [8,9], the obtained general macroscopic model was not completely homogenized. A completely homogenized model was obtained in the same papers [8,9] by the linearization of the local problem in the block, which allowed obtaining an explicit integrodifferential expression for the memory. Consequently the differential microscopic flow equations transform into integrodifferential macroscopic equations. Another completely homogenized model (always without Barenblatt's nonequilibrium) was obtained in [10,11] for the case of a lower contrast between the block and fracture permeability. This case leads to a lower degree of the disequilibrium, described by a model with short memory. The macroscopic model was completely homogenized without any linearization.
If now we introduce the microscopic Barenblatt capillary nonequilibrium in a double porosity medium, then the microscopic and the macroscopic nonequilibrium effects will interact with each other. The physical mechanism of this interaction is still unclear, and its mathematical model is unknown. The analysis of this interaction is the main objective of the present paper. More precisely, our aim is to develop a macroscopic model for two-phase incompressible flow in a double porosity medium including the local nonequilibrium.
The two-phase flow in double porosity medium with local nonequilibrium has been studied in [12] where upscaling was performed numerically.
Usually the local capillary nonequilibrium is introduced through the dynamic capillary pressure relationship. Another approach is based on Kondaurov's phenomenology developed for viscoelasticity. This approach was applied to twophase flow in double porosity media in [13,14].
Mathematical study of the immiscible two-phase flow with local nonequilibrium in porous media concerns, for the moment, only homogeneous media. The local nonequilibrium is introduced through the dynamic capillary pressuresaturation relationship. Existence and uniqueness results can be found in [15][16][17][18], traveling wave solutions have been obtained in [19][20][21], and the numerical analysis of the nonequilibrium models is presented in [22,23].
In this paper, the homogenization method is used to derive the upscaled model for two-phase flow with local Barenblatt's nonequilibrium in double porosity media with periodic microstructure. More precisely, we have obtained the following new results: (i) a completely homogenized model for the case of a linearized flow in blocks; (ii) a completely homogenized model for the case of low contrast in the block-fracture properties; (iii) an explicit relationship for calculating the memory functions and the relaxation times; (iv) the qualitative effects of the nonequilibrium effects described by double memory operators (double convolution integrals).
The rest of the paper is organized as follows. In the next section, we formulate the two-phase flow problem with local nonequilibrium effects on the microscale and formulate the main assumptions on the data. Section 3 is devoted to the presentation of the homogenization method. First, in Section 3.1 we present the two-scale asymptotic expansions. Secondly, the macroscopic variables are defined in Section 3.2. Thirdly, in Section 3.3 we obtain a homogenized model. Finally, in Section 3.4 we derive, in addition, a relationship between the average capillary pressures in blocks and fractures. In Section 4 we are dealing with the first completely homogenized model obtained for the case of a linearized cell problem in the blocks. The second completely homogenized model is Advances in Mathematical Physics 3 formulated in Section 5 for the case of low contrast in blockfracture properties. In Section 6, we exhibit the numerical results for two test cases. We present numerical simulations to analyze the impact of the nonequilibrium effects for twophase flow in a fractured reservoir. Additional conclusions are drawn in Section 7.

Formulation of the Problem at the Microscale and Main Assumptions
2.1. Medium Geometry and Properties. We consider twophase immiscible flow of oil and water at the Darcy scale in an inhomogeneous porous medium consisting of two types of rocks having different porosities and permeability. The heterogeneities are organized in the form of a periodic system of low permeable blocks immersed into a connected highly permeable subdomain called "fractures." The heterogeneity size is small with respect to the entire domain size. Such a system has an intrinsic macroscopic behavior, independent of the boundary conditions. Consequently, we will distinguish two different scales of consideration: (i) the microscale: the mathematical point corresponds to an elementary representative volume of porous medium; the microscale flow is described by Darcy's law; (ii) the macroscale: the mathematical point corresponds to a period of medium heterogeneity. The macroscopic flow equations may be obtained by the homogenization of the microscale equations.
To describe the microscopic geometry, we will introduce the following notations: the overall domain Ω ⊂ R ( = 2, 3) is bounded and connected and has periodic structure. The heterogeneity period is a small parameter. The low permeable blocks are Ω m , and the connected highly permeable medium (the fracture) is Ω f . The interface between the blocks and the fracture is Γ fm . Let fl (0, 1) be the period of the medium extended by −1 times. Then m , f , and Γ fm are the images of Ω m , Ω f , and Γ fm respectively; see Figure 1. We also introduce the notations: Ω fl Ω × (0, ) and Ω ℓ, fl Ω ℓ, × (0, ), for ℓ = f, m. Note that the fractures are not an empty space but represent a porous medium with permeability much higher than that of the blocks.
The fractures are not thin; that is, their thickness is not an additional small parameter, as in [24].
The double porosity medium is defined through the following properties of the porosity ( ) and the absolute permeability tensor ( ): where the values 0 < f , m < 1, 0 < f , m < max < +∞ are positive constants that do not depend on ; and I is the identity matrix.

Microscale Nonequilibrium Effects.
We assume that the flow in the fractures is in equilibrium, while it is in local disequilibrium in the blocks. This means that the relative permeability and the capillary pressures may be considered as functions only of the local saturation in fractures, while in blocks they depend on a nonlocal functional of saturation which keeps the memory of the process. We adopt the model of Barenblatt with short memory, by introducing the functional m ( m , m ), called the efficient saturation and defined similar to (1): where m is the relaxation time in blocks Ω m , and m is the water saturation in blocks. We assume that the relative permeability and the capillary pressure in blocks depend only on this functional. The degree 4 Advances in Mathematical Physics of the disequilibrium is determined by the relaxation time m : the larger m , the higher disequilibrium in the system.
The relaxation time in fractures is much lower than in blocks, and, consequently, we can assume the local equilibrium in fractures.

Microscale Flow Equations.
Let us introduce the following notations, for ℓ = f, m: (i) ℓ = ℓ, is the water saturation in the medium Ω ℓ (the oil saturation is then 1 − ℓ ); (ii) ℓ, ( , ) and ℓ, ( , ) are the pressures of water and oil in medium Ω ℓ ; (iii) f, ( f ) and f, ( f ) are the relative permeability of water and of oil in fractures Ω f , and they are given functions and depend only on the local saturation; (iv) m, ( m ) and m, ( m ) are the relative permeability of water and of oil in blocks Ω m , and they are given functions and depend only on the nonlocal functional m ; ( (vi) and are the dynamic viscosities which are assumed to be positive constants; (vii) f, ( f ) and m, ( m ) are the capillary pressures in fractures and in blocks.
In what follows, we assume that the capillary pressure functions f, ( f ), m, ( m ) are positive decreasing Lipschitz functions of their arguments. We also assume that there exist constants , ∈ (0, ∞) such that the relative permeability satisfies the bounds:

Microscale Flow Equations.
For the sake of simplicity, we neglect the gravitational effects and the source/sink terms.

Mass Balance in Fractures.
The mass conservation of each phase can be written in the framework of the standard Darcy-Muskat equations [25][26][27][28]. In Ω f, :

Conditions at the Block-Fracture Interface.
We assume the continuity of the phase fluxes and pressures at the interface Γ fm : where ⃗ ] is the unit outer normal on Γ fm . Note that the continuity of phase pressures leads to the continuity of the capillary pressure. In several situations the capillary pressure may be discontinuous, which is the case of the so-called capillary trapping of the nonwetting phase in low permeable blocks [32,33]. These effects are not analyzed in the present paper.

Initial Conditions.
The initial condition in the fracture system read: In the present paper we consider only the situation of the initial equilibrium state, when init m ( ) = init m ( ). In more general case the dependence of init m on init m is determined from the last equation in (7), [31]. Then the initial conditions in the matrix part read: These equations will be called "the microscopic flow model." The conditions at the reservoir boundary are assumed to not influence the homogenized model.

Upscaled Model
In this section we present the upscaled model derived by the two-scale asymptotic expansions. The technique of derivation is presented in Appendix.

Two-Scale Asymptotic Expansions.
We assume that the solutions to problem (6)-(8) depend on the space variable Advances in Mathematical Physics 5 in the following sense: they are functions of the macroscopic (or "slow") variable ∈ Ω and the microscopic (or "fast") variable ∈ . The macroscopic and microscopic scales are related by the small parameter up to a translation as = / . This implies that where ∇ , ∇ are the gradients with respect to the variables , .
We are interested in finding the set ⟨ f, , f, , m, , m, , f , m , m ⟩, which is the solution to (6)-(8) having the following asymptotic expansion: where all the coefficients of these expansions are -periodic in .
For the phase mobility, we obtain the asymptotic expansion (for = , ): Note that the first terms f, , f, , f do not depend on the fast variable , which is explained by the fact that the saturation field in fractures is instantaneously established and, consequently, is uniform within a cell; see, for instance, [7,34].

Macroscopic Variables.
The physical meaning of the first terms of the asymptotic expansions is , and ( , ) are the macroscopic water saturation, water pressure, and oil pressure in fractures; (ii) m ( , , ), ( , , ), and ( , , ) are the nonhomogenized first terms of the asymptotic expansions for the local water saturation, water pressure, and oil pressure in blocks; (iii) m ( , , ) is the nonhomogenized first term for the nonlocal Barenblatt saturation in blocks; (iv) | ℓ | is the volume of the domain ℓ (ℓ = f, m); is the macroscopic tensor of the absolute permeability, homogenized over the fractures, where , 1 ≤ ≤ , is a solution of the following auxiliary cell problem: where ⃗ 1 , . . . , ⃗ is the canonical basis of R .

Macroscopic Model.
The technique of homogenization is presented in Appendix.

The homogenized model reads in
with the initial condition f ( , 0) = init f ( ). This model is determined through the function m ( , , ) that is the local saturation in blocks, which is a nonlocal (in time) functional of the effective saturation m : where ( ⋆ )( ) = ∫ 0 ( − ) ( ) . This is the explicit solution to (4) reformulated for the zero terms of the asymptotic expansion (12).
The function m , in turn, is the solution to the following problem within a block ∈ m : The properties of the phase mobility should ensure the positive definition of the capillary diffusion parameter m .
Remark 1. The model (15)-(17) is not completely homogenized. Indeed, the macroscopic equations (15) contain the microscopic function m ( , , ), which can be determined only after solving the cell problem (17). On the other hand, the cell problem (17) contains the macroscopic function f ( , ) through the boundary condition, which cannot be eliminated. This means that this problem must be solved for all the cells.

Remark 2.
In the case of the local equilibrium in the blocks we have m = S m = m , according to Section 2.2. Then the macroscopic model (15)-(17) takes the form of the double porosity model for two-phase flow obtained by many authors; see, for example, [11,24,35,36].  (17); thus, m , in turn, represents a nonlocal functional. Consequently, the source term Q in (15) has a double nonlocality.

Doubly Nonequilibrium
Relationship for the Macroscopic Capillary Pressure. Although the macroscopic model is incompletely homogenized, it enables obtaining a general nonequilibrium relationship for the macroscopic capillary pressure, which is sufficient in practice to analyze qualitatively the overall process. It has the following form: where the function = ( , , ) is the solution of the following cell problem: where To prove (19), let us multiply (17) by where ( , ) is a nonzero infinitely differentiable function with a compact support in the cylinder Ω . Integrating over m × Ω we obtain the following relationship: Integrating by parts, using Gauss's theorem, taking into account (20) and the fact that m, which is equivalent to (19). Note that (19) generalizes the nonequilibrium capillary relationship obtained earlier in [9] for the case of local equilibrium.

Completely Homogenized Model I: Linearized Flow in Blocks
A linearization of the cell problem (17) is a way to obtain an explicit analytical relationship for the function m , then for m , and then for the source term Q in (15) and finally to obtain a completely homogenized model.

Linearized Nonequilibrium Cell Problem in Blocks.
The idea of the linearization and, thus, simplification of the local problem comes from [8]. This approach already enabled obtaining a completely homogenized two-phase flow model in double porosity media in [9,24]. In this paper we make use of the constant linearization; that is, we set where the function m is defined in (18). The linearized imbibition equation (17) reads where we keep the same notation m as in (17) for the function satisfying the linearized problem (26). can be presented in the following form with separated fast variable and slow variable :

Separation of the Fast and Slow Variables in the Linearized
where the function w = w( , ) is independent of the slow variable and is the solution of the following auxiliary cell problem: In order to prove (27), we proceed as follows. We will plug the nonequilibrium parameter m given by (27) in (26) and show that (26) is satisfied. First, we note that it follows from (16) that and from (26), we get Now, taking into account the definition of the function w given in (28) we observe that the nonequilibrium parameter m given by (27) is the solution of the boundary value problem (26). This completes the proof of the representation (27).

Completely Homogenized Model I.
Using the obtained solution for the cell problems (27) and (28), as well as the relationship (30), we obtain immediately for the exchange term Q in the macroscopic model (15) (33) where the auxiliary function w( , ) is the solution of the local problem (28). It is sufficient to substitute (28) into (15) to obtain (33). Then the macroscopic model (15) takes the form with the initial condition f ( , 0) = init f ( ). The local saturation in blocks m ( , , ) is calculated as The averaged saturation in blocks ⟨ m ⟩ m and the averaged Barenblatt pseudo-saturation ⟨ m ⟩ m are calculated as The averaged capillary pressure in the blocks can be calculated directly as ⟨ m, ( m )⟩ m since the function m, ( ) is given.
The general relationship between the capillary pressures (19) becomes then useless. In addition, it remains noncompletely homogenized, since the function ( , , ) defined by problem (20) remains dependent on m and so forth.

Remark 4.
(1) Although the problem is completely homogenized, the calculation of the averaged functions and coefficients requires solving the cell problem (28) for the integrodifferential equation.
(2) The double convolution integral in the macroscopic model (34) clearly shows the appearance of the double memory phenomenon, one of which is caused by the blockfracture interaction while the second one is produced by the microscale nonequilibrium in blocks.

Completely Homogenized Model II: Low Contrast in Block-Fracture Properties
The second way to obtain a completely homogenized model is to assume that the block permeability 2 m is still very low 8 Advances in Mathematical Physics with respect to fractures, but the contrast in the permeability is not so high as previously. Then the parameter m in (3) should be sufficiently high. Its order, however, should not disturb the asymptotic expansions (12), which is ensured if The fact that the contrast in block-fracture properties is not very high means that the delay in the block behavior is expected to be not very high. Thus we deal with the case of a low macroscopic nonequilibrium.
The macroscopic model in this case may be obtained from the general model (15), (16), (17) by asymptotic expansion with respect to parameter .
The justification to the fact that the expansion over two small parameters , is commutative may be done in the following way. According to [7], the block's contribution to the macroscopic model is double: (i) they form a source flow from the block center towards fracture, which determines the exchange terms Q in (15); and (ii) they can contribute to the effective permeability if the block permeability is not too low. For the case of the block permeability equal to 2 m with m ∼ 1, the macroscopic model (15) contains only the source term but does not contain the block contribution to the effective permeability (ctep). If we progressively increase the block permeability, then the source flow will decrease, but ctep will grow. The source flow and ctep become of the same order when m ∼ −1 . If 1 ≪ m ≪ −1 (or ≪ ≪ 1), the source flow is more important than ctep. In this case, the ctep may be still neglected; then the macroscopic model is analogous to (15), but the source terms will have smaller order and may be expand asymptotically over . This corresponds exactly to the case (37).

Asymptotic Expansions
wherêf,̂m, and so forth are the average values over the fractures and the block, respectively. Consequently to such a structure of the expansion, the functions of the first approximation satisfy the following property: According to [11], in this case the first term̂m( , ) does not depend on the fast variable , differently from the case which is justified by the low contrast in the medium heterogeneity and, consequently, by the sufficiently fast stabilization of the saturation field in blocks.

Remark 5.
The expression (40) is a generalization of the nonequilibrium relationship (38c) of Bourgeat and Panfilov obtained in [11], which takes into account the double nonequilibrium: the microscopic phase redistribution in the blocks characterized by the microscopic relaxation time m and the macroscopic relaxation time macr characterized by the delay in phase redistribution between the blocks and the fractures. Thus, the low contrast in medium properties, which corresponds to sufficiently weak nonstationary effects within the blocks, leads to the appearance of a short macroscopic memory in the upscaled system, with small relaxation time macr . (15) and (40) we obtain the following completely upscaled model: with the initial condition̂f( , 0) = init f ( ). This system of five equations represents the closed completely homogenized model with respect to five functionŝf, m ,̂m,̂, and̂.

Completely Homogenized Model II. From
Note that the latter relationship for the capillary pressure in (47) represents the equation that defines the saturation in blockŝm (as the solution of the nonlinear ordinary differential equation of the first order).

Double Nonequilibrium Capillary Relationship for Small
where macr is defined through (41). Then, in the case of low contrast in medium properties and low microscopic disequilibrium, we obtain the kinetic relationship for the capillary pressure in blocks, in which the general relaxation time represents a weighted sum of the microscopic and the macroscopic relaxation times.

Analysis of the Impact of the Double Nonequilibrium
To better understand the impact of the local nonequilibrium on the macroscopic behavior, we consider two test cases based on the completely homogenized model (47). Let us assume that the saturation in fractures is constant in space and is a given function of time, whose behavior reflects the physical process we are analyzing. The capillary pressures are described by the functions We analyzed two tests, in which the initial state was in local disequilibrium within blocks and in macroscale disequilibrium between blocks and fractures. We maintained a given behavior of the saturation in fractures and analyzed the answer of the blocks on such a perturbation. The saturation in blocks,̂m, was calculated from the differential equation (48): Test Case 1. In the first test, the saturation in fractures, f ( ), was constant in time and equal to 0.155. The corresponding equilibrium saturation in blocks was 0.455; however the initial saturation in blocks, init m , was equal to 0.2, which was in high disequilibrium.
Such a disequilibrium causes the spontaneous evolution of the saturation in blocks, which tends to the equilibrium value, as shown in Figure 4.
As expected, the time of reaching the total equilibrium is higher in the case of double disequilibrium.
Test Case 2. In the second test, the saturation in fractures, f ( ), periodically oscillated in time, which corresponds in practice to the functioning of an underground storage of injected gas in an aquifer. In this case, the saturation is that of the injected gas. The increase of the injected saturation means the period of gas injection, while the saturation minimums correspond to gas extraction. If the system was in total equilibrium, then its behavior would be such that is shown in Figure 5.
The true behavior of the saturation in blocks is shown in Figure 6, where the initial saturation in blocks, init m , was in the nonequilibrium with fractures.  As seen, the saturation in blocks is in permanent delay with respect to the equilibrium curve. The delay is expressed in two main forms: (i) the minimum and the maximum values of saturation are very different from the equilibrium case; (ii) the phase-shift is observed for the nonequilibrium case.
Both effects are higher for the case of the double nonequilibrium.
Therefore, the local Barenblatt disequilibrium certainly increases the degree of the nonequilibrium in the system but also can lead to very different qualitative behavior of the saturation field, such that the equilibrium behavior can be never reached, as in the example with underground gas storage.

Conclusion
We have developed, using homogenization, a new effective model of two-phase flow in double porosity media which takes into account the local disequilibrium of phase redistribution, introduced first by Barenblatt. This process can be very special in low permeable media. The obtained model is more general than those developed previously in [8,9] for the case of the local equilibrium.
In the case of low contrast in block-fracture properties, the model may be completely homogenized and reduced to system (47), in which the last relationship relates the average capillary pressures in blocks and fractures. The total relaxation time is explicitly calculated as a function of the average saturation in blocks. It includes two terms: the microscale time of Barenblatt and the macroscale time of delay between block and fractures. These two components of the relaxation times underline the double nonequilibrium in the system: one of them appears on the microscale in the form of the noninstantaneous phase redistribution, while the second one is the permanent delay between the saturations within the blocks and the fractures.
The simulations we have performed show that the nonequilibrium can have very important impact on the behavior of the fluids, especially in the case of oscillating regimes typical, for instance, for underground gas storages. In this case, the average saturation follows a strict nonequilibrium behavior and can never reach the equilibrium state. where the functions are not defined but will disappear from the average model, while the functions ( = 1, . . . , ) are the solutions to the cell problems (14) resulting from (A.1).
(4) Substituting (A.9) into (A.7) we obtain the final form of the averaged equations: which is the macroscopic model (15). Herein K ⋆ is defined in Section 3.2. which is (17).

Competing Interests
The authors declare that they have no competing interests.