Two-phase flow through fractured porous media

I. I. Bogdanov, V. V. Mourzenko, J. -F. Thovert, P. M. Adler

► To cite this version:

I. I. Bogdanov, V. V. Mourzenko, J. -F. Thovert, P. M. Adler. Two-phase flow through fractured porous media. Physical Review E, 2003, 68, pp. 767-782. 10.1103/PhysRevE.68.026703. insu-03598413

HAL Id: insu-03598413

https://insu.hal.science/insu-03598413

Submitted on 16 Mar 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Two-phase flow through fractured porous media

I. I. Bogdanov,¹,² V. V. Mourzenko,¹ J.-F. Thovert,¹ and P. M. Adler ²

¹Laboratoire de Combustion et de Détonique, SP2MI, BP 30179, 86962 Futuroscope Cedex, France

²Institut de Physique du Globe de Paris, tour 24, 4 Place Jussieu, 75252 Paris Cedex 05, France

(Received 24 October 2002; revised manuscript received 9 April 2003; published 15 August 2003)

Two-phase flow in fractured porous media is investigated by means of a direct and complete numerical solution of the generalized Darcy equations in a three-dimensional discrete fracture description. The numerical model applies to arbitrary fracture network geometry, and to arbitrary distributions of permeabilities in the porous matrix and in the fractures. It is used here in order to obtain the steady-state macroscopic relative permeabilities of random fractured media. Results are presented as functions of the mean saturation and are discussed in comparison with simple models.

DOI: 10.1103/PhysRevE.68.026703

PACS number(s): 47.11.+j, 47.55.Mh, 91.45.Vz

# I. INTRODUCTION

The description of two-phase flows in fractured porous media is a challenging problem, because of the multiple scales that are involved and of the nonlinearity of the governing equations. Yortsos [1] discusses the various regimes that may take place in different situations.

Immiscible binary flows in homogeneous porous media are commonly described by generalized Darcy's equations, which are based on the concept of relative permeability (Marle [2]; Bear and Bachmat [3]; Dullien [4]). A two-dimensional (2D) equivalent of this formulation can be applied to two-phase flows through fractures (e.g., Wang and Narasimhan [5]). This description results from the upscaling of the Stokes equations that govern the fluid motion at the pore level, and it applies on a local scale, large compared to the microscopic pore scale, but small compared to the macroscopic scale on which the medium properties may vary. It assumes that the microscopic distribution of the two fluids is controlled by the capillary forces. The capillary pressure and the relative permeabilities are functions of the microstructure of the porous material, and depend on the local saturation.

Just like the classical Darcy law, generalized Darcy's equation was first introduced on an empirical basis (Muskat and Meres [6]; Muskat et al. [7]), and later received a justification in various theoretical frameworks (see Bourgeat [8]), if the classical requirement of scale separation of the homogenization theory and a few conditions relative to the fluid interfacial properties are satisfied.

Two-phase flows in heterogeneous or fractured porous media are described by the generalized Darcy laws for the two fluids, coupled by a global continuity equation, and supplemented with constitutive equations for the relative permeabilities and saturation-capillary pressure relationships, which are generally nonlinear. The difficulty of the simulation stems from this nonlinearity, from the sharp contrast of the matrix and fracture properties and from the random character of the medium geometry. Thus, most of the earlier conceptual or numerical approaches replaced this complex system by simpler idealized models.

A common approach is based on the extension of the double porosity model introduced for single-phase flow by Barenblatt and Zheltov [9] and Barenblatt et al. [10], fol-

lowed by Warren and Root [11]. The review article of Kaczemi and Gilman [12] describes many variants of such models. Panfilov [13] combines the double porosity model with a double scale expansion in order to derive the effective macroscopic phase permeabilities, and numerically solves the closure problems in a few simple two-dimensional situations.

Capillary models have also been used for direct simulation of drainage and imbibition (Fenwick and Blunt [14]; Patzek [15]; Hughes and Blunt [16]) or to evaluate effective macroscopic relative permeabilities ([17]).

The present work is based on a three-dimensional discrete description of the fracture network and of the embedding matrix. Any fracture network geometry, any type of boundary condition, and any distribution of the fracture and matrix properties can be addressed, without simplifying approximations. Therefore, this description can be used to investigate any type of flow or transport problems, as well as other processes, which may be considered in future developments, such as mechanical deformation and hydromechanical coupling, that simplified models like equivalent pipe networks are at a loss to address. The main drawback of this direct approach used to be its computational requirements—but thanks to the progress of the computers it is not a real issue anymore. The numerical tools described here are able to handle several hundreds of fractures.

Other discrete fracture numerical models have been described in the literature. Many are reviewed by Bogdanov et al. [18]. They generally do not incorporate a full 3D description of both fracture and matrix flow. For instance, the simulation package Rockflow (Kaiser et al. [19]) apparently does not include a 3D mesh of the matrix rock. The package TOUGH2 (see e.g., Wu and Pruess [20]) does not incorporate fractures as discrete elements; fractures are covered by volume elements, and a dual-porosity model is applied. The package FRACMAN distributed by Golder is apparently limited to a few fractures and to single-phase flow. The numerical model CompFlow (Unger et al. [21]) and its extension by Slough et al. [22] impose that all the fracture planes are orthogonal. In all these cases, a major limitation seems to be the lack of an appropriate 3D mesh generator. A possible exception is the model of Bastian et al. [23], coupled with the mesh generator of Schöberl [24], but we are not aware of any systematic application of this software.
The physical model relies on generalized Darcy's equations, which are assumed to apply on the intermediate scale of our description, i.e., typically on a metric scale, with nonlinear and possibly hysteretic constitutive models for the capillary function and the relative permeabilities. In addition, we use a local equilibrium hypothesis, including between the fractures and the immediately surrounding matrix, which means that the capillary pressure is continuous.

There are of course many complex, yet realistic and interesting situations where such a description is inadequate. For instance, fractures generally have self-affine geometries and scaling properties, which make an homogenized description very difficult, if not impossible, on a range of scales. Impossible fluid displacement with buoyancy effects in self-affine fractures was investigated theoretically by Schmittbuhl et al. [25] and experimentally by Auradou et al. [26]. It was shown that long range correlations induce peculiar phase distribution patterns and control many aspects of the flow process. Many other recent studies addressed two-phase flows in fractures or in porous media on a microscopic scale (e.g., Aker et al. [27]; Flekkøy et al. [28]; Kundsen and Hansen [29]; Méheust et al. [30]). Emphasis is often put on flow channelization and fingering, and many effects that cannot be easily homogenized are reported.

However, we do not address this first upscaling problem in this paper. Instead, we assume that our intermediate description scale is large enough to encompass the range of local self-affinity, so that the Darcy equations are applicable. Note, to support this, that even though the fracture surfaces are often self-affine over a very wide range, their two surfaces are generally mated on the large scale. Hence, a cutoff length exists beyond which a regular behavior of the fracture aperture and of its transport properties is recovered (see, e.g., Brown et al. [31]).

From this local description, we may proceed with a second upscaling up to the field scale, when the standard statistical homogeneity requirements are fulfilled. Most of the present calculations address this situation, in Secs. V and VI. There are also cases where a second upscaling is inapplicable, for at least three reasons. The first reason stems from the nonlinearity of the equations, which causes the apparent flow properties to become rate dependent. In the following, a criterion is introduced, which is basically a large scale capillary number, below which this effect is negligible. The two other reasons are related to structural features. They are the two classical inhibitors for any homogenization procedure, namely, large scale heterogeneity, such as a fracture network with scaling properties, and some types of flows where the boundary conditions induce significant gradients over ranges that are comparable with the typical length scales of the structure. In such situations, the numerical model can still be applied for direct simulations. Large scale heterogeneity or self-affinity are not addressed in this paper, but a detailed example of simulations with a doublet of injection-production wells is provided.

The purpose of this paper is to briefly present the methodology and the first results obtained in the determination of the two-phase flow properties of fractured porous media. It is a significant extension of our previous papers on single-

phase permeability of fracture networks (Koudina et al. [32]) and on incompressible or compressible single-phase flow in fractured porous media (Bogdanov et al. [18]; Bogdanov et al. [33]). A general description of the approach is given by Adler and Thovert [34]. As already stated, our simulation code for unsteady two-phase flow, based on a fully three-dimensional discrete fracture model, can be applied to case studies, in an arbitrary setting, or to determine the macroscopic properties of a fractured medium to be used in an upscaled description.

The paper is organized as follows. Section II provides the mathematical framework for two-phase flow in fractured porous media, including the transport and the constitutive equations that are eventually reformulated in dimensionless form. Dimensionless parameters and criteria are also introduced to quantify various physical regimes. In particular, an *a priori* criterion for the possibility of upscaling generalized Darcy's equation is devised, which is later confirmed by the numerical simulations.

Section III addresses the numerical aspects. The 3D meshing of randomly fractured media is described first. Then, the spatial and temporal discretizations of the equations and the solution algorithm are presented.

Networks of parallel fractures are considered in Sec. IV. The illustrative simple example of an array of infinite parallel fractures is treated first. Since an analytical solution can be obtained in this case, it provides a direct check of the numerical codes. A generalization for steady-state two-phase flow of the classical result of Snow [35] for single-phase flow in networks of infinite plane fractures is given. Finally, the simulation of the flow in a closed regularly compartmented reservoir is presented in details, and discussed in comparison with homogenized models.

Complex realistic situations are addressed in Sec. V. We consider here networks of random fractures in a permeable rock matrix. Detailed results are given first in Sec. V A for an illustrative case. Then, more systematic results are presented in Sec. V B; the steady-state macroscale phase relative permeabilities are determined as functions of saturation, for typical situations with percolating or nonpercolating fracture networks. The influence of the other parameters is briefly considered in Sec. V C.

The previous results are discussed in Sec. VI, in comparison with two simple models.

Finally, a few concluding remarks are gathered in Sec. VII.

## II. MATHEMATICAL FORMULATION

### A. Transport equations

Recall that the governing equations for the flow are written at some intermediate scale, small compared to the fracture extension but large compared to the typical pore size in the matrix and to the typical fracture aperture. Hence, they result from the homogenization of the microscopic Navier-Stokes equations, and the standard requirements of statistical homogeneity for this preliminary upscaling are supposed to be fulfilled.
Let the porous rock matrix have a porosity $\epsilon_m$ and a bulk permeability $K_m [L^2]$ that can vary with space. The flow in the matrix is described by a generalized Darcy law for each phase, with relative permeabilities $K_{r,i}$ ($i=w,n$). Subscripts $w$ and $n$ refer to the wetting and nonwetting fluids, respectively. The local seepage velocities $\overline{v}_i$ are given by

$$ \overline{v}_i = - \frac{K_m K_{r,i}}{\mu_i} \nabla (P_i - \rho_i gz) \quad (i=w,n), \quad (1a) $$

where $\mu_i$ is the viscosity, $\rho_i$ is the density, and $P_i$ is the pressure for fluid $i$. For concision, denote $\Phi_i$ the potential $P_i - \rho_i gz$ and $\Lambda_i = K_{r,i}/\mu_i$ the phase mobilities. Then, Eq. (1a) reads

$$ \overline{v}_i = -K_m \Lambda_i \nabla \Phi_i \quad (i=w,n). \quad (1b) $$

The fluids are considered as incompressible. Hence, two continuity equations and a global condition on the saturations $S_i$ can be written as

$$ S_n + S_w = 1, \quad (2a) $$

$$ \epsilon_m \frac{\partial S_i}{\partial t} + \nabla \cdot \overline{v}_i = 0 \quad (i=w,n). \quad (2b) $$

Additional constitutive equations are required to relate the capillary pressure $P_c = P_n - P_w$ and the relative permeabilities $K_{r,i}$ to the fluid saturations. They are discussed in Sec. II B.

Note that in the limit of an inviscid and weightless non-wetting phase ($\mu_n=0, \rho_n=0$), $P_n$ is a constant, which can be taken equal to 0 without loss of generality, and Eqs. (1) and (2) reduce to

$$ \epsilon_m \frac{\partial S_w}{\partial t} + \nabla \cdot [K_m \Lambda_w \nabla (P_c + \rho_w gz)] = 0. \quad (3) $$

This is the so-called Richard's equation (Richards [36]), which is commonly used in hydrological studies to describe unsaturated ground water flow.

Equations similar to Eqs. (1) and (2) are applied for the flow through the fractures. We assume that the hydraulic properties of a fracture can be described by an effective conductivity $\sigma [L^3]$. The in-plane flow rates $j_{s,i}$ per unit width are related to the surface pressure gradients $\nabla_s P_i$ by the two-dimensional generalized Darcy laws:

$$ j_{s,i} = - \frac{\sigma \sigma_{r,i}}{\mu_i} \nabla_s \Phi_i \quad (i=w,n), \quad (4a) $$

where $\sigma_{r,i}$ are the relative permeabilities of the fractures. If $\lambda_i$ denotes the fluid mobility in the fracture $\sigma_{r,i}/\mu_i$, Eq. (4a) reads

$$ j_{s,i} = -\sigma \lambda_i \nabla_s \Phi_i \quad (i=w,n). \quad (4b) $$

The conductivity $\sigma$ can be position and fracture dependent. For a fracture that can be viewed locally as a plane channel

of aperture $b$, filled with a porous material (e.g., gouge or breccia) with permeability $K_f$, $\sigma$ is given by

$$ \sigma = bK_f. \quad (5) $$

It is assumed that the fractures oppose a negligible resistance to flow normal to their plane. Hence, the pressures $P_i, P_c$, and the potentials $\Phi_i$ are continuous across the fractures.

Again, constitutive equations, discussed in the following subsection, are required to relate the capillary pressure and the relative permeabilities $\sigma_{r,i}$ to the fluid saturations in the fractures.

Conservation equations similar to Eq. (2) could be written for the fractures, which should include exchange terms with the surrounding matrix. In view of the finite volume scheme used for the numerical solution, it is more convenient to write a global conservation equation that accounts for both matrix and fracture flow in a control volume. Suppose that the volume $\Omega$, with boundary $\partial\Omega$, contains part of one or several fractures, denoted by $F$. By applying the divergence theorem, conservation of phase $i$ in $\Omega$ can be written as

$$ \int_{\Omega-F} \epsilon_m \frac{\partial S_i}{\partial t} dv + \int_{\Omega \cap F} \epsilon_f \frac{\partial S_i}{\partial t} dv + \int_{\partial \Omega-F} \mathbf{n} \cdot \overline{\mathbf{v}_i} ds \\ + \int_{\partial \Omega \cap F} \mathbf{n} \cdot \mathbf{j}_{s,i} dl = 0 \quad (i=w,n), \quad (6) $$

where $\mathbf{n}$ is the unit vector normal to $\partial\Omega$. The volume of the fractures in $F$ is supposed to be negligible compared to the pore volume in the matrix in $\Omega-F$. Therefore, Eq. (6) can be simplified into

$$ \int_{\Omega} \epsilon_m \frac{\partial S_i}{\partial t} dv + \int_{\partial \Omega} \mathbf{n} \cdot \overline{\mathbf{v}_i} ds + \int_{\partial \Omega \cap F} \mathbf{n} \cdot \mathbf{j}_{s,i} dl = 0 \quad (i=w,n). \quad (7) $$

Thus, the fractures introduce a singular contribution to the mass balance equation (2b).

The flow equations can be rewritten in terms of one of the potentials and of the capillary pressure, thanks to condition (2a). In the rock matrix, the sum of Eqs. (2b) for the two phases yields

$$ \begin{aligned} & \nabla \cdot [K_m (\Lambda_w + \Lambda_n) \nabla \Phi_n] \\ & = \nabla \cdot (K_m \Lambda_w \nabla P_c) - \frac{\partial}{\partial z} (K_m \Lambda_w) \Delta \rho g \mathbf{e}_z, \end{aligned} \quad (8a) $$

