Discretization schemes for the two simplified global double porosity models of immiscible incompressible two-phase flow

We present the discretization schemes for the two simplified homogenized models of immiscible incompressible two-phase flow in double porosity media with thin fractures. The two models were derived previously by the authors by different linearizations of the nonlinear local problem called the imbibition equation which appears in the homogenized model after passage to the limit as ε → 0. The models are fully homogenized with the matrix-fracture source terms expressed as a convolution.


Introduction
Naturally fractured reservoirs are often met in petroleum engineering and hydrology applications, for instance in secondary oil recovery in hydrocarbon reservoirs, CO 2 sequestration or the nuclear waste management.Naturally fractured reservoirs form a particular type of porous media which are made of two superimposed continua: a network of disconnected periodically staggered middle-sized matrix blocks of ordinary porous media is intertwisted on a fine scale by a interconnected system of thin fractures.In such system the fractures are considerably more permeable compared to the matrix blocks.Therefore the effective permeability of a reservoir is higher than the permeability of the rock matrix solely.The fluid flow through the fractured reservoir occurs foremost inside the fracture system where the fluids are much more mobile than in the porous rock; on the other hand the porous matrix stores most of the fluid.Furthermore, the typical width of the fractures is much smaller than the size of the reservoir.The mentioned discrepancies give rise to large difficulties in numerical modelling of multiphase flow in fractured porous media since at the same time fractures and matrix as well as their interactions must be incorporated in a flow model.
In double porosity or dual-continuum models of multiphase flow in fractured porous media one does not consider individual fractures but instead a network of tiny interconnected fractures.In such systems two large differences in scale are present: fracture width is very small compared to the scale of the domain and fracture permeability is much larger than the permeability of surrounding porous rock.Certain averaging or homogenization processes are applied in order to obtain simplified description of matrix-fracture system and its interactions on a global scale of a reservoir.
The method of two small parameters ε, δ for homogenization of the two-phase flow in double porosity media was proposed in [1,2].Here ε is the periodicity parameter of the matrix blocks system while δ is the width of the fractures relative to the size of a rock matrix block.The homogenization procedure of this method includes two steps: first perform homogenization as ε → 0 and then pass to the limit as δ → 0. As further references on this procedure we refer to [3,4] and the references therein.The first result on the homogenization of the two-phase flow in double porosity media was obtained in [5] by the method of two small parameters.The method of two small parameters was later used also in [6] where a homogenized model of the two-phase non-equilibrium flow in fractured porous media was derived.
In the global double porosity δ-problem obtained in [5] after passing to the limit as ε → 0 there appear additional matrix-fracture source terms that are defined implicitly via solutions of a nonlinear local boundary value problem known as the imbibition equation.The nonlinearity of the imbibition equation presents difficulties in numerical simulations of the model since no analytic expression of the matrix-fracture source term is available.In order to overcome this issue, one can linearize the imbibition equation and then express the matrixfracture source term explicitly from the linearized equation.In [5] the imbibition equation is linearized by using an appropriate constant, as suggested by [7].Then in [8] we present a new, variable and more general linearization of the imbibition equation which gives a new simplified double porosity model.In this paper we provide the dicretization schemes for the two simplified double porosity models obtained by the two linearizations.
The rest of the paper is organized as follows.In Section 2 we present the two fully homogenized simplified double porosity models of the two-phase flow in the case of thin fractures.The main contributions of the paper are contained in Sections 3 and 4 in which we propose the discretization schemes of the effective system including the matrix-fracture source terms in the convolution form obtained by constant linearization and by variable linearization, respectively.