where $\Delta\rho = \rho_n - \rho_w$ is the density contrast and $\mathbf{e}_z$ is the vertical unit vector. This linear stationary equation relates $\Phi_n$ to the instantaneous capillary pressure field. On the other hand, Eq. (2b) for the wetting phase can be written as

$$ \epsilon_m \frac{\partial S_w}{\partial t} = \nabla \cdot [K_m \Lambda_w \nabla (\Phi_n - P_c + \Delta \rho gz)]. \quad (8b) $$

It describes the temporal evolution of $S_w$ (and thus of $P_c$ and $S_n$), as a function of the fields $\Phi_n$ and $P_c$.
The balance equations (7) over a control volume, which possibly contains fractures, can be transformed in the same way

$$ 
\begin{aligned}
& \int_{\partial \Omega} K_m (\Lambda_w + \Lambda_n) \mathbf{n} \cdot \nabla \Phi_n ds + \int_{\partial \Omega \cap F} \sigma (\lambda_w + \lambda_n) \mathbf{n} \cdot \nabla \Phi_n dl \\
& = \int_{\partial \Omega} K_m \Lambda_w (\mathbf{n} \cdot \nabla P_c - \Delta \rho g \mathbf{n} \cdot \mathbf{e}_z) ds \\
& \quad + \int_{\partial \Omega \cap F} \sigma \lambda_w (\mathbf{n} \cdot \nabla P_c - \Delta \rho g \mathbf{n} \cdot \mathbf{e}_z) dl,
\end{aligned}
\quad (9a) $$

$$ 
\begin{aligned}
\int_{\Omega} \epsilon_m \frac{\partial S_w}{\partial t} dv = & \int_{\partial \Omega} K_m \Lambda_w \mathbf{n} \cdot \nabla (\Phi_n - P_c + \Delta \rho g z) ds \\
& + \int_{\partial \Omega \cap F} \sigma \lambda_w \mathbf{n} \cdot \nabla (\Phi_n - P_c + \Delta \rho g z) dl.
\end{aligned}
\quad (9b) $$

This formulation is the starting point for the numerical scheme described in Sec. III.

## B. Constitutive equations

Constitutive equations are required for the closure of the set of transport equations in the preceding subsection. Note that their choice is partly arbitrary, and different models may be suited for various types of rock matrix or fractures. Those implemented here are among the most common, and they could easily be substituted with others to address specific situations.

First, the stress balance at the fluid interface at the microscopic scale has to be taken into account. Due to interfacial tension, a pressure jump $P_c$ takes place across the interface, which is called the capillary pressure

$$ P_c = P_n - P_w = \Phi_n - \Phi_w + \Delta \rho g z. \quad (10) $$

Experimental evidence (see, e.g., Marle [2]; Adler and Brenner [37]; Bear and Bachmat [3]; Reitsma and Kueper [38]) shows that $P_c$ is related to the saturations $S_i$. Numerous phenomenological or semiempirical formulas exist to model this relationship (see Chen et al. [39] and the references above). The most widely used is the van Genuchten equation (van Genuchten [40]), which can be written as

$$ S_w = \left[ 1 + \left( \frac{P_c}{P_0} \right)^n \right]^{\frac{(1-n)}{n}}, \quad (11) $$

where $P_0$ is a characteristic pressure, $n$ is an index. Typical values of $n$ range from 1 to 4. Estimates of $P_0$ are given below.

Note that experimental measurements of capillary pressures are generally performed with fluids at rest, when the interface freely settles at a location that minimizes its area. Viscous forces may deform the interface if the fluids are flowing, which may affect the capillary pressure-saturation relationship. However, this effect should remain reasonably

small if the capillary number $C$ that compares the magnitudes of the viscous and interfacial stresses is small enough. It is defined as

$$ C = \frac{\mu U}{\gamma}, \quad (12) $$

where $U$ is a typical fluid velocity and $\gamma$ is the interfacial tension. This number is generally very small for underground flows (typically less than $10^{-6}$).

Equation (11) was applied both in the rock matrix and in the fractures. Recall that owing to the local equilibrium hypothesis, the pressures $P_n$ and $P_w$ and the capillary pressure $P_c$ are continuous; they are equal in a fracture and in the matrix rock along its surface. However, the parameters $P_0$ and $n$ are generally different in the two domains. In particular, at the microscopic scale, the pressure jump across the interface between the two fluids is inversely proportional to the meniscus radius, which is of the order of the typical pore size. Therefore, one may expect that

$$ P_0 \propto \frac{\gamma}{\sqrt{K}}, \quad \frac{P_{0,f}}{P_{0,m}} \approx \sqrt{\frac{K_m}{K_f}} = \kappa. \quad (13) $$

Unless otherwise stated, the computations in the following correspond to $P_{0,f} = \kappa P_{0,m}$.

On the other hand, the relative permeabilities appearing in generalized Darcy's equations (1) and (4) also depend on the fluid saturations. Many models have been considered for porous media (see, e.g., Bear and Bachmat [3]; Chen et al. [39]). Again, the most widely used model for the wetting phase relative permeability was proposed by van Genuchten [40]:

$$ K_{r,w} = S_w^{1/2} \left[ 1 - \left( 1 - S_w^{n/(n-1)} \right)^{(n-1)/n} \right]^2. \quad (14) $$

Note that Eqs. (11) and (14) suppose that $S_w$ can vary over the whole range from 0 to 1. If its practical variations are limited by irreducible and maximal values $S_{wr}$ and $S_{ws}$, Eqs. (11) and (14) are generally written in terms of the effective saturation $\tilde{S}_w$:

$$ \tilde{S}_w = \frac{S_w - S_{wr}}{S_{ws} - S_{wr}}. \quad (15) $$

Residual saturations were not considered in the present simulations, but could easily be included.

The relative permeability for the nonwetting phase is also sometimes modeled according to Eq. (14), with $K_{r,w}$ and $S_w$ replaced by $K_{r,n}$ and $S_n$, respectively. However, the relative permeability curves for the two phases are generally not mirror images of one another. Thus, a different model was used here, which is discussed below.

Two-phase flows in fractures have given rise to comparatively less experimental studies than three-dimensional porous media, but a few references can be found in the literature. They are reviewed, for instance, by Persoff and Pruess [41] and Fourar [42]. The relative permeabilities $\sigma_{r,i}$ are sometimes found proportional to the saturations $S_i$, with
$\sigma_{r,w} + \sigma_{r,n} = 1$ (e.g., Romm [43]), but more generally lower values with $\sigma_{r,w} + \sigma_{r,n} \ll 1$ are observed (see Fourar and Boies [44]). This is especially true for rough walled fractures. In the present simulations, we used a simple model for $\sigma_{r,n}$:

$$ \sigma_{r,n} = S_n^q, \quad (16) $$

with the exponent $q$ equal to 2. This model was also applied for $K_{r,n}$ in the rock matrix. Morover, $\sigma_{r,w}$ was described by an equation of the type of Eq. (14). It can be noted that in most of the situations considered here, the wetting phase saturation in the fractures is very small and, therefore, $\sigma_{r,w}$ is very small and $\sigma_{r,n}$ is of the order of unity, for any reasonable choice of constitutive equations.

In summary, the capillary pressure $P_c$, the relative permeabilities for the wetting fluid $K_{r,w}$ and $\sigma_{r,w}$ and the relative permeabilities for the nonwetting fluid $K_{r,n}$ and $\sigma_{r,n}$ are described both in the matrix and in the fractures by Eqs. (11), (14), and (16), respectively.

Note that although the capillary pressure $P_c$ is identical in the fracture and in the adjacent matrix, the saturations may differ if $P_{0,f} \neq P_{0,m}$ [see Eq. (13)] or if the exponent $n$ has different values $n_f \neq n_m$. Therefore, the relative permeabilities $K_{r,i}$ and $\sigma_{r,i}$ may also be different. Sharp saturation contrasts are actually the rule, due to the small value of $\kappa$. For instance, with the constitutive parameters in Secs. V A and V B, $S_w$ in the fractures is smaller than 0.05 for any saturation smaller than 0.85 in the matrix nearby.

### C. Dimensionless formulation

Dimensionless variables, denoted by primes, can be introduced by using characteristic quantities of the problem. A possible choice is based on the typical values of the matrix porosity and permeability, $\bar{\epsilon}_m$ and $\bar{K}_m$, on the pressure $P_{0,m}$, on the wetting fluid viscosity $\mu_m$ and on a characteristic length scale $L$:

$$ \mathbf{r}' = \frac{1}{L} \mathbf{r}, \quad \nabla' = L\nabla, \quad \epsilon'_m = \frac{\epsilon_m}{\epsilon_m}, \quad (17a) $$

$$ P'_i = \frac{P_i}{P_{0,m}}, \quad \Phi'_i = \frac{\Phi_i}{P_{0,m}}, \quad \rho'_i = \frac{gL}{P_{0,m}}\rho_i, \quad (17b) $$

$$ K' = \frac{K}{\bar{K}_m}, \quad \sigma' = \frac{\sigma}{L\bar{K}_m}, \quad \Lambda'_i = \mu_w \Lambda_i, \quad (17c) $$

$$ t' = \frac{\bar{K}_m P_{0,m}}{\bar{\epsilon}_m \mu_w L^2}, \quad t, \quad \mathbf{v}'_i = \frac{\mu_w L}{\bar{K}_m P_{0,m}} \mathbf{v}_i, \quad j'_i = \frac{\mu_w}{\bar{K}_m P_{0,m}} j_i. \quad (17d) $$

Hence, the dimensionless transport equations (1), (2), (4), (8), and (9) read

$$ \overline{\mathbf{v}}'_i = -K'_m \Lambda'_i \nabla \Phi'_i \quad (i=w,n), \quad (18) $$

$$ S_n + S_w = 1, \quad (19a) $$

$$ \epsilon'_m \frac{\partial S_i}{\partial t'} + \nabla' \cdot \overline{\mathbf{v}}'_i = 0 \quad (i=w,n), \quad (19b) $$

$$ j'_{s,i} = -\sigma' \lambda'_i \nabla'_s \Phi'_i \quad (i=w,n), \quad (20) $$

$$ \begin{aligned} & \nabla' \cdot [K'_m(\Lambda'_w + \Lambda'_n)\nabla' \Phi'_n] \\ &= \nabla' \cdot (K'_m \Lambda'_w \nabla' P'_c) - \frac{\partial}{\partial z'} (K'_m \Lambda'_w) \Delta \rho' \mathbf{e}_z, \end{aligned} \quad (21a) $$

$$ \epsilon'_{m} \frac{\partial S_w}{\partial t'} = \nabla' \cdot [K'_{m} \Lambda'_{w} \nabla' (\Phi'_{n} - P'_{c} + \Delta \rho' z')], \quad (21b) $$

$$ \begin{aligned} & \int_{\partial\Omega} K'_{m}(\Lambda'_{w} + \Lambda'_{n})\mathbf{n} \cdot \nabla' \Phi'_{n} ds' \\ & \qquad + \int_{\partial\Omega \cap F} \sigma'(\lambda'_{w} + \lambda'_{n})\mathbf{n} \cdot \nabla' \Phi'_{n} dl' \\ &= \int_{\partial\Omega} K'_{m}\Lambda'_{w}(\mathbf{n} \cdot \nabla' P'_{c} - \Delta\rho'\mathbf{n} \cdot \mathbf{e}_{z})ds' \\ & \qquad + \int_{\partial\Omega \cap F} \sigma'\lambda'_{w}(\mathbf{n} \cdot \nabla' P'_{c} - \Delta\rho'\mathbf{n} \cdot \mathbf{e}_{z})dl', \end{aligned} \quad (22a) $$

$$ \begin{aligned} \int_{\Omega} \epsilon'_{m} \frac{\partial S_w}{\partial t'} dv' = & \int_{\partial\Omega} K'_{m} \Lambda'_{w} \mathbf{n} \cdot \nabla' (\Phi'_{n} - P'_{c} + \Delta\rho' z) ds' \\ & + \int_{\partial\Omega \cap F} \sigma' \lambda'_{w} \mathbf{n} \cdot \nabla' (\Phi'_{n} - P'_{c} + \Delta\rho' z) dl'. \end{aligned} \quad (22b) $$

The constitutive equations (10), (11), (13), (14), and (16) yield

$$ P'_{c} = P'_{n} - P'_{w} = \Phi'_{n} - \Phi'_{w} + \Delta\rho' z', \quad (23) $$

$$ S_w = [1 + P_c'^{n_m}]^{(1-n_m)/n_m} \quad (\text{matrix}), \quad (24a) $$

$$ S_w = \left[ 1 + \left( \frac{P_c'}{\kappa} \right)^{n_f} \right]^{(1-n_f)/n_f} \quad (\text{fractures}), \quad (24b) $$

$$ K_{r,w} = S_w^{1/2} \left[ 1 - \left( 1 - S_w^{(n_m/(n_m-1))^{(n_m-1)/n_m}} \right)^{n_f} \right] $$

$$ K_{r,n} = S_n^q \quad (\text{matrix}), \quad (25a) $$

$$ \sigma_{r,w} = S_w^{1/2} \left[ 1 - \left( 1 - S_w^{n_f/(n_f-1)} \right)^{(n_f-1)/n_f} \right]^2 $$

$$ \sigma_{r,n} = S_n^q \quad (\text{fractures}). \quad (25b) $$
## D. Governing dimensionless parameters

In view of the numerous physical parameters, many dimensionless numbers play a role in the problem. Some of them describe the fractured medium intrinsic properties, such as $\sigma'$, which compares the single-phase fracture and matrix transmissivities. Others characterize the fluid properties such as the viscosity ratio. Finally, dimensionless numbers can be introduced to compare the magnitude of the three types of forces acting on the fluids, namely, viscous forces, interfacial tension, and buoyancy.

The capillary number $C$ (12) determines whether the fluid interface is locally deformed by the flow. It can be estimated in the matrix from the typical wetting phase velocity [see Eq. (1a)]. By using Eqs. (13) and (17),

$$C_m \sim \frac{\sqrt{K_m}}{L} K_{r,w} G, \quad (26)$$

where $G = \nabla' P'$ is the magnitude of the dimensionless macroscopic pressure gradient. The ratio $\sqrt{K_m}/L$ of the microscopic to macroscopic length scales is generally of the order of $10^{-6}$ or less. The counterpart of $C_m$ for the fractures is $C_f \sim C_m/\kappa$, which is larger but still much smaller than 1. Thus, the shape of the fluid interface is determined at the pore scale by the capillary forces only, with the important consequence that the local constitutive equations (24) and (25) do not depend on the flow rate. Note that this does not rule out a possible dependence on history, for example, with hysteretic capillary function and relative permeability.

On a larger scale, however, the flow can induce changes in the saturation field. An upper bound $\Delta S_M$ of the order of magnitude of the local saturation variations between a rest equilibrium state and a stationary flow (or by extension between two stationary flows with different macroscopic driving pressure gradients) can be estimated as

$$\Delta S_M \sim L |\nabla P| \left| \frac{\partial S}{\partial P_c} \right| = G \left| \frac{\partial S}{\partial P'_c} \right|. \quad (27)$$

The derivative in Eq. (27) is at most of the order of unity. For instance, it ranges between 0.1 and 0.4 for $n_m=2$ in Eq. (24a) when $0.2 \le S \le 0.95$. Hence, the dimensionless pressure gradient $G$ is a direct measure of the possible saturation variations, and it can be regarded as a macroscopic capillary number. In particular, these variations are negligible if $G$ satisfies

$$G \ll \frac{1}{\left| \frac{\partial S}{\partial P'_c} \right|}. \quad (28)$$

More precisely, the saturation variations should have a small effect if they are small compared to both $S$ and $1-S$. Since $\Delta S_M$ in Eq. (27) is already a loose overestimate, these two conditions can be combined into the following criterion:

$$G \le G_c = \frac{\min(S, 1-S)}{\left| \frac{\partial S}{\partial P'_c} \right|}. \quad (29)$$

This defines the physical range where two-phase flow can safely be described in terms of macroscopic relative permeabilities independent of the flow rates. It will be shown in Sec. V C 2 that Eq. (29) is a conservative criterion; the macroscopic relative permeabilities are actually found fairly constant for pressure gradients as large as $(3/2)G_c$.

The contribution of gravity forces can be quantified by the Bond number, which compares the buoyancy forces to the interfacial forces, and is classically defined as

$$B_0 = \frac{g \Delta \rho L^2}{\gamma}, \quad (30)$$

where $\Delta \rho$ is the difference in density of the two fluids. Although gravity effects are not addressed in the following, it is worth noting that in view of Eq. (13) the vertical capillary pressure gradient at rest, which controls the vertical variations of saturation, is related to $B_0$ by the ratio of the microscopic to macroscopic length scale

$$\frac{\partial P'_c}{\partial z'} = \frac{\sqrt{K_m}}{L} B_0. \quad (31)$$

All the dimensionless numbers in the above describe the balance of forces in a rest or in a steady state. The transient dynamics between two different states can be described by introducing the capillary diffusion coefficient $D_c$ [see Eq. (8b)],

$$D_c = - \frac{K K_r}{\epsilon \mu} \frac{\partial P_c}{\partial S}. \quad (32)$$

In view of the numerical results obtained in the following, it seems that the limiting step for the transition from one state to another with a different saturation field is the flow of the nonwetting fluid through the matrix, at least in the investigated range of parameters. Therefore, $K$, $K_r$, $\epsilon$, and $\mu$ should be replaced in Eq. (32) by $K_m$, $K_{r,n}$, $\epsilon_m$, and $\mu_n$, respectively. It is natural to build a Péclet number $Pe_c$, using this diffusion coefficient and the typical nonwetting phase velocity $U = K_m K_{r,n} \nabla P / \mu_n$, which yields

$$Pe_c = \frac{UL}{D_c} = L \left| \frac{\nabla P}{\partial S} \right| = \frac{G}{\left| \frac{\partial P'_c}{\partial S} \right|}. \quad (33)$$

It is worth noting that $Pe_c$ is actually independent of which fluid is used as a reference, since $K_r$ and $\mu$ in $U$ and $D_c$ cancel in Eq. (33). $Pe_c$ only depends on the pressure gradient magnitude and on the capillary properties. In view of Eq. (28), this Péclet number is smaller than 1 in the range defined by criterion (29).
Furthermore, a grid Péclet number can be obtained simply as $(\delta_M/L)Pe_c$, which is even smaller. This ensures that the flow should not induce sharp saturation fronts with a thickness smaller than the grid resolution, which would be poorly accounted for by the finite volume formulation described in Sec. III B.

A typical transient time $T$ can be defined as follows:

$$T = \frac{L^2}{D_c} = \frac{\epsilon_m \mu_n L^2}{P_{0,m} K_m} \frac{1}{K_{r,n} \left| \frac{\partial P_c'}{\partial S} \right|} \quad (34a)$$

In dimensionless forms, this time constant reads [see Eq. (17d)]

$$T' = \frac{\mu_n}{\mu_w} \frac{1}{K_{r,n} \left| \frac{\partial P_c'}{\partial S} \right|} \quad (34b)$$

In the range of parameters considered in the following, $T'$ varies from 1 to 100. In dimensional terms, this roughly corresponds from one day to one year. Such long transition times can have important consequences for the interpretation of well tests and for the exploitation of oil reservoirs.

### III. NUMERICAL MODEL

#### A. Three-dimensional meshing of fractured porous media

The first step of the numerical solution is to discretize the fracture network and then the porous medium surrounding the fractures in a consistent way. The geometry of the mesh to be generated is constrained by many randomly located fractures. In addition, meshes should be routinely built for large statistical sets of stochastically generated samples. Therefore, a very robust and fully automated meshing algorithm is required.

A literature review was done by Bogdanov et al. [18]. The meshing technique was also extensively described in this paper; therefore, it is briefly schematized here. The fracture network is triangulated first, as described by Koudina et al. [32]. Then the space between the fractures is paved by an unstructured boundary-constrained tetrahedral mesh, according to an advancing front technique. Typical performances and computational requirements are discussed by Bogdanov et al. [18]. The grid resolution can be quantified by the typical size $\delta_M$ of the surface and volume elements.

Three-dimensional views of two triangulated fractured media are shown in Fig. 1. In both cases, the fractures are hexagonal, with circumscribed radius $R$, and the cell size is $L=4R$. Figures 1(a) and 1(b) are the samples used in Sec. VA. The cell contains $N_{fr}=16$ fractures, $\delta_M=R/3$, and the mesh contains about 1800 node points, 22 000 triangles, and 11 000 tetrahedra. Figure 1(c) is one of the samples used in Sec. VB, with $N_{fr}=32$. The discretization is slightly finer ($\delta_M=R/4$). The mesh contains about 3700 node points, 46 000 triangles, and 23 000 tetrahedra.

<figure><img src="image_2.png" /><figcaption>FIG. 1. The network of $N_{fr}=16$ fractures in the sample used for the simulations in Figs. 12 and 13 (a). The saturation maps in Fig. 12 correspond to the horizontal marked plane $\Pi$. The three-dimensional meshes of the same fractured medium (b) and of another sample with $N_{fr}=32$ (c). Distances are normalized by the fracture-circumscribed radius $R$. Both samples are spatially periodic, with cell size $L=4R$. The tetrahedral volume elements in the cubic unit cell $-2 \le x, y, z \le +2$ are displayed. The protuding fractures in (b,c) sit astride the boundaries with the neighboring cells. Several periodic replicas are shown for some of them. For the sake of clarity, the edges and the intersection lines of the fractures have been thickened.</figcaption></figure>

The influence of the discretization parameter $\delta_M$ on the flow calculation results is illustrated and quantified in Sec. VA.

#### B. Spatial discretization of the equations

We describe in this subsection the spatial discretization of Eqs. (22). Time discretization is addressed in the following subsection, as part of the description of the solution algorithm.

The rock matrix is represented by tetrahedral volume elements and the fractures are represented by triangular surface elements. The transport coefficients $K_m$ and $\sigma$, as well as the porosity $\epsilon_m$, are considered as uniform over these elements.

The nonwetting phase potential $\Phi_n$ and the capillary pressure $P_c$ are evaluated at the mesh points located at the vertices of the tetrahedra and triangles. Since the fractures are viewed as vanishingly thin, empty, or very permeable layers,
<figure><img src="image_3.png" /><figcaption>FIG. 2. Control volume $\Omega$ for the mass balance equations at a node; in this illustration, the node is on a fracture.</figcaption></figure>

there is no pressure jump between two points facing each other on the two opposite sides of a fracture. Thus, single values of $\Phi_n$ and $P_c$ can be used in the numerical formulation per vertex of fracture element.

A finite volume formulation of the problem is obtained by applying the balance equations (22) to control volumes $\Omega$ surrounding each of the mesh points, as shown in Fig. 2.

For the evaluation of the surface integrals in Eq. (22), the gradients $\nabla' \Phi'_n$ and $\nabla' P'_c$ are considered as constant over each mesh element and, therefore, they are linear functions of $\Phi'_n$ and $P'_c$ at the element vertices. The mobilities $\Lambda'_i$ and $\lambda'_i$ are also supposed to be piecewise constant, per tetrahedron or triangle. For each mesh volume or surface element, a mean saturation $\langle S_w \rangle$ is defined as the volume average of the saturations $S_w$ evaluated at the element vertices. Then, the mobilities in the elements are deduced from the mean saturation $\langle S_w \rangle$. For instance, Eq. (22a) can be written as a set of linear equations relating the values of $\Phi'_n$ and $\Psi' = P'_c - \Delta\rho'z'$ at the node points

$$ \mathbf{B} \cdot \mathbf{\Phi}'_n - \mathbf{A} \cdot \mathbf{\Psi}' = 0, \quad (35) $$

where the matrices $\mathbf{B}$ and $\mathbf{A}$ depend on the absolute and relative permeabilities, in addition to the mesh geometry.

The volume integral in Eq. (22b) can be evaluated in two different ways. On one hand, one may use in each tetrahedron the value of $S_w$ deduced from the capillary pressure at the central grid point. Alternatively, one may use in the tetrahedra the mean saturation $\langle S_w \rangle$ defined above, which also depends on the capillary pressures in the neighboring mesh points. A similar treatment of the time derivative, described in the following subsection, applies in both cases.

### C. Time discretization and solution algorithm

The strong nonlinearity of the coefficients in Eqs. (22) requests an implicit time formulation. The one used here is an extension of the modified Picard scheme described by Celia et al. [45] for the solution of Richards equation. The main idea of the method is the linearization of the time derivative in a mass conservative form.

The derivative $\partial S_w / \partial t'$ is approximated by the first-order difference over a time step $\delta t'$,

$$ \frac{\partial S_w}{\partial t'} = \frac{1}{\delta t'} (S_w^p - S_w^{p-1}), \quad (36) $$

where the superscript $p$ refers to the successive times $t'_0 + p \delta t'$. $S_w^p$ is obtained as the asymptotic value of successive estimations $S_w^{p,k}$

$$ S_w^p = \lim_{k \to \infty} S_w^{p,k}. \quad (37) $$

Finally, $S_w^{p,k}$ results from the iterative calculation of

$$ S_w^{p,k+1} = S_w^{p,k} + \left( \frac{\partial S_w}{\partial P'_c} \right)^{p,k} (P_c'^{p,k+1} - P_c'^{p,k}). \quad (38) $$

Hence, the estimate of $\partial S_w / \partial t$ at iteration $k+1$ is

$$ \left( \frac{\partial S_w}{\partial t'} \right)^{p,k+1} = \frac{S_w^{p,k+1} - S_w^{p-1}}{\delta t'} = \frac{S_w^{p,k+1} - S_w^{p,k}}{\delta t'} + \frac{S_w^{p,k} - S_w^{p-1}}{\delta t'} \\ = \frac{1}{\delta t'} \left( \frac{\partial S_w}{\partial P'_c} \right)^{p,k} y'^{p,k+1} + \frac{S_w^{p,k} - S_w^{p-1}}{\delta t'}, \quad (39) $$

with $y'^{p,k+1} = P_c'^{p,k+1} - P_c'^{p,k}$. Therefore, the volume integral in Eq. (22b) can be evaluated at iteration $k+1$ as

$$ \left( \int_{\Omega} \epsilon'_{m} \frac{\partial S_w}{\partial t'} dv \right)^{p,k+1} = \int_{\Omega} \epsilon'_{m} \left[ \frac{1}{\delta t'} \left( \frac{\partial S_w}{\partial P'_c} \right)^{p,k} y'^{p,k+1} + \frac{S_w^{p,k} - S_w^{p-1}}{\delta t'} \right] dv. \quad (40) $$

The discretized equation (22b) for the capillary pressure at current iteration $k+1$ given the potential field can be written in vectorial form as

$$ \frac{1}{\delta t'} \mathbf{C}^{p,k} \cdot \mathbf{y}'^{p,k+1} + \frac{1}{\delta t'} \mathbf{\Omega} \cdot (\mathbf{S}_w^{p,k} - \mathbf{S}_w^{p-1}) = \mathbf{A}^{p,k} \cdot (\mathbf{\Phi}_n'^{p,k} - \mathbf{\Psi}'^{p,k+1}). \quad (41) $$

The matrices $\mathbf{C}$ and $\mathbf{\Omega}$ are diagonal if $S_w$ in the volume integral in Eq. (22b) is evaluated from the capillary pressure at the central grid point only (see end of Sec. III B).

Since $\mathbf{\Psi}'^{p,k+1} = \mathbf{\Psi}'^{p,k} + \mathbf{y}'^{p,k+1}$, Eqs. (35) and (41) yield finally the following set of equations, which summarizes the whole algorithm

$$ \mathbf{B}^{p,k} \cdot \mathbf{\Phi}_n'^{p,k} = \mathbf{A}^{p,k} \cdot \mathbf{\Psi}'^{p,k} \quad (42a) $$
$$ \left( \frac{1}{\delta t'} C^{p,k} + A^{p,k} \right) \cdot \mathbf{y}'^{p,k+1} = A^{p,k} (\Phi_n^{p,k} - \Psi^{p,k}) - \frac{1}{\delta t'} \mathbf{\Omega} \cdot (S_w^{p,k} - S_w^{p-1}). \quad (42b) $$

In a first implementation, the global mass conservation equation (35) was solved once per time step for the potential field $\Phi'_n$, with the phase distribution quantified by the capillary pressure $\Psi'$ at the beginning of the time step; then Eq. (42b) was solved iteratively until $\mathbf{y}'^{p,k+1}$ converges toward 0, to determine the change in saturation occurring during $\delta t$. However, it appeared that reevaluating the matrices A and B and the field $\Phi'_n$ after each updating of $\Psi'$ improves the numerical performances. Thus, Eqs. (42a) and (42b) are solved alternatively until convergence of $\mathbf{y}'^{p,k+1}$ toward 0, which makes the solution scheme fully implicit. Both elementary problems (42a) for $\Phi'_n$ and (42b) for $\mathbf{y}'^{p,k+1}$ are solved by use of a CGSTAB conjugate gradient algorithm.

Note that in the modified Picard formulation, the solution for the capillary pressure at time step $p$ is searched for in terms of its increment with respect to time step $p-1$, which is the sum of the successive $\mathbf{y}'^{p,k+1}$. Hence, the magnitude of the right-hand side in Eq. (42b) can be used as a natural stopping criterion for the Picard iterative loop (see Celia et al. [45]).

As already mentioned, the right-hand side of Eq. (42b) should decrease along the $k$ iterations. It sometimes happens, however, that its norm increases from stage $k$ to stage $k+1$. In such cases, a backtracking technique can be used (see Press et al. [46]), i.e., $\Psi'^{p,k}$ is incremented by a fraction $\lambda$ of $\mathbf{y}'^{p,k+1}$, with $\lambda_{min} \le \lambda < 1$, in order to obtain $\Psi'^{p,k+1}$. Alternatively, a Newton step can be applied to determine from the fields at stage $k$ the increment of $\Psi'^{p,k}$ that yields the smallest norm of the right-hand side of Eq. (42b) at stage $k+1$. In practice, maximum efficiency is obtained with a combination of the two techniques.

The value of $\delta t'$ is automatically recomputed at each time step. Its choice is based on a desired maximum saturation increment $\Delta S$. The new step $\delta t'_{p+1}$ is deduced from the maximum of the instantaneous local time derivative $\partial S_w / \partial t'$. In view of Eq. (36), this derivative is exactly the rate of variation over the previous step. Hence,

$$ \delta t'_{p+1} = \frac{\Delta S}{\max|\partial S_w / \partial t'|} = \frac{\Delta S \delta t'_{p}}{\max|S_w^p - S_w^{p-1}|} \quad (43) $$

Other choices are possible, such as a maximum difference of total flow rate through the cell over the time step. This latter criterion might be more relevant for the computations of steady-state global medium properties, especially for large contrasts between the properties of matrix and fractures.

<figure><img src="image_4.png" /></figure>

a

<figure><img src="image_5.png" /><figcaption>FIG. 3. Homogeneous rock matrix containing an array of infinite parallel fractures (a) and its relative permeabilities (b). The solid lines are the curves for the matrix and for the fractures, which are identical, with $n_m = n_f = q = 2$. The other lines correspond to Eq. (47) for $\sigma' = 1/4$, $\kappa = 10^{-3/2} (\dots\dots\dots)$ and $\sigma' = 25$, $\kappa = 10^{-5/2} (\dots\dots\dots)$. The symbols are results of numerical simulations for $\bar{K}_{r,n}$ (□) and $\bar{K}_{r,w}$ (○).</figcaption></figure>

## IV. REGULAR FRACTURE NETWORKS

### A. Array of infinite parallel plane fractures

Consider the situation sketched in Fig. 3(a) of a porous matrix that contains an array of infinite parallel plane fractures, with a spacing $L$ that is taken as the unit length [see Eq. (17)]. A flow is imposed by a pressure gradient parallel to the fracture planes. The matrix permeability is $K_m$, the fracture conductivity is $\sigma$ and its aperture $b \ll L$. The global absolute permeability of this fractured medium, in the direction of the fractures, is

$$ \bar{K} = \frac{(L-b)K_m + \sigma}{L} \approx K_m \left( 1 + \frac{\sigma}{LK_m} \right). \quad (44) $$

The capillary pressure is supposed to be initially uniform,
equal to $P_c$, and it obviously remains so for a stationary flow parallel to the fractures. The wetting fluid saturations in the matrix and in the fractures, $S_{wm}$ and $S_{wf}$, are different because of different capillary properties. The mean global saturation is dominated by the matrix saturation

$$ \langle S_w \rangle = \frac{(L-b) \epsilon_m S_{wm} + b \epsilon_f S_{wf}}{(L-b) \epsilon_m + b \epsilon_f} \approx S_{wm}. \quad (45) $$

The global relative permability for fluid i, $\bar{K}_{r,i}$, is easily obtained from

$$ L\bar{K}\bar{K}_{r,i} = (L-b)K_m K_{r,i} + \sigma\sigma_{r,i}, \quad (i=n,w). \quad (46) $$

This yields finally, in dimensionless terms,

$$ \bar{K}' = 1 + \sigma', \quad \bar{K}_{r,i} = \frac{K_{r,i} + \sigma' \sigma_{r,i}}{1 + \sigma'} \quad (i=n,w). \quad (47) $$

Results (44)–(47) are also valid if the fractures are not evenly spaced, with L equal to the mean spacing, i.e., to the inverse of the volumetric fracture area.

These formulas are compared to the results of the numerical code in Fig. 3. A perfect agreement is obtained. The influence of the fractures is much stronger on $\bar{K}_{r,n}$ than on $\bar{K}_{r,w}$, since they offer a preferential path for the nonwetting fluid, due to the small values of $\kappa$.

## B. Multiple families of parallel plane fractures

In order to generalize Eq. (47) to media containing several families of fractures, it is convenient to rewrite it in a slightly different form.

Note first that $\sigma'$ is simply the ratio of the absolute fracture transmissivity $\sigma$ to the absolute transmissivity $LK_m$ of the slab of matrix between two fractures. For a fractured medium in capillary equilibrium, it is natural to define the corresponding ratio $\sigma'_i$ for each of the phases

$$ \sigma'_i = \frac{\sigma \sigma_{r,i}}{LK_m K_{r,i}} = \frac{\sigma_{r,i}}{K_{r,i}} \sigma' \quad (i=n,w). \quad (48) $$

The prefactor of $\sigma'$ in Eq. (48) only depends on the capillary pressure. Then, Eq. (47) can be reformulated in the nicely symmetric form

$$ \bar{K}' = 1 + \sigma', \quad \bar{K}'_i = K_{r,i}(1 + \sigma'_i) \quad (i=n,w), \quad (49) $$

where $\bar{K}'_i$ is the global effective permeability for phase i, $\bar{K}'_i = \bar{K}' \bar{K}_{r,i}$.

Suppose now that the medium contains N families of fractures, with conductivities $\sigma_p$, spacings or reciprocal volumetric areas $L_p$ and normal vectors $n_p$ ($p=1$ to $N$), with $\sigma'_p = \sigma_p / L_p K_m$. When $S_w=1$ and a macroscopic pressure gradient $\nabla P$ is applied, the seepage velocity in the matrix is $\bar{v} = -1/\mu K_m \nabla P$ and the flow rate per unit width in the fractures is $j_p = -1/\mu \sigma_p (I - n_p n_p) \cdot \nabla P$. Thus, the absolute permeability of the fractured medium is given by

<figure><img src="image_6.png" /><figcaption>FIG. 4. A sugar-box reservoir divided into 5×7 square blocks by an array of fractures. An injection well (left) and a production well (right) are located at the centers of two blocks. Wetting fluid is injected with a flow rate $Q_w$, and the nonwetting fluid fraction in the produced fluid is $F_n$.</figcaption></figure>

$$ \bar{K}' = 1 + \sum_{p=1}^{N} \sigma'_{p} (I - n_{p} n_{p}). \quad (50a) $$

Except for the unit term in the right-hand side, which stands for the contribution of the matrix, this is the classical result of Snow [35].

A similar expression can be obtained when considering each of the phases $i=n$ or $w$, by replacing the matrix and fracture absolute permeabilities by their product with the corresponding relative permeabilities. This yields the global effective permeability tensor for the fluid i

$$ \bar{K}'_i = K_{r,i} \left( I + \sum_{p=1}^{N} \sigma'_{p,i} (I - n_p n_p) \right) \quad (i=n,w). \quad (50b) $$

In the particular case of a continuous and isotropic orientation distribution of identical fractures, with total volumetric area $1/L$, the tensors $\bar{K}'$ and $\bar{K}'_i$ are spherical, and (50) reduces to

$$ \bar{K}' = 1 + \frac{2}{3} \sigma', \quad \bar{K}'_i = K_{r,i} \left(1 + \frac{2}{3} \sigma'_i\right). \quad (51) $$

Note that both Eqs. (50a) and (50b) result from superposition principles. Hence, Eq. (50b) is valid only as long as the flow equations for each of the two fluids are linear. This requirement is fullfilled in a stationary flow in capillary equilibrium, but it is generally not satisfied in transient flows.

## C. Sugar-box reservoir

We consider here the two-dimensional problem sketched in Fig. 4. A rectangular reservoir with impermeable boundaries is divided into 5×7 square blocks with permeability $K_m$ and size $L^2$. L is taken as the unit length. The blocks are separated by an array of fractures with uniform properties $\sigma' = 1$ and $\kappa = 10^{-3/2}$. The rheological parameters are $\mu_n = 10\mu_w$ and $n_m = n_f = q=2$. The reservoir is initially at rest,
with a uniform capillary pressure $P'_c = 10$, i.e., $S_w \approx 0.10$ in the matrix and $S_w \approx 0.0032$ in the fractures.

Wetting fluid is injected with a constant flow rate $Q'_w = 1$ through a well (on the left in Fig. 4), and fluid is produced at another well (on the right in Fig. 4), which contains a fraction $F_n$ of nonwetting fluid. Both wells are located at the centers of matrix blocks. The nonwetting phase pressure at the production well is taken as the pressure reference, i.e., $P_n^{out} = 0$.

Although the problem is two-dimensional, the calculations are performed in three dimensions, in a $5 \times 7 \times 1$ L³ volume, with periodicity conditions along the z axis. The blocks are subdivided into smaller cubes with size L/5, which contain 24 tetrahedral volume elements; hence, the mesh is periodic with the nodes arranged on a face-centered cubic lattice. A dimensionless time $t'$ is defined as the ratio of the injected volume to the total volume V of the reservoir

$$t' = \frac{Q_w t}{V}. \qquad (52)$$

Saturation maps in the rock matrix are displayed for various times in Fig. 5. The progression of a relatively steep saturation front is clearly visible. Downstream of the front, the initial saturation is undisturbed ($S_w \approx 0.1$). Then, saturation rises up to about 0.7 over a distance about one block size L.

The front separates two regions where the mobility of the wetting fluid in the matrix is much smaller (downstream) or much larger (upstream) than the mobility of the nonwetting fluid. In the initial conditions that prevail downstream of the front, $K_{r,w} \approx 7.7 \times 10^{-6}$, whereas $K_{r,n} \approx 0.81$; hence, $\Lambda_w < 10^{-5} \Lambda_n$. Upstream of the front, $S_w > 0.7$, $K_{r,w} > 0.072$, and $K_{r,n} < 0.086$, which yields $\Lambda_w > 10 \Lambda_n$.

The flow globally consists of the displacement of the nonwetting fluid by the wetting fluid, as illustrated by the velocity maps for both fluids in Fig. 6, when the saturation front is about halfway through the reservoir. The wetting fluid flow rate is imposed at the injection well. Its streamlines diverge from the well towards the front, whereas the nonwetting phase velocity is very small. Downstream of the saturation front, the nonwetting fluid streamlines converge from the front towards the production well. The wetting fluid is nearly immobile, and the produced fluid is almost only the nonwetting fluid. This pattern lasts until the saturation front reaches the production well.

However, not the whole nonwetting fluid is displaced, since its mobility in the matrix becomes smaller than the mobility of the wetting fluid when $S_w$ is larger than 0.7. Thus, a significant fraction of the nonwetting fluid is left behind the front in the bulk of the matrix blocks. By capillary diffusion, it slowly reaches the surrounding fractures, through which it eventually flows to the production well. This entrapped nonwetting fluid is clearly visible in the saturation maps of Fig. 5, and in the mean saturation profiles in Fig. 7.

When the front reaches the production well, the mobility of the wetting fluid in the matrix around the well suddenly increases, and it exceeds that of the nonwetting fluid when

<figure><img src="image_7.png" /><figcaption>FIG. 5. Wetting phase saturation at various times (top to bottom). The left column corresponds to the sugar-box reservoir sketched in Fig. 4. The right column corresponds to an unfractured reservoir with uniform equivalent properties; the broken lines are shown only to make the comparison easier with the left column. In both cases, the walls are impermeable.</figcaption></figure>

$S_w \approx 0.7$. The proportion of produced wetting fluid rises dramatically (see Fig. 8). Simultaneously, the overall pressure drop strongly increases (see Fig. 9), since the increase of $\Lambda_w$ does not compensate the decrease of $\Lambda_n$.

It should be noted that the fractures are never a preferential path for the wetting fluid. In the initial conditions ($P'_c = 10, S_{w,m} \approx 0.10$), the transmissivity ratios $\sigma'_w$ [see Eq. (48)] is about $10^{-7}$; behind the front ($P'_c \approx 1.0, S_{w,m} \approx 0.71$), it is of the order of $10^{-6}$, and in a late stage not reached in this simulation ($P'_c = 0.2, S_{w,m} \approx 0.98$), it is still smaller than $10^{-4}$. Nevertheless, larger velocities are observed in Fig. 6 in the matrix near to the fractures, since $S_w$ (and thus $K_{r,w}$) is larger there than in the bulk of the blocks, due to the depletion of nonwetting fluid that migrates towards the fractures.

Conversely, the matrix and the fractures initially have contributions of the same order of magnitude to the nonwetting fluid flow, with $\sigma'_n \approx 1.2$ when $S_{w,m} \approx 0.10$, and the fractures become dominant as soon as the saturation in the matrix
<figure><img src="image_8.png" /><figcaption>FIG. 6. The velocity fields in the sugar-box matrix for the wetting (a) and nonwetting (b) fluids at time $t'=0.345$ in the simulations of Fig. 5. A few large vectors near to the wells are not plotted. The dotted lines correspond to the fractures.</figcaption></figure>

trix increases (e.g., $\sigma'_n > 10$ when $S_{w,m} > 0.7$).

Figure 10 shows in more details the flow field of both phases at various times in the block located at the middle of the reservoir. The saturation front enters the block from the left in the first frame ($t'=0.27$). Then, it progresses until it roughly covers the block at $t'=0.345$. The flows of the two fluids are roughly pistonlike upstream and downstream of the front. Within the front, the nonwetting fluid flow starts diverging towards the surrounding fractures. Finally, when the front has left the block ($t' \ge 0.6$), a quasisteady pattern is observed, where the nonwetting fluid diffuses from the bulk of the block to the fractures. Note that the velocities are much smaller than in the earlier stages, and that they have been rescaled in the figure.

For comparison, another simulation was performed by replacing the fractured rock by an homogeneous medium with uniform equivalent properties. Effective permeabilities obtained from Eq. (50b) are applied for each of the phases, and the capillary pressure-saturation relationship of the matrix is

<figure><img src="image_9.png" /><figcaption>FIG. 7. Saturation profiles, in average over the y sections, at various times that correspond to the maps in Fig. 5. The abscissae are in block size units. The data correspond to the full calculations (—) and to the simplified problem in a homogeneous equivalent medium (----).</figcaption></figure>

used, according to Eq. (45). The saturation maps obtained by this homogenized description are shown in the right column in Fig. 5, in comparison with the results of the full calculations.

<figure><img src="image_10.png" /><figcaption>FIG. 8. The fraction $F_n$ of nonwetting fluid in the produced fluid, as a function of the dimensionless time $t'$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 9. Overall pressure drop $\Delta P'_n$ between the injection and production wells, as a function of the dimensionless time $t'$.</figcaption></figure>

Although Eq. (50b) is a priori not valid in this situation, since the saturation varies and the flow equations are really nonlinear, the two sets of results are in good qualitative agreement on the global scale.

However, significant differences remain even after long times. In the late stages in Fig. 5, the nonwetting fluid is still present in the blocks in the upstream part of the reservoir, whereas this area is fully saturated with wetting fluid in the results of the homogenized model. This is also apparent in the mean saturation profiles in Fig. 7. This is due to the underestimation in the homogenized model of the transfer time of the nonwetting fluid from the bulk of the blocks to the fractures.

In the late stages, the rate of nonwetting fluid transfer from the matrix blocks to the fractures can be estimated by the following simple argument. Let $\langle P_n \rangle_B$ and $\langle P_n \rangle_F$ denote the nonwetting fluid pressures, averaged over a matrix block and over the fractures surrounding it, respectively; moreover, $\langle S_w \rangle_B$ denotes the mean saturation in the block, and $\langle \Lambda_n \rangle_B$ the mean nonwetting fluid mobility. The typical pressure gradient is $2(\langle P_n \rangle_B - \langle P_n \rangle_F)/L$, and the block perimeter is $4L$. Therefore, the transfer rate to the fractures is of the order of $8K_m\langle \Lambda_n \rangle_B(\langle P_n \rangle_B - \langle P_n \rangle_F)$,

$$ \frac{d\langle S_w \rangle_B}{dt'} \approx 8\langle \Lambda_n \rangle_B(\langle P_n \rangle_B - \langle P_n \rangle_F). \quad (53a) $$

This is a typical exchange term in dual-porosity models for compressible single-phase flows (Warren and Root [11]). It can be generalized in a more symmetric form

$$ \frac{d\langle S_w \rangle_B}{dt'} \approx 8 \frac{\langle \Lambda_w \rangle_B \langle \Lambda_n \rangle_B}{\langle \Lambda_w \rangle_B + \langle \Lambda_n \rangle_B} (\langle P_c \rangle_B - \langle P_c \rangle_F). \quad (53b) $$

When $S_w$ is large, Eq. (53b) reduces to Eq. (53a), since $\Lambda_n \ll \Lambda_w$ and the gradient of $P_w$ is much smaller than the gradient of $P_n$, but Eq. (53b) can also apply in the opposite limit when the wetting fluid mobility is small.

<figure><img src="image_12.png" /><figcaption>FIG. 10. Wetting (left column) and nonwetting (right column) velocity fields in the matrix block at the center of the reservoir, at various times.</figcaption></figure>

The time derivative of the nonwetting fluid content of the central block in the reservoir, as obtained from the full numerical calculations, is compared to predictions (53) in Fig. 11. A good agreement, within about 15%, is observed for $t' \ge 0.6$, which corresponds to the time when the saturation front has crossed and left the block (see Fig. 5), and when the nonwetting fluid diverging flow pattern (see Fig. 10) has established. At earlier times, when $\Lambda_n$ is not much smaller than $\Lambda_w$, the saturation variations mostly result from the displacement mechanism, which can obviously not be described by Eq. (53), but they are reasonably well accounted for by the homogeneous medium model.

Consequently, a fairly good description of the flow both
<figure><img src="image_13.png" /><figcaption>FIG. 11. The time derivative $d\langle S_w \rangle_B / dt'$ of the mean saturation in the matrix block at the center of the reservoir as a function of the dimensionless time $t'$. The data correspond to the numerical simulations (—) and to the predictions of Eq. (53a) (---) and Eq. (53b) (---).</figcaption></figure>

downstream and upstream of the saturation front could be achieved by a dual-porosity model, with the capillary pressure-saturation relationship of the matrix in the primary porosity, the fluid effective permeabilities obtained from Eq. (50b) in the secondary porosity, and a coupling equation of the form of Eq. (53b). Such a model is presented by Panfilov [13].

However, it should be remembered that the present test case is the ideal sugar-box model, where such a description can indeed be expected to be successfully and practically applicable. In the general case of complex random networks of finite fractures, there is no equivalent of result (50b). The determination of the steady-state global relative permeabilities is the object of Sec. V. The modelization of the exchange term is another difficult problem, due to the randomness of the matrix domains shape and size, and it should probably involve a broad range of time scales. Finally, in view of Eq. (53b), such a model cannot account for the spatial variations of saturation that exist on the fracture scale, even in a steady-state flow, as demonstrated in Sec. V A.

## V. RANDOMLY FRACTURED POROUS MEDIA

We address in this section complex situations, where the rock matrix contains a network of randomly located fractures. The first subsection provides detailed illustrative results in a single random realization, including the evolution of the saturation maps, starting from various initial phase distributions. Then, more systematic results are presented in Sec. V B. The relative permeabilities of randomly fractured rocks are computed for two fracture densities that correspond to nonpercolating and percolating networks, respectively. In view of the large number of parameters, only the mean saturation is varied, and all other quantities are kept constant and equal to typical values. The influence of these parameters is briefly considered in Sec. V C.

Throughout this section, the fractured medium is supposed to be macroscopically homogeneous. Therefore, for the purpose of determining macroscale relative permeabilities, an unbounded medium is represented by the periodic juxtaposition of identical cubic unit cells. Hence, the geometry, the local fluxes, the saturation fields, and the pressure gradient are spatially periodic. The flow is induced by a large scale pressure gradient. The global phase saturations are given as initial conditions and are conserved along time.

### A. Illustrative examples

A detailed set of results is presented in this section, relative to an homogeneous matrix rock containing the fracture network shown in Figs. 1(a) and 1(b). The fractures are plane regular hexagons, with a circumscribed radius $R$, which is taken as the length unit $L$ in Eq. (17), and they have a constant permeability $\sigma' = 1$ with $\kappa = 10^{-3/2}$. The cell size is $(4R)^3$, and $\delta_M = R/3$. The network does not percolate. The constitutive parameters are set as described in Sec. II B, with $n_m = n_f = q = 2$. The fluid densities are equal and their viscosity ratio is $\mu_n / \mu_w = 10$.

All the frames in Fig. 12 are wetting phase saturation maps $S_w$ in the matrix in the horizontal section $\Pi$ in Fig. 1(a). The three columns of Fig. 12 show the evolution of the saturation field when a macroscopic flow is induced by a pressure gradient $\|\nabla P\| = P_{0,m}/R$, starting from the three different initial saturation fields in the top row, with identical global saturation $\bar{S}_w = 0.371$.

The leftmost initial state corresponds to rest equilibrium, with uniform capillary pressure. The two others are arbitrary initial phase distributions, with bands normal or parallel to the applied pressure gradient $\nabla P$.

In all cases, an identical steady regime is reached, where saturation is not uniform. $S_w$ ranges from about 0.32 to 0.43, and it is different on the inlet and outlet sides of the fractures. The disturbances in the saturation field introduced by the presence of the fractures during a steady flow with respect to the rest state are due to the different capillary properties of the fractures and rock matrix [see Eq. (13)]. They are observed in all our simulations, and increase with the fracture permeability $\sigma'$ and with the mean flow rate (or pressure gradient). This effect will be specifically investigated in a future work.

The evolutions of the global mean seepage velocities $\bar{v}_w'$ and $\bar{v}_n'$ are plotted for the three initial phase distributions in Fig. 13. $\bar{v}_w'$ corresponds to the instantaneous value of the global relative permeability $\bar{K}_{r,w}$ defined in Sec. V B, and $\bar{v}_n'$ corresponds to $\mu_w / \mu_n \bar{K}_{r,n}$. Identical phase flow rates are reached in all cases, although the convergence is not monotonous. The dimensionless time required to reach the limit is of the order of 10, which is consistent with the time constant $T' \approx 3$ predicted by Eq. (34b).

Before we proceed with systematic calculations in the following section, it may be the right place to discuss the effects
<figure><img src="image_14.png" /><figcaption>FIG. 12. Wetting phase saturation at various times (top to bottom), in the plane II marked in Fig. 1(a). The mean flow is oriented from the left to the right of the figure. The mean saturation is always $\bar{S}_w = 0.371$. Each column corresponds to the evolution of the saturation, which starts from the initial condition shown at the top of the column.</figcaption></figure>

of the grid resolution on the results of the simulations. This is done here by considering a second example, which we often used in the past as an illustrative benchmark (see Koudina et al. [32]). The sample has a size $L=3R$ and it contains $N_f=10$ fractures. All the transport and constitutive parameters are set exactly as in the previous example. The fracture network percolates in the x direction only. The spanning cluster contains six fractures, and four fractures are isolated. Hence, a variety of situations can be tested by setting the macroscopic pressure gradient along the x, y, and z axes.

Three meshes were built, with different discretization parameters $\delta_M=R/3$, $R/5$, and $R/8$. Sections through these meshes are shown in the left column of Fig. 14. As in the previous case, calculations were run from an arbitrary initial phase distribution, until a steady state was reached, where the phase permeabilities were measured. The results are given in Table I.

Note first that nearly identical effects are obtained in the three directions, regardless of the percolation status. The computed absolute permeability $\bar{K}$ varies very little; it is larger by 1% ($\delta_M/R=5$) and by 3% ($\delta_M/R=3$) than for $\delta_M/R=8$. The wetting phase flow rate, i.e., $\bar{K}\bar{K}_{r,w}$, is nearly insensitive to the grid resolution, which is natural since it is essentially the same as for intact matrix, as discussed in Sec. V B. The largest differences are observed for the nonwetting phase flow rate. $\bar{K}\bar{K}_{r,n}$ is larger by 3% ($\delta_M/R=5$) and by 8% ($\delta_M/R=3$) than for $\delta_M/R=8$. The error scales roughly as a quadratic function $0.7(\delta_M/R)^2$. Hence, an overestimate of about 5% is a fair guess of the error associated with $\delta_M = R/4$, which is used throughout the rest of the paper.

Saturation maps for a steady flow along the x direction are shown in the right column of Fig. 14, for the three discretizations. They are in excellent agreement.
<figure><img src="image_15.png" /><figcaption>FIG. 13. The dimensionless seepage velocities $\bar{v}'_w$ and $\bar{v}'_n$ of the wetting and nonwetting phases as functions of the dimensionless time $t'$, for the simulations in Fig. 12. 1, 2, and 3 correspond to the left, middle, and right columns of Fig. 12.</figcaption></figure>

## B. Steady-state relative permeabilities as functions of mean saturation

The previous example has shown that identical saturation fields and phase flow rates are eventually reached when a fractured medium is submitted to a macroscopic pressure gradient $\nabla P$, starting from very different initial phase distributions. It is, therefore, possible to define steady state macroscopic phase relative permeabilities for this medium at a given mean saturation $\bar{S}_w$. These relative permeabilities $\bar{K}_{r,i}$ are intrinsic in the sense that they do not depend on the initial conditions; it will be shown on a few examples in Sec. V C that they do not depend either on the magnitude of the applied macroscopic pressure gradient, nor on the viscosity ratio, at least in a reasonable range.

Thus, for the steady flow of two given fluids in a fractured medium, it is possible to relate the mean global phase flow rates $\bar{v}_i$ to the phase pressure gradients by a generalized Darcy's law of form (1)

$$ \bar{v}_i = - \frac{\bar{K} \bar{K}_{r,i}}{\mu_i} (\nabla \bar{P}_i - \rho_i g \bar{e}_z) \quad (i=w,n), \quad (54) $$

where $\bar{K}$ is the macroscopic absolute permeability of the medium. Recall that $\bar{K}$ was investigated by Bogdanov *et al.* [18], for the same type of fractured porous media as considered here.

Equation (54) is the first step towards the upscaling of the steady state two-phase flow problem. It corresponds to the transport equation (1) in an homogeneous material. The conservation equation (2) applies also to $\bar{S}_i$ and $\bar{v}_i$, for steady flows. Finally, the data for the relative permeabilities $\bar{K}_{r,i}$

<figure><img src="image_16.png" /><figcaption>FIG. 14. Illustration of the grid resolution effects. Sections through the mesh of a fractured sampled with $L/R=3$, $N_{fr}=10$, and $\delta_M/R=1/3$, $1/5$, and $1/8$ (top to bottom, left column), and corresponding wetting phase saturation maps for a steady flow oriented from the left to the right of the figure (right column).</figcaption></figure>

from the present numerical calculations are the macroscale counterpart of the constitutive equation (14). The last missing element is an upscaled capillary pressure-saturation relationship. This topic will be discussed at the end of this subsection.

<table><caption>TABLE I. Steady-state permeabilities computed in the same sample along three directions for three meshes with different resolutions (see text in Sec. V A).</caption><thead><tr><th>Direction</th><th>δ<sub>M</sub>/R=3</th><th>δ<sub>M</sub>/R=5</th><th>δ<sub>M</sub>/R=8</th></tr></thead><tbody><tr><td>x</td><td>1.611</td><td>1.579</td><td>1.558</td></tr><tr><td>y</td><td>1.455</td><td>1.431</td><td>1.411</td></tr><tr><td>z</td><td>1.392</td><td>1.369</td><td>1.354</td></tr><tr><td>x</td><td>0.003136</td><td>0.003125</td><td>0.003218</td></tr><tr><td>y</td><td>0.003121</td><td>0.003129</td><td>0.003126</td></tr><tr><td>z</td><td>0.003084</td><td>0.003097</td><td>0.003099</td></tr><tr><td>x</td><td>0.8928</td><td>0.8530</td><td>0.8290</td></tr><tr><td>y</td><td>0.7491</td><td>0.7191</td><td>0.6950</td></tr><tr><td>z</td><td>0.6898</td><td>0.6617</td><td>0.6449</td></tr></tbody></table>
<figure><img src="image_17.png" /></figure>

FIG. 15. Macroscopic relative permeabilities $\bar{K}_{r,i}$ as functions of the mean saturation $\bar{S}_w$. Data are for samples containing 16 (a) or 32 (b) hexagonal fractures. The cell size is $L=4R$. The fractures have a permeability $\sigma'=1$, and $\kappa=10^{-3/2}$. The fluids have equal densities. The symbols are the averages of $\bar{K}_{r,n}$ (□) and $\bar{K}_{r,w}$ (○) over 27 calculations, conducted in $N_r=9$ random realizations, with a pressure gradient with magnitude $P_{0,m}/R$ set along the x, y, and z axes. The horizontal lines show the full variation range of the individual data. The solid lines are the relative permeabilities for the fractures and for the rock matrix, with $n_m=n_f=q=2$. The broken line in (b) is prediction (51) for $\bar{K}_{r,n}$ for infinite plane fractures with the same characteristics and the same global intrinsic permeability.

In view of the large number of parameters, which include the rock and fractures geometry and transport properties, the fluid characteristics, and the coefficients in the constitutive equations, in addition to the initial and boundary conditions, we chose to study only the influence of the mean saturation $\bar{S}_w$ on $\bar{K}_{r,i}$, for two fracture densities and a single typical value of the other parameters. Specifically, the exponents $n_m$, $n_f$, and $q$ are all equal to 2; the viscosity ratio is $\mu_n/\mu_w=10$ and the fluids have the same density; the fracture permeability is $\sigma'=1$ with $\kappa=10^{-3/2}$.

However, we considered samples of size $L=4R$ containing either 16 or 32 randomly located fractures. An example with $N_{fr}=32$ is shown in Fig. 1(c). $N_r=9$ realizations were generated in each case, and the flow equations were solved with a pressure gradient of magnitude $P_{0,m}/R$ set along the x, y, and z axes, successively. In the first case, the percolation probability of the fracture network in a prescribed direction is about 20%, whereas in the latter, it is about 80%. The examples displayed in Figs. 1(b) and 1(c) belong to the first and second families, respectively.

The computations were run starting from initial rest state equilibrium conditions, i.e., with a uniform capillary pressure corresponding to various mean saturations $\bar{S}_w$ in the range 0.1-0.9, until convergence of the saturation field. In addition, a single-phase calculation was performed in order to determine the sample absolute permeability. The relative permeabilities $\bar{K}_{r,i}$ for each case were then deduced from the phase flow rates via Eq. (54).

The results are shown in Figs. 15(a) and 15(b), for the 16- and 32-fracture samples, respectively. The symbols correspond to the statistical averages over 27 calculations, and the

error bars to the full range of variation of the individual data. The solid lines are the relative permeabilities for the fractures and for the rock matrix, which are identical functions of the saturation, in the present case. The mean results are also summarized in Table II.

In spite of the difference in percolation probability between the two cases, the general aspects of the results are similar. The presence of the fractures increases the relative permeability for the nonwetting phase and decreases the relative permeability for the wetting phase, with respect to the intact matrix material. However, the amplitude of these variations is larger for the denser fracture networks.

The strongest effects are observed for the largest saturations and for the nonwetting fluid permeability $\bar{K}_{r,n}$. This is a consequence of the different capillary functions of the fractures and rock matrix [see Eq. (13)]. For the same value of $P_c$, the nonwetting phase saturation is much larger in the fractures than in the surrounding matrix, and the relative permeability $\sigma_{r,n}$ is larger than $K_{r,n}$. Thus, the fractures are preferential paths for the nonwetting phase.

Conversely, $\bar{K}_{r,w}$ is smaller than $K_{r,w}$ in the rock matrix, but this is mostly a consequence of the increase of the absolute permeability induced by the presence of the fractures, $\bar{K}>K_m$. The products $\bar{K}\bar{K}_{r,w}$ and $K_mK_{r,w}$ are identical, which means that the fractures do not significantly affect the wetting phase flow rate, with respect to the intact rock. This is confirmed in Fig. 16, where the normalized phase flow rates $F_i=\bar{K}\bar{K}_{r,i}/K_m=\bar{K}'\bar{K}_{r,i}$ are plotted as functions of $\bar{S}_w$.

Let us finally consider the macroscopic capillary pressure-saturation relationship. The mean saturation $\bar{S}_w$ corresponds
<table><caption>TABLE II. Numerical results in the base case $\sigma'=1$, $G=1$, $\mu_n/\mu_w=10$, $n_f=n_m=q=2$, $\kappa=10^{-3/2}$. The data are averages over the x, y, and z directions for nine random realizations of fractured media containing $N_{fr}=16$ of 32 hexagonal fractures. $P'_{c,r}$ is the initial (rest state) capillary pressure, and $\overline{S}_w$ is the corresponding mean saturation. $\overline{P}'_c$ is the volume-averaged capillary pressure in the steady regime. $K_{r,n}$ and $K_{r,w}$ are the matrix relative permeabilities corresponding to $\overline{S}_w$. $\overline{K}'\overline{K}_{r,i}$ is the normalized flow rate of phase i in the final steady regime. $\overline{K}_{r,i}$ is the corresponding global relative permeability.</caption><thead><tr><th rowspan="2">P'<sub>c,r</sub></th><th colspan="7">N<sub>fr</sub>=16, σ'=1, G=1, μ<sub>n</sub>/μ<sub>w</sub>=10, n<sub>f</sub>=2, n<sub>m</sub>=2, κ=10<sup>-3/2</sup></th></tr><tr><th>S&#773;<sub>w</sub></th><th>P&#773;<sup>'</sup><sub>c</sub></th><th>K<sub>r,n</sub></th><th>K<sub>r,w</sub></th><th>K&#773;<sup>'</sup><sub>r,n</sub></th><th>K&#773;<sub>r,n</sub></th><th>K&#773;<sup>'</sup><sub>r,w</sub></th><th>K&#773;<sub>r,w</sub></th></tr></thead><tbody><tr><td>0.00</td><td>1.0000</td><td>0.000</td><td>0.000</td><td>1.000</td><td>0.000</td><td>0.000</td><td>1.316</td><td>1.000</td></tr><tr><td>0.75</td><td>0.8000</td><td>0.746</td><td>0.040</td><td>0.143</td><td>0.121</td><td>0.092</td><td>0.152</td><td>0.116</td></tr><tr><td>1.00</td><td>0.7071</td><td>1.007</td><td>0.086</td><td>0.0721</td><td>0.209</td><td>0.158</td><td>0.0741</td><td>0.0564</td></tr><tr><td>1.50</td><td>0.5547</td><td>1.509</td><td>0.198</td><td>0.0210</td><td>0.379</td><td>0.288</td><td>0.0211</td><td>0.0161</td></tr><tr><td>2.00</td><td>0.4472</td><td>2.007</td><td>0.306</td><td>0.00745</td><td>0.521</td><td>0.396</td><td>0.00745</td><td>0.00567</td></tr><tr><td>2.50</td><td>0.3714</td><td>2.505</td><td>0.395</td><td>0.00312</td><td>0.632</td><td>0.480</td><td>0.00311</td><td>0.00237</td></tr><tr><td>3.00</td><td>0.3162</td><td>3.004</td><td>0.468</td><td>0.00148</td><td>0.719</td><td>0.546</td><td>0.00148</td><td>0.00113</td></tr><tr><td>5.00</td><td>0.1961</td><td>5.030</td><td>0.646</td><td>0.000167</td><td>0.927</td><td>0.704</td><td>0.000165</td><td>0.000126</td></tr><tr><td>10.00</td><td>0.0995</td><td>10.000</td><td>0.811</td><td>0.000008</td><td>1.109</td><td>0.843</td><td>0.000008</td><td>0.000006</td></tr></tbody></table><table><thead><tr><th rowspan="2">P'<sub>c,r</sub></th><th colspan="7">N<sub>fr</sub>=32, σ'=1, G=1, μ<sub>n</sub>/μ<sub>w</sub>=10, n<sub>f</sub>=2, n<sub>m</sub>=2, κ=10<sup>-3/2</sup></th></tr><tr><th>S&#773;<sub>w</sub></th><th>P&#773;<sup>'</sup><sub>c</sub></th><th>K<sub>r,n</sub></th><th>K<sub>r,w</sub></th><th>K&#773;<sup>'</sup><sub>r,n</sub></th><th>K&#773;<sub>r,n</sub></th><th>K&#773;<sup>'</sup><sub>r,w</sub></th><th>K&#773;<sub>r,w</sub></th></tr></thead><tbody><tr><td>0.00</td><td>1.0000</td><td>0.000</td><td>0.000</td><td>1.000</td><td>0.000</td><td>0.000</td><td>1.657</td><td>1.000</td></tr><tr><td>0.50</td><td>0.8945</td><td>0.465</td><td>0.011</td><td>0.289</td><td>0.165</td><td>0.100</td><td>0.338</td><td>0.204</td></tr><tr><td>0.75</td><td>0.8001</td><td>0.747</td><td>0.040</td><td>0.143</td><td>0.299</td><td>0.180</td><td>0.153</td><td>0.0926</td></tr><tr><td>1.00</td><td>0.7071</td><td>1.009</td><td>0.086</td><td>0.0721</td><td>0.428</td><td>0.258</td><td>0.0739</td><td>0.0446</td></tr><tr><td>1.50</td><td>0.5547</td><td>1.512</td><td>0.198</td><td>0.0210</td><td>0.630</td><td>0.380</td><td>0.0211</td><td>0.0128</td></tr><tr><td>2.00</td><td>0.4472</td><td>2.009</td><td>0.306</td><td>0.00745</td><td>0.795</td><td>0.479</td><td>0.00745</td><td>0.00450</td></tr><tr><td>3.00</td><td>0.3162</td><td>3.005</td><td>0.468</td><td>0.00148</td><td>1.017</td><td>0.614</td><td>0.00148</td><td>0.00089</td></tr><tr><td>5.00</td><td>0.1961</td><td>4.997</td><td>0.646</td><td>0.000167</td><td>1.241</td><td>0.748</td><td>0.000168</td><td>0.000102</td></tr><tr><td>10.00</td><td>0.0995</td><td>10.000</td><td>0.811</td><td>0.000008</td><td>1.438</td><td>0.868</td><td>0.000008</td><td>0.000005</td></tr></tbody></table>

in a rest state to a capillary pressure $P_{c,r}$. Since the interstitial volume in the medium is widely dominated by the pore volume in the rock matrix, $P_{c,r}$ is related to $\overline{S}_w$ by the law (11) for the matrix, with $n=n_m$ and $P_0=P_{0,m}$. When a steady flow takes place through the fractured medium, the saturation and capillary pressure fields are not uniform (see Fig. 12). The volume average of the fluctuations of $S_w$ is 0, since $\overline{S}_w$ is conserved. However, due to the nonlinearity of Eq. (11), the volume average $\overline{P}_c$ is not necessarily equal to $P_{c,r}$. It was calculated for the 27 steady states obtained in the two types of fractured media. It appears that the difference between $\overline{P}_c$ and $P_{c,r}$ is negligible (see Table II). Hence, the volume averaged capillary pressure does not differ between rest state and steady flow.

Recall that in the present case the global volume average is equivalent to an average over the matrix only, since the fracture volume is negligible. The previous statement might not be valid if the matrix contained three-dimensional heterogeneousities such as lenses of a more permeable material. In this case, the concentration of the nonwetting fluid in the permeable region would induce a noticeable decrease of its saturation in the matrix.

### C. Influence of the other parameters

The influence of various parameters, namely, the fracture permeability $\sigma'$, the magnitude G of the driving pressure

gradient, the viscosity ratio $\mu_n/\mu_w$, and the exponents $n_m$ and $n_f$, was briefly tested by varying a single parameter at a time with respect to the base case considered in Sec. V B. Calculations were run for a few values of the mean saturation $\overline{S}_w$, on a single sample of fractured medium, containing $N_{fr}=16$ or 32 fractures. The results are given in Tables III–VI.

#### 1. Fracture permeability

It should be recalled first that the fracture permeability $\sigma$ depends both on the aperture and on the filling permeability [see Eq. (5)]. In the following, the aperture was kept constant. Hence, variations of $\sigma$ correspond to variations of $\kappa$, which, in turn, modify the capillary properties of the fracture [see Eqs. (11) and (13)].

The numerical results for $\sigma'=1$ to 1000 are given in Table III. Recall that $\sigma'=1$ corresponds to the base case in Sec. V B.

Note at first that the normalized flux $\overline{K}\overline{K}_{r,w}$ for $\overline{S}_w=1$ is the fractured medium absolute permeability $\overline{K}$. As observed by Bogdanov et al. [18], for large fracture permeabilities, it increases asymptotically as a linear function of $\sigma'$ for the percolating network ($N_{fr}=32$), whereas it tends toward a finite limit for the nonpercolating network ($N_{fr}=16$). In the latter case, the flow rate is controlled when the fractures are very permeable by the gaps that the fluid has to cross through
<figure><img src="image_18.png" /></figure>

a

<figure><img src="image_19.png" /></figure>

b

<caption>FIG. 16. Normalized phase flow rates $F_n = \bar{K}' \bar{K}_{r,n}$ (□) and $F_w = \bar{K}' \bar{K}_{r,w}$ (○) corresponding to the data in Fig. 15.</caption>

the matrix between the fractures.

When $\sigma$ increases, the wetting fluid relative permeability $\bar{K}_{r,w}$ decreases, but this is merely a consequence of the increase of $\bar{K}$; indeed, the product $\bar{K}\bar{K}_{r,w}$ remains nearly constant. It was already noted in Sec. V B that the fractures do not significantly contribute to the wetting phase flow when $\sigma' = 1$. For larger $\sigma'$, $S_w$ in the fractures, and thus $\sigma_{r,w}$, are smaller. It can easily be checked that with the present parameters, $\sigma_{r,w}$ decreases much faster than $\sigma'$ increases. Hence, the fractures contribute even less to the wetting fluid flow

when $\sigma'$ increases. Recall that this applies only when the variations of $\sigma$ result from a change of the fracture filling permeability.

On the other hand, the fractures are more and more saturated with the nonwetting fluid when $\sigma'$ increases. Hence, when the fracture network percolates and dominates the global absolute permeability $\bar{K}$, the total permeability $\bar{K}\bar{K}_{r,n}$ approaches $\bar{K}$, and $\bar{K}_{r,n}$ tends toward unity. When the network does not percolate, the non-wetting fluid flux increases with $\sigma'$, but it tends toward a finite limit, as does $\bar{K}$. The global relative permeability $\bar{K}_{r,n}$ first increases, and then slightly decreases, with a maximum around $\sigma' = 10$.

## 2. Pressure gradient

Denote by $G$ the dimensionless macroscopic pressure gradient, i.e.,

$$ \overline{\nabla P} = G \frac{P_{0,m}}{R} \qquad (55) $$

It was decreased to 0.1 or increased to 10 in a few cases, as shown in Table IV. Of course, the global absolute (single phase) permeability $\bar{K}$ does not depend on $G$, and it is a constant for each of the two parts of the Table, for $N_{fr} = 16$ or 32.

For low saturations ($\bar{S}_w \le 0.2$), the phase relative permeabilities are nearly unchanged by a change of $G$ from 1 to 10. The variations of $\bar{K}_{r,i}$ never exceed 3%.

For large saturations ($\bar{S}_w \ge 0.7$), the relative permeabilities vary more significantly when $G$ changes from 0.1 to 1. Generally, both $\bar{K}_{r,n}$ and $\bar{K}_{r,w}$ slightly increase when $G$ increases by a factor of 10. The variation is maximum (19%) for the nonwetting phase with a large saturation ($\bar{S}_w = 0.8$). However, the nonwetting phase permeability behaves in the opposite way for $\bar{S}_w = 0.9$, with a decrease by 7%.

Several mechanisms may play a role in these variations. The spatial saturation variations increase with $G$, causing local variations of relative permeability and modifying the hydraulic interactions between the fractures. Strong gradients may even cause the intrusion of wetting fluid in the fractures. This will be investigated in detail in a future work.

Still, these observations are quite consistent with criterion (29). With the present parameters, $G_c$ is equal to $\sqrt{(1-S)/(1+S)}/S$. For $S_w = 0.0995$ ($G_c \approx 11$), $G=1$ and $G=10$ yield identical results. For $S_w = 0.196$ ($G_c \approx 6$), the results for $G=1$ and $G=10$ differ at most by a few percents. The same applies to the data for $G=0.1$ and $G=1$, when $S_w = 0.7$ ($G_c \approx 0.8$). However, $G_c \approx 0.5$ and $G_c \approx 0.3$ for $S_w = 0.8$ and $S_w = 0.9$, respectively. Then, $G=1$ exceeds $2G_c$, and significant saturation changes may occur, thereby modifying the apparent relative permeabilities with respect to the smaller gradient $G=0.1$, by 10–20%.

It is gratifying to note that Eq. (29) indeed provides a reliable a priori criterion stating whether intrinsic macroscopic relative permeabilities can or cannot be defined. It will be seen in Sec. VI A that the same criterion also deter-
<table><caption>TABLE III. Numerical results for individual realizations with a nonpercolating ($N_{fr}=16$) or a percolating ($N_{fr}=32$) fracture network, when the fracture permeability $\sigma'$ varies with respect to the base case in Table II. Same notations as in Table II.</caption><caption>G = 1, μ<sub>n</sub>/μ<sub>w</sub> = 10, n<sub>f</sub> = 2, n<sub>m</sub> = 2, κ = 10<sup>-3/2</sup></caption><thead><tr><th rowspan="2">σ'</th><th rowspan="2">P'<sub>c,r</sub></th><th rowspan="2">S&#773;<sub>w</sub></th><th colspan="4">N<sub>fr</sub> = 16</th><th colspan="4">N<sub>fr</sub> = 32</th></tr><tr><th>K&#773;<sup>'</sup> K&#773;<sub>r,n</sub></th><th>K&#773;<sub>r,n</sub></th><th>K&#773;<sup>'</sup> K&#773;<sub>r,w</sub></th><th>K&#773;<sub>r,w</sub></th><th>K&#773;<sup>'</sup> K&#773;<sub>r,n</sub></th><th>K&#773;<sub>r,n</sub></th><th>K&#773;<sup>'</sup> K&#773;<sub>r,w</sub></th><th>K&#773;<sub>r,w</sub></th></tr></thead><tbody><tr><td>1</td><td>0.00</td><td>1.000</td><td>0.000</td><td>0.000</td><td>1.318</td><td>1.000</td><td>0.000</td><td>0.000</td><td>1.66</td><td>1.000</td></tr><tr><td>10</td><td>0.00</td><td>1.000</td><td>0.000</td><td>0.000</td><td>2.239</td><td>1.000</td><td>0.000</td><td>0.000</td><td>4.45</td><td>1.000</td></tr><tr><td>100</td><td>0.00</td><td>1.000</td><td>0.000</td><td>0.000</td><td>3.521</td><td>1.000</td><td>0.000</td><td>0.000</td><td>19.6</td><td>1.000</td></tr><tr><td>1000</td><td>0.00</td><td>1.000</td><td>0.000</td><td>0.000</td><td>4.298</td><td>1.000</td><td>0.000</td><td>0.000</td><td>154.</td><td>1.000</td></tr><tr><td>1</td><td>1.50</td><td>0.555</td><td>0.383</td><td>0.290</td><td>0.0212</td><td>0.0161</td><td>0.630</td><td>0.380</td><td>0.0211</td><td>0.0128</td></tr><tr><td>10</td><td>1.50</td><td>0.555</td><td>0.686</td><td>0.303</td><td>0.0213</td><td>0.0096</td><td>2.36</td><td>0.532</td><td>0.0215</td><td>0.00483</td></tr><tr><td>100</td><td>1.50</td><td>0.555</td><td>0.968</td><td>0.270</td><td>0.0211</td><td>0.0065</td><td>16.0</td><td>0.806</td><td>0.0223</td><td>0.00116</td></tr><tr><td>1000</td><td>1.50</td><td>0.555</td><td>1.078</td><td>0.251</td><td>0.0211</td><td>0.0057</td><td>149.</td><td>0.967</td><td>0.0225</td><td>0.00016</td></tr><tr><td>1</td><td>3.00</td><td>0.316</td><td>0.721</td><td>0.547</td><td>0.00148</td><td>0.00113</td><td>1.02</td><td>0.613</td><td>0.00148</td><td>0.000895</td></tr><tr><td>10</td><td>3.00</td><td>0.316</td><td>1.266</td><td>0.563</td><td>0.00147</td><td>0.00066</td><td>3.17</td><td>0.713</td><td>0.00148</td><td>0.000333</td></tr><tr><td>100</td><td>3.00</td><td>0.316</td><td>1.910</td><td>0.533</td><td>0.00143</td><td>0.00044</td><td>17.3</td><td>0.876</td><td>0.00149</td><td>0.000077</td></tr><tr><td>1000</td><td>3.00</td><td>0.316</td><td>2.229</td><td>0.513</td><td>0.00142</td><td>0.00038</td><td>151.</td><td>0.978</td><td>0.00149</td><td>0.000010</td></tr></tbody></table>

mines whether the macroscopic two-phase flow properties of the fractured porous medium can be accurately evaluated by using a simple first-order model.

### 3. Viscosity ratio

A few calculations were run with a viscosity contrast $\mu_n/\mu_w=1$ instead of 10 in the base case of Sec. V B. The

results in Table V show a total absence of influence of this parameter. This also supports the macroscopic description (54).

### 4. Capillary functions

Finally, the exponents $n_f$ and $n_m$ were changed from 2 to 3 in a few cases. The results are presented in Table VI. The

<table><caption>TABLE IV. Numerical results for individual realizations with a nonpercolating ($N_{fr}=16$) or a percolating ($N_{fr}=32$) fracture network, when the macroscopic pressure gradient is changed by a factor of 10 with respect to the base case in Table II. Same notations as in Table II. Lines starting with a star correspond to the predictions of model (57) and (59).</caption><caption>σ' = 1, μ<sub>n</sub>/μ<sub>w</sub> = 10, n<sub>f</sub> = 2, n<sub>m</sub> = 2, κ = 10<sup>-3/2</sup></caption><thead><tr><th rowspan="2">G</th><th rowspan="2">P'<sub>c,r</sub></th><th rowspan="2">S&#773;<sub>w</sub></th><th colspan="4">N<sub>fr</sub> = 16</th><th colspan="4">N<sub>fr</sub> = 32</th></tr><tr><th>K&#773;<sup>'</sup> K&#773;<sub>r,n</sub></th><th>K&#773;<sub>r,n</sub></th><th>K&#773;<sup>'</sup> K&#773;<sub>r,w</sub></th><th>K&#773;<sub>r,w</sub></th><th>K&#773;<sup>'</sup> K&#773;<sub>r,n</sub></th><th>K&#773;<sub>r,n</sub></th><th>K&#773;<sup>'</sup> K&#773;<sub>r,w</sub></th><th>K&#773;<sub>r,w</sub></th></tr></thead><tbody><tr><td>1.0</td><td>0.50</td><td>0.895</td><td></td><td></td><td></td><td></td><td>0.172</td><td>0.104</td><td>0.332</td><td>0.200</td></tr><tr><td>0.1</td><td>0.50</td><td>0.895</td><td></td><td></td><td></td><td></td><td>0.183</td><td>0.111</td><td>0.289</td><td>0.175</td></tr><tr><td>*</td><td>0.50</td><td>0.895</td><td></td><td></td><td></td><td></td><td></td><td>0.110</td><td></td><td>0.174</td></tr><tr><td>1.0</td><td>0.75</td><td>0.800</td><td>0.131</td><td>0.0992</td><td>0.152</td><td>0.115</td><td>0.305</td><td>0.1838</td><td>0.151</td><td>0.0911</td></tr><tr><td>0.1</td><td>0.75</td><td>0.800</td><td>0.110</td><td>0.0829</td><td>0.143</td><td>0.109</td><td>0.283</td><td>0.1710</td><td>0.143</td><td>0.0864</td></tr><tr><td>*</td><td>0.75</td><td>0.800</td><td></td><td>0.0817</td><td></td><td>0.109</td><td></td><td></td><td>0.1716</td><td>0.0863</td></tr><tr><td>1.0</td><td>1.00</td><td>0.707</td><td>0.215</td><td>0.1628</td><td>0.0743</td><td>0.0564</td><td>0.420</td><td>0.253</td><td>0.0739</td><td>0.0446</td></tr><tr><td>0.1</td><td>1.00</td><td>0.707</td><td>0.197</td><td>0.1493</td><td>0.0721</td><td>0.0548</td><td>0.400</td><td>0.242</td><td>0.0721</td><td>0.0435</td></tr><tr><td>*</td><td>1.00</td><td>0.707</td><td></td><td>0.1485</td><td></td><td>0.0548</td><td></td><td>0.241</td><td></td><td>0.0435</td></tr><tr><td>1.0</td><td>5.00</td><td>0.196</td><td>0.927</td><td>0.7032</td><td>0.000168</td><td>0.000128</td><td>1.240</td><td>0.748</td><td>0.000171</td><td>0.000103</td></tr><tr><td>10.</td><td>5.00</td><td>0.196</td><td>0.933</td><td>0.7081</td><td>0.000174</td><td>0.000132</td><td>1.252</td><td>0.755</td><td>0.000175</td><td>0.000105</td></tr><tr><td>*</td><td>5.00</td><td>0.196</td><td></td><td>0.7078</td><td></td><td>0.000127</td><td></td><td>0.748</td><td></td><td>0.000101</td></tr><tr><td>1.0</td><td>10.00</td><td>0.0995</td><td>1.111</td><td>0.8430</td><td>0.000008</td><td>0.000006</td><td>1.439</td><td>0.868</td><td>0.000008</td><td>0.000005</td></tr><tr><td>10.</td><td>10.00</td><td>0.0995</td><td>1.111</td><td>0.8430</td><td>0.000008</td><td>0.000006</td><td>1.441</td><td>0.869</td><td>0.000008</td><td>0.000005</td></tr><tr><td>*</td><td>10.00</td><td>0.0995</td><td></td><td>0.8464</td><td></td><td>0.000006</td><td></td><td>0.868</td><td></td><td>0.000005</td></tr></tbody></table>
<table><caption>TABLE V. Numerical results for individual realizations with a nonpercolating ($N_{fr}$=16) or a percolating ($N_{fr}$=32) fracture network, when the viscosity contrast $\mu_n/\mu_w$ is changed by a factor of 10 with respect to the base case in Table II. Same notations as in Table II.</caption><caption>σ' = 1, G = 1, n<sub>f</sub> = 2, n<sub>m</sub> = 2, κ = 10<sup>-3/2</sup></caption><thead><tr><th rowspan="2">μ<sub>n</sub>/μ<sub>w</sub></th><th rowspan="2">P'<sub>c,r</sub></th><th rowspan="2">S&#773;<sub>w</sub></th><th colspan="3">N<sub>fr</sub> = 16</th><th colspan="3">N<sub>fr</sub> = 32</th></tr><tr><th>K&#773; K&#773;<sub>r,n</sub></th><th>K<sub>r,n</sub></th><th>K&#773; K&#773;<sub>r,w</sub></th><th>K<sub>r,w</sub></th><th>K&#773; K&#773;<sub>r,n</sub></th><th>K<sub>r,n</sub></th><th>K&#773; K&#773;<sub>r,w</sub></th><th>K<sub>r,w</sub></th></tr></thead><tbody><tr><td>10</td><td>1.00</td><td>0.707</td><td>0.215</td><td>0.163</td><td>0.0743</td><td>0.0564</td><td>0.420</td><td>0.253</td><td>0.0739</td><td>0.0446</td></tr><tr><td>1</td><td>1.00</td><td>0.707</td><td>0.215</td><td>0.163</td><td>0.0743</td><td>0.0565</td><td>0.420</td><td>0.253</td><td>0.0739</td><td>0.0446</td></tr><tr><td>10</td><td>3.00</td><td>0.316</td><td>0.721</td><td>0.547</td><td>0.00148</td><td>0.00113</td><td></td><td></td><td></td><td></td></tr><tr><td>1</td><td>3.00</td><td>0.316</td><td>0.721</td><td>0.547</td><td>0.00148</td><td>0.00112</td><td></td><td></td><td></td><td></td></tr></tbody></table>

change of $n_f$ induces negligible changes in the relative permeabilities. This is because the fractures are nearly saturated with nonwetting fluid both for $n_f=2$ and $n_f=3$. A similar behavior can be expected for any fracture capillary function that yields $S_w \ll 1$.

A change in the capillary properties of the rock matrix has more dramatic effects. First, for a given initial capillary pressure $P_{c,r}$, it changes the mean saturation $S_w$, as shown by Table VI. However, for the same saturation, it also modifies the relative permeabilities in the matrix [see Eq. (14)]. Therefore, $K_{r,w}$ for $n_m=3$ does not correspond to the value obtained with the same mean saturation when $n_m=2$. For instance, with $N_{fr}=16$, $n_m=2$ and $S_w=0.63$, an interpolation of the data in Table II yields $K_{r,w} \approx 0.038$, instead of 0.086 for $n_m=3$ in Table VI. The nonwetting phase global relative permeabilities are in better agreement, but this is because the model for $K_{r,n}$ does not involve the exponent $n_m$ [see Eq. (16)].

## VI. DISCUSSION

### A. Comparison with a capillarity dominated model

Underlying the macroscopic description of steady two-phase flows in terms of relative permeabilities independent of the driving pressure gradient magnitude is the assumption that the spatial phase distribution is not significantly influenced by the flow, with respect to the equilibrium rest state. This corresponds to the small value of the capillary number (12) and to criterion (29).

In this approximation, the local distribution of the phases in the pore volume is assumed to be determined by capillary forces only, and each phase flows through its own system of channels as if the other phase were immobile. Thus, the two flows are actually decoupled, and can be treated as two single-phase flows in a medium where the local permeability is determined by saturation, via the local relative permeabilities.

These considerations, together with the additional feature that the fractures are nearly saturated with the nonwetting fluid, can be applied to devise a simple model for the prediction of the global relative permeabilities.

Recall first that the macroscopic single-phase permeability of a fractured medium can be written as (see Bogdanov et al. [18])

$$ \bar{K} = K_m \chi(\rho', \sigma', \omega'), \quad (56) $$

where $\rho'$ is a measure of the network density, and $\omega'$ is the hydraulic resistance of the fractures to cross flow, which is taken equal to 0 in the present case. In the following, we write in short $\chi(\sigma')$, since the two other parameters are fixed. The function $\chi$ was tabulated by Bogdanov et al. [18] for various values of the density and fracture permeability; it can also easily be computed for any given fractured medium, since it involves only single-phase flow.

Suppose that in the rest state the mean saturation $S_w$ corresponds to a capillary pressure $P_c^0$. The saturation in the matrix $S_{w,m}$ is nearly equal to $S_w$, and the associated relative

<table><caption>TABLE VI. Numerical results for individual realizations with a nonpercolating $N_{fr}=16$) or a percolating ($N_{fr}=32$) fracture network, when the exponents $n_f$ or $n_m$ are modified with respect to the base case in Table II. Same notations as in Table II.</caption><caption>σ' = 1, G = 1, μ<sub>n</sub>/μ<sub>w</sub> = 10, κ = 10<sup>-3/2</sup></caption><thead><tr><th rowspan="2">n<sub>f</sub></th><th rowspan="2">n<sub>m</sub></th><th rowspan="2">P'<sub>c,r</sub></th><th rowspan="2">S&#773;<sub>w</sub></th><th colspan="3">N<sub>fr</sub> = 16</th><th colspan="3">N<sub>fr</sub> = 32</th></tr><tr><th>K&#773; K&#773;<sub>r,n</sub></th><th>K<sub>r,n</sub></th><th>K&#773; K&#773;<sub>r,w</sub></th><th>K<sub>r,w</sub></th><th>K&#773; K&#773;<sub>r,n</sub></th><th>K<sub>r,n</sub></th><th>K&#773; K&#773;<sub>r,w</sub></th><th>K<sub>r,w</sub></th></tr></thead><tbody><tr><td>2</td><td>2</td><td>1.00</td><td>0.707</td><td>0.215</td><td>0.163</td><td>0.0743</td><td>0.0564</td><td>0.420</td><td>0.253</td><td>0.0739</td><td>0.0446</td></tr><tr><td>3</td><td>2</td><td>1.00</td><td>0.707</td><td>0.219</td><td>0.166</td><td>0.0744</td><td>0.0565</td><td>0.436</td><td>0.263</td><td>0.0740</td><td>0.0446</td></tr><tr><td>2</td><td>3</td><td>1.00</td><td>0.630</td><td>0.308</td><td>0.233</td><td>0.1127</td><td>0.0857</td><td>0.540</td><td>0.325</td><td>0.1116</td><td>0.0673</td></tr><tr><td>2</td><td>2</td><td>3.00</td><td>0.316</td><td>0.721</td><td>0.547</td><td>0.00148</td><td>0.00113</td><td></td><td></td><td></td><td></td></tr><tr><td>3</td><td>2</td><td>3.00</td><td>0.316</td><td>0.724</td><td>0.550</td><td>0.00148</td><td>0.00113</td><td></td><td></td><td></td><td></td></tr><tr><td>2</td><td>3</td><td>3.00</td><td>0.109</td><td>1.091</td><td>0.828</td><td>0.00019</td><td>0.00015</td><td></td><td></td><td></td><td></td></tr></tbody></table>
<table><caption>TABLE VII. Comparison of the full solution of the two-phase flow equations in single realizations of porous media containing a non-percolating ($N_{fr}$=16) or a percolating ($N_{fr}$=32) fracture network, with models (57) and (59). Data are for $G=1$, $\mu_n/\mu_w=10$, $n_f=n_m=q=2$, $\kappa=10^{-3/2}$.</caption><thead><tr><th rowspan="3">S&#773;<sub>w</sub></th><th rowspan="3">σ'</th><th colspan="4">N<sub>fr</sub> = 16</th><th colspan="4">N<sub>fr</sub> = 32</th></tr><tr><th colspan="2">K&#773;<sub>r,n</sub></th><th colspan="2">10<sup>3</sup> K&#773;<sub>r,w</sub></th><th colspan="2">K&#773;<sub>r,n</sub></th><th colspan="2">10<sup>3</sup> K&#773;<sub>r,w</sub></th></tr><tr><th>Num.</th><th>Model</th><th>Num.</th><th>Model</th><th>Num.</th><th>Model</th><th>Num.</th><th>Model</th></tr></thead><tbody><tr><td rowspan="4">0.316</td><td>1</td><td>0.5467</td><td>0.5497</td><td>1.125</td><td>1.126</td><td>0.6129</td><td>0.6102</td><td>0.8945</td><td>0.8935</td></tr><tr><td>10</td><td>0.5626</td><td>0.5515</td><td>0.6620</td><td>0.6693</td><td>0.7126</td><td>0.7140</td><td>0.3325</td><td>0.3330</td></tr><tr><td>100</td><td>0.5333</td><td>0.5049</td><td>0.4417</td><td>0.4572</td><td>0.8759</td><td>0.8721</td><td>0.07670</td><td>0.07641</td></tr><tr><td>1000</td><td>0.5126</td><td>0.4765</td><td>0.3841</td><td>0.4016</td><td>0.9784</td><td>0.9765</td><td>0.01018</td><td>0.01009</td></tr><tr><td rowspan="4">0.555</td><td>1</td><td>0.2905</td><td>0.2834</td><td>16.08</td><td>15.97</td><td>0.3800</td><td>0.3708</td><td>12.72</td><td>12.68</td></tr><tr><td>10</td><td>0.3032</td><td>0.2749</td><td>9.606</td><td>9.495</td><td>0.5315</td><td>0.5244</td><td>4.827</td><td>4.725</td></tr><tr><td>100</td><td>0.2699</td><td>0.2272</td><td>6.490</td><td>6.486</td><td>0.8062</td><td>0.7998</td><td>1.155</td><td>1.084</td></tr><tr><td>1000</td><td>0.2509</td><td>0.2042</td><td>5.674</td><td>5.698</td><td>0.9671</td><td>0.9644</td><td>0.1548</td><td>0.1431</td></tr></tbody></table>

permeabilities are denoted by $K_{r,i}^0$. The saturation in the fractures results from Eqs. (11) and (13), and is generally close to 1, except for very large mean saturation $\bar{S}_w$. The corresponding relative permeabilities in the fractures are $\sigma_{r,i}^0$, with $\sigma_{r,w}^0 \ll 1$.

Consider first the wetting-phase flow. As already noted in Sec. V B, the fractures have a negligible contribution to the flow, since $\sigma_{r,w}^0 \ll 1$. In the present model, the fractures do not present any resistance to cross-flow, and thus, they are neutral with respect to the wetting fluid. Hence, the wetting fluids flows through a uniform medium with apparent permeability $K_m K_{r,w}^0$, and the global wetting-phase relative permeability is

$$ \bar{K}_{r,w} = \frac{K_m K_{r,w}^0}{\bar{K}} = \frac{K_{r,w}^0}{\chi(\sigma')} \quad (57) $$

On the other hand, the nonwetting fluid flows through a fractured porous medium with apparent matrix permeability $K_m K_{r,n}^0$ and fracture permeability $\sigma \sigma_{r,n}^0$. Denote by $\sigma'_n$ the dimensionless ratio defined similarly to $\sigma'$ for the nonwetting phase [see Eq. (48)]

$$ \sigma'_{n} = \frac{\sigma \sigma_{r,n}^{0}}{L K_{m} K_{r,n}^{0}} = \sigma' \frac{\sigma_{r,n}^{0}}{K_{r,n}^{0}} \quad (58) $$

Then, the global apparent permeability of the fractured medium is given by Eq. (56) as $K_m K_{r,n}^0 \chi(\sigma'_n)$, and the corresponding relative permeability for the nonwetting phase is

$$ \bar{K}_{r,n} = \frac{K_m K_{r,n}^0 \chi(\sigma'_n)}{\bar{K}} = K_{r,n}^0 \frac{\chi(\sigma'_n)}{\chi(\sigma')} \quad (59) $$

Predictions (57) and (59) are compared with the full solution of the two-phase flow equations in Table VII, for single realizations of porous media containing a nonpercolating ($N_{fr}$=16) or percolating ($N_{fr}$=32) fracture network, for various fracture permeabilities $\sigma' = 1$ to $10^3$ and two moderate values of $\bar{S}_w$, 0.32 and 0.55. The other parameters are

set as in Sec. V B. The critical gradient $G_c$ is about 4 and 2 for $\bar{S}_w=0.32$ and 0.55, respectively; therefore, criterion (29) is fairly well satisfied. The function $\chi$ was computed beforehand for a few values of $\sigma'$, and the values required for the application of Eqs. (57) and (59) were obtained by a spline interpolation of the data.

Note first that the agreement is always very good for the wetting phase. The only significant differences are an underestimation by about 7% for the networks of very permeable fractures ($\sigma \ge 100$) and for the largest saturation $\bar{S}_w=0.55$. In this situation, the fracture network seems to slightly contribute to the wetting fluid flow.

The agreement is also very good for the nonwetting phase in the sample containing a percolating fracture network. The differences between the numerical solution and the prediction (59) is always less than 1%, except for the least permeable fractures $\sigma' = 1$, where it reaches 3%.

The two previous cases are the most simple, since each phase flows in a single domain, namely, the wetting phase in the matrix and the nonwetting fluid in the percolating fracture network. Nonwetting phase flow when the network does not percolate is more complex, since the fluid has to flow through both the fractures and the matrix.

Nevertheless, a good agreement is also observed for the nonwetting phase for $\bar{S}_w=0.32$. Prediction (59) slightly underestimates $\bar{K}_{r,n}$, but the difference is small. It increases with the fracture permeability, but remains smaller than 8% for $\sigma' = 10^3$.

The largest errors occur for very permeable fractures and $\bar{S}_w=0.55$. The underestimation of $\bar{K}_{r,n}$ by Eq. (59) reaches 15% and 18% for $\sigma' = 10^2$ and $10^3$, respectively.

Another comparison is provided in Table IV, with numerical simulations for small and large values of the mean saturation, and different magnitudes of the macroscopic pressure gradient $G$. In all the cases, predictions (57) and (59) are in excellent agreement with the full numerical solution for the smallest value of $G$, as expected in view of the simplifying assumption of the model.
In summary, the predictions of the model are quite satisfactory, in the range of pressure gradient defined by criterion (29).

## B. Comparison with the infinite plane model

It was shown by Koudina et al. [32] that in the limit of large densities, the single phase permeability of fracture networks is reasonably well approximated by a simple model of infinite plane fractures, where the volumetric fracture area is conserved. The same applies for fractured porous media when the fractures are very permeable, since the contribution of the matrix to the flow is negligible. Therefore, one may wonder whether a similar property applies for the nonwetting phase in the case of two-phase flow through fractured porous media, since the fracture transmissivity is generally much larger than the matrix transmissivity, even for moderate $\sigma'$, due to the larger relative permeability in the fractures [see Eq. (58)].

This was tested by comparing in Fig. 15(b) the results of the full solution of the two-phase flow equations with the analytical expression (51), for identical fracture aperture $b$ and permeability ratio $\kappa$. The volumetric area of the fractures was set in Eq. (51) in order to match the global intrinsic permeability $\bar{K}$ of the random fractured media. It appears that Eq. (51) does not agree with the numerical data. For $\bar{S}_w=0.8$, Eq. (51) overestimates $\bar{K}_{r,n}$ by a factor larger than two, even though $\sigma'_n$ is equal to about 20. For $\bar{S}_w=0.9$, with $\sigma'_n \approx 85$, the overestimation is a factor larger than 3.

Paradoxically, a better agreement is observed for small saturations $\bar{S}_w$, where the model is less justified, since $\sigma'_n$ is not large compared to 1 (e.g., $\sigma'_n \approx 2\sigma'=2$ for $\bar{S}_w=0.3$). This is probably because the fractures do not introduce strong heterogeneities, and any reasonable mixing rule would yield reasonable estimates of the transport coefficient.

Hence, it can be concluded that the infinite plane fracture model fails to account for the two-phase flow properties of random fractured porous media, even when the fracture density is fitted to match the single-phase permeability.

## VII. CONCLUDING REMARKS

We presented in this paper a numerical tool for the simulation of two-phase flows in fractured porous media, together with a set of applications that demonstrated its ability to handle steady or transient flows in complex random media.

Owing to the large number of physical parameters, the scope of the simulations was necessarily restricted to a few representative situations. Still, a systematic study of the macroscopic flow properties was conducted as a function of the global mean saturation, for two types of media containing percolating or nonpercolating fracture networks, and the influence of many other parameters was briefly addressed.

Among the main results is the demonstration that the large scale steady-state flow properties can be estimated from the rest state phase distribution, if criterion (29) is fullfilled. However, this still requires the solution of the flow equations in the actual fractured medium complex geometry, since neither the absolute nor the relative permeabilities are accurately predicted by an infinite plane model for small or moderate fracture density. For faster flows, the phase distribution is influenced by the fluid motion, and a full solution of the coupled two-phase flow equations cannot be avoided. This is the macroscopic counterpart of the property that the local capillary function measured in a rest state applies under flow if, but only if, the capillary number is small enough.

Although the results correspond to certain types of constitutive equations (capillary function and relative permeabilities), many of the main features are expected to apply in a more general setting. For instance, criterion (29) may have to be applied either to the matrix or to the fractures or in a mixed form, but it probably always applies in some sense. Similarly, time constant (34b) may actually have to be evaluated for the wetting or nonwetting phase, in the fractures or in the matrix, depending on which flow is physically expected to be limitant, but probably always have the same general form.

The numerical model should be extended by including a few additional features. Some are minor changes, such as the introduction of residual saturations. Hysteretic local capillary functions can also, in principle, be easily introduced, except for possible numerical convergence difficulties, and would probably induce new types of behaviors on the macroscopic scale. Mixed-wettability cannot be accounted for in the present implementation of the model, but could probably be introduced, provided that consistent constitutive equations are used (see, e.g., Kjosavik et al. [47]). Finally, fractures can act as obstacles to the flow in their normal directions, for several reasons such as capillary effects, filling with a material less permeable than the matrix rock, or alteration of the matrix material along the fracture caused, for instance, by solute precipitation. This would result in apparent pressure discontinuities on the intermediate scale of the present description. Technically, accounting for these effects requires to decouple the pressures on either side and in the middle of the fractures, as done by Bogdanov et al. [18] for single-phase flow. However, capillary barriers would require additional constitutive equations, which are apparently not yet clearly stated.

## ACKNOWLEDGMENTS

Most computations were performed at CINES (subsidized by the MENESR) whose support is gratefully acknowledged.
[1] Y.C. Yortsos, in *Dynamics of Fluids in Fractured Rock*, Geophysical Monograph Series, Vol. 122, edited by B. Faybishenko, P. A. Witherspoon, and S. M. Benson (AGU, Washington, D.C., 2000).

[2] C.M. Marle, *Multiphase Flows in Porous Media* (Technip, Paris, 1981).

[3] J. Bear and Y. Bachmat, *Introduction to Modeling of Transport Phenomena in Porous Media* (Kluwer, Dordrecht, 1991).

[4] F.A.L. Dullien, *Porous Media. Fluid Transport and Pore Structure* (Academic Press, San Diego, 1992).

[5] J.S.Y. Wang and T.N. Narasimhan, Water Resour. Res. **21**, 1861 (1985).

[6] M. Muskat and M.W. Meres, *Physics (N.Y.)* **7**, 346 (1936).

[7] M. Muskat, R.D. Wyckoff, H.G. Botset and M.W. Meres, *Trans. AIME* **123**, 69 (1937).

[8] A. Bourgeat, in *Homogenization and Porous Media*, edited by U. Hornung (Springer, New York, 1997).

[9] G.I. Barenblatt and Yu.P. Zheltov, Dokl. Akad. Nauk SSSR **132**, 545 (1960) [Sov. Phys. Dokl. **5**, 522 (1960)].

[10] G.I. Barenblatt, Yu.P. Zheltov, and I.N. Kochina, J. Appl. Math. Mech. **24**, 852 (1960).

[11] J.R. Warren and P.J. Root, SPEJ **3**, 245 (1963).

[12] H. Kazemi and J. R. Gilman, in *Flow and Contaminant Transport in Fractured Rock*, edited by J. Bear, C.F. Tsang, and G. de Marsily (Academic Press, San Diego, 1993).

[13] M. Panfilov, *Macroscale Models of Flow Through Highly Heterogeneous Porous Media* (Kluwer, Dordrecht, 2000).

[14] D.H. Fenwick and M.J. Blunt, Adv. Water Resour. **21**, 121 (1998).

[15] T.W. Patzek, SPEJ **6**, 144 (2001).

[16] R.G. Hughes and M.J. Blunt, Adv. Water Resour. **24**, 409 (2001).

[17] M.J. Blunt, M.D. Jackson, M. Piri, and P.H. Valvatne, Adv. Water Resour. **25**, 1069 (2002).

[18] I.I. Bogdanov, V.V. Mourzenko, J.-F. Thovert and P.M. Adler, *Water Resour. Res.* **39**, 10.1029/2001WR000756 (2003).

[19] R. Kaiser, O. Kolditz, and T. Rother, in *Proceedings of the 28th IAHR Conference on Hydraulic Engineering for Sustainable Water Resources Management, Graz, Austria, 1999*, edited by H. Bergmann, R. Krainer, and H. Breinhaelter (Techn. University Graz, Graz, 1999).

[20] Y.S. Wu and K. Pruess, *Adv. Water Resour.* **23**, 699 (2000).

[21] A.J.A. Unger, P.A. Forsyth, and E.A. Sudicky, *Adv. Water Resour.* **19**, 1 (1996).

[22] K.J. Slough, E.A. Sudicky, and P.A. Forsyth, *J. Contam. Hydrol.* **40**, 107 (1999).

[23] P. Bastian, Z. Chen, R. Ewing, R. Helmig, H. Jakobs, and V. Reichenberger, in *Numerical Treatment of Multiphase Flows in Porous Media*, edited by Z. Chen, R.E. Ewing, and Z.C. Shi (Springer, Heidelberg, 2000).

[24] J. Schöberl, *Comput. Vis. Sci.* **1**, 41 (1997).

[25] J. Schmittbuhl, A. Hansen, H. Auradou, and K.J. Måløy, *Phys. Rev. E* **61**, 3985 (2000).

[26] H. Auradou, K.J. Måløy, J. Schmittbuhl, and A. Hansen, *Transp. Porous Media* **50**, 267 (2003).

[27] E. Aker, K.J. Måløy, and A. Hansen, *Phys. Rev. E* **61**, 2936 (2000).

[28] E.G. Flekkøy, J. Schmittbuhl, F. Løvholt, U. Oxaal, K.J. Måløy, and P. Aagaard, *Phys. Rev. E* **65**, 036312 (2002).

[29] H.A. Knudsen and A. Hansen, *Phys. Rev. E* **65**, 056310 (2002).

[30] Y. Méheust, G. Lóvoll, K.J. Måløy, and J. Schmittbuhl, *Phys. Rev. E* **66**, 051603 (2002).

[31] S.R. Brown, R.L. Kranz, and B.P. Bonner, *Geophys. Res. Lett.* **13**, 1430 (1986).

[32] N. Koudina, R. Gonzalez Garcia, J.-F. Thovert, and P.M. Adler, *Phys. Rev. E* **57**, 4466 (1998).

[33] I.I. Bogdanov, V.V. Mourzenko, J.-F. Thovert, and P.M. Adler, *Water Resour. Res.* **39**, 10.1029/2000WR000080 (2003).

[34] P.M. Adler and J.-F. Thovert, *Fractures and Fracture Networks* (Kluwer, Dordrecht, 1999).

[35] D.T. Snow, *Water Resour. Res.* **5**, 1273 (1969).

[36] L.A. Richards, *Physics (N.Y.)* **1**, 318 (1931).

[37] P.M. Adler and H. Brenner, *Annu. Rev. Fluid Mech.* **20**, 35 (1988).

[38] S. Reitsma and B.H. Kueper, *Water Resour. Res.* **30**, 865 (1994).

[39] J. Chen, J.W. Hopmans, and M.E. Grismer, *Adv. Water Resour.* **22**, 479 (1999).

[40] M.T. van Genuchten, *Soil Sci. Soc. Am. J.* **44**, 892 (1980).

[41] P. Persoff and K. Pruess, *Water Resour. Res.* **31**, 1175 (1995).

[42] M. Fourar and R. Lenormand, *Oil Gas Sci. Technol. Rev. IFP* **55**, 259 (2002).

[43] E.S. Romm, *Fluid Flow in Fractured Rocks* (Nedra, Moscow, 1966).

[44] M. Fourar and S. Bories, *Int. J. Multiphase Flow* **21**, 621 (1995).

[45] M.A. Celia, E.T. Bouloutas, and R.L. Zarba, *Water Resour. Res.* **26**, 1483 (1990).

[46] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, *Numerical Recipes in Fortran: The Art of Scientific Computing*, 2nd ed. (Cambridge University Press, Cambridge, 1994).

[47] A. Kjosavik, J.K. Ringen, and S.M. Skjaeland, *SPE paper 59314*, presented at the 2000 SPE/DOE Improved Oil Recovery Symposium, Tulsa, Oklahoma, 2000.