Global simplified double porosity models I and II
In this section we present two global double porosity -type models of incompressible two-phase flow in fractured porous media derived previously by the authors in [5] and [8], respectively.In Sections 3 and 4 we will present the numerical schemes for these models.
Let the reservoir Ω ⊂ R d , d = 1, 2, 3, and we denote Ω T = Ω × (0, T ), with T > 0 fixed.In the derivation of our models it is assumed that the reservoir Ω is composed of the matrix blocks and the highly permeable connected network of fractures.The periodicity of the fractured porous medium is depicted by a small parameter ε ≪ 1 representing the characteristic matrix block size which is small with regards to the size of the flow domain.A second small parameter δ, satisfying 0 < ε ≪ δ < 1, describes the relative thickness or opening of the fractures so that the fractures' thickness is assumed to be of order εδ.The reference cell is Y = (0, 1) d and we take the domain Ω to be paved by cells εY (see Figures 1, 2).The porosities of the blocks and the fractures are supposed to be constant and are denoted by Φ m and Φ f , respectively.The permeabilities of the blocks and the fractures are highly contrasted.In the double porosity model one assumes that the fracture porosity is k f , while the matrix porosity is (εδ) 2 k m , where k f and k m are of the same order.
The starting point in derivation of our global simplified models is the system describing the two-phase incompressible flow in Ω T , which depends on two small parameters ε and δ.In order to obtain the effective model, the first step is to pass to the limit as ε → 0, which has been proved rigorously in [9,10,11].The resulting global δ -double porosity model is not fully homogenized (by the definition of [12]) since the effective permeability and the matrixfracture source term are defined via solutions of certain coupled local problems.However, by letting the small parameter δ to zero one can obtain full decoupling of the global and the local problem.For the case of the two-phase flow this problem was first studied in [5] by performing the linearization of the imbibition equation by a constant, as suggested by [7] (see the details in [5]).It is shown in [5] that in the limit as δ → 0 the effective fracture equations coupled with the linearized imbibition equation obtained by this constant linearization reduce to the system: Here S, P w and P n are the wetting phase saturation, the wetting phase pressure and the non wetting phase pressure in the system of fractures, respectively; P c,f , λ w,f and λ n,f are the capillary pressure function and the phase mobility functions in the fractures; is the reduced fracture permeability.The terms Q w and Q n are the wetting phase and the non wetting phase source terms modeling the phase mass transfer from the matrix to the fracture system governed by the capillary imbibition and are given by Here P c,m , λ w,m and λ n,m are the capillary pressure function and the phase mobility functions in the matrix and the total mobility in the matrix λ m (s) = λ w,m (s) + λ n,m (s).
Let us note that in this final effective model the effective porosity is the fracture porosity, the effective permeability is reduced fracture permeability, and the matrix-fracture source term Q w (t) is expressed as the convolution (2) with the kernel The model ( 1), ( 2) (double porosity model I) is fully homogenized, in other words the local and the global problems are completely decoupled, since the effective permeability and the matrix-fracture source term are given explicitly and there is no need to solve local problems.
The new, more general approach to the linearization of the imbibition equation presented in [8] leads to our second effective model (double porosity model II).In this case a variable linearization of the imbibition equation was proposed.The effective model consists of the system (1) with the effective matrix-fracture source term Here with We note that linearization of the imbibition equation by a constant is a special case of this new linearization by using a variable coefficient.
Remark 1 Numerical simulations (for the details see [8]) have shown that the effective matrix-fracture exchange source term (4) obtained by the new variable linearization in [8] provides a much better approximation of the exact matrix-fracture source term compared to the corresponding effective matrix-fracture source term (2) obtained previously by a constant linearization in [5].

Discretization of double porosity model I
The model ( 1)-(2) will be discretized by the cell centered finite volume method on a structured grid with the two-point flux approximation.First we present the time discretization.
Assume that we have a sequence of time steps: 0 = t 0 < t 1 < • • • < t n < • • • and denote δt n = t n+1 − t n and also I n = (t n−1 , t n ].All unknowns are supposed to be piecewise constant in time so that S(x, t) = k S k (x)χ I k (t), where S k (x) = S(x, t k ), and similarly for other variables.Implicit Euler discretization gives for t ∈ I n+1 , The source term is discretized in the following way: where we denoted for 1 ≤ k ≤ n, Obviously, If the time grid is equidistant, then we have leading to a convolution-like representation: Generally we have: where, for n > 0, Also, note that in the case n = 0 we have, Therefore, for n = 0 we have is monotone decreasing.We also introduce Let us note that in the equidistant time stepping we have I n+1 k+1 − I n k = 0 and then In the non equidistant case the terms I n+1 k+1 − I n k can have any sign, so we will introduce the assumption that the time discretization is such that This will always be the case if the time stepping is close to the equidistant one.
with D n k > 0 for k = 0, 1, . . ., n.Let us also note that We use standard finite volume discretization of the two phase system written in the phase formulation [13] (see [14] for notations): In this discretization we use phase by phase upstream choice: the value of the mobility of each phase on the edge K|L is determined by the sign of the difference of the discrete phase pressure. with where E n+1 w and E n+1 n are two subsets of E such that

Discretization of double porosity model II
Second model differs from the first one only in the matrix-fracture exchange term which takes the form: We have chosen expression for the matrix-fracture exchange term given by (18) but it is also possible to use other forms, for example (4).
We will discretize this model using the same approach as in the constant linearization model.Assume that we have a sequence of time instances 0 = t If the saturation S is constant in time on each interval (t k , t k+1 ) then P(S(x, (τ x ) −1 (u))) − P(S(x, 0)) is constant on each interval (τ k x , τ k+1 x ) and we can write x 0 P(S(x, (τ x ) −1 (u))) − P(S(x, 0)) x − u (P(S k (x)) − P(S 0 (x))) .
As before we have Note that for s ∈ (t k−1 , t k ) we have where α l m = α m (t l ), so that Therefore, we have for k ≤ n, For notational simplicity we will introduce for k ≤ n and U n k = 0 for k > n.Then we can write: We finally obtain the following scheme: Note that . where Using standard finite volume discretization of the two phase system written in the phase formulation (see [13]) we get

Conclusion
We study the double porosity models (1), ( 2) and ( 1), ( 4) of the incompressible two-phase flow in fractured porous media which are obtained after a standard periodic homogenization as ε → 0 and then linearizing the nonlinear imbibition equation in a two different ways: by replacing the nonlinear function (3) by a constant (as in [5], based on the idea of [7]), and by replacing (3) by a variable function, as in [8].The need for the linearization of the imbibition equation comes from the fact that its nonlinearity causes difficulties for numerical simulations since the source terms at each point of the domain Ω are given through solutions of the imbibition equation so the associated boundary value problem needs to be solved many times.In models (1), ( 2) and ( 1), (4) the coefficients of the system are given explicitly and for their calculation solving local problems is no more required.In this work we present discretization schemes by the cell centered finite volume method for the effective system in both cases of linearization.Special care is taken of the matrix-fracture source terms which are given in a form of a convolution.The numerical simulations of the simplified double porosity model are a part of our future research.