Control-Volume Model for Simulation of

Water Injection in Fractured Media:

Incorporating Matrix Heterogeneity and
Reservoir Wettability Effects
J.E.P. Monteagudo* and A. Firoozabadi, SPE, Reservoir Engineering Research Institute

Summary increased (Durlofsky 1994). In the standard CV method, Delaunay

The control-volume discrete-fracture (CVDF) model is extended triangles are locally homogeneous and the polygonal CV cell may
to incorporate heterogeneity in rock and in rock-fluid properties. A be heterogeneous (see Fig. 1a). For accurate streamlines, several
novel algorithm is proposed to model strong water-wetting with authors (Verma 1996; Edwards 2002; Edwards and Rogers 1998;
zero capillary pressure in the fractures. The extended method is Prevost 2000; Aavatsmark et al. 1998a) have proposed that the
used to simulate: (1) oil production in a layered faulted reservoir, polygonal CV cell must be locally homogeneous, implying het-
(2) laboratory displacement tests in a stack of matrix blocks with erogeneous Delaunay triangles (see Fig. 1b). The latter configura-
a large contrast in fracture and matrix capillary pressure functions, tion, however, generates additional problems in the simulation of
and (3) water injection in 2D and 3D fractured media with mixed- multiphase flow in porous media. Basically, from mesh generation
wettability state. Our results show that the algorithm is suitable for standpoint, it may not be possible to generate an unstructured mesh
the simulation of water injection in heterogeneous porous media where the boundaries of the CV median-dual cell conform to het-
both in water-wet and mixed-wettability states. The novel ap- erogeneous interfaces in the domain. Conforming mesh is impor-
proach with zero fracture capillary and nonzero matrix capillary tant for the discrete-fracture approach. Therefore, it would be nec-
pressure allows the proper prediction of sharp fronts in the fractures. essary to first generate a standard CV cell mesh, and later a ho-
mogenization procedure would be required to obtain CV cells with
Introduction constant permeability. The homogenization or upscaling of per-
This work is focused on the numerical treatment of two main meability is somehow possible, but the same is not true for rock-
physical aspects of multiphase flow in fractured porous media: fluid properties; most challenging is capillary pressure with dif-
heterogeneity in rock-fluid properties and reservoir wettability. ferent endpoints. Therefore, the approach with the homogeneous
In a previous work (Monteagudo and Firoozabadi 2004), a CV cell may be suitable for single-phase simulation where rock-
CVDF method was used to discretize the system of equations fluid interactions are not part of the problem. However, rock-fluid
governing water injection in fractured media with strong-water- interactions have to be taken into account for simulation of mul-
wettability state and homogeneous matrix and rock-fluid proper- tiphase flow in fractured porous medium. Frequently, capillary
ties. The method was restricted to a finite contrast in matrix- pressure is disregarded in two-phase flow simulations; however,
fracture capillary pressure. In this work, we extend the CVDF capillary pressure is of importance for simulation of multiphase
model for simulation of water injection in fractured media com- flow in fractured porous media (Monteagudo and Firoozabadi
prised of heterogeneous rocks and wettability conditions from 2004; Karimi-Fard and Firoozabadi 2003). Predictions of flow
strong-water-wetting to mixed-wetting conditions. We also present pattern and oil recovery may be severely affected if capillary pres-
a formulation for infinite contrast in capillary pressures of matrix sure effect is neglected.
and fractures (zero capillary pressure in the fracture and finite With these considerations, we have opted for the use of the
capillary pressure in the matrix). standard CV method with a novel approach for rock-fluid hetero-
The control volume (CV) method, first proposed by Baliga and geneity based on the capillary pressure continuity concept. The
Patankar (1980), is a finite-volume formulation over dual cells numerical inaccuracies of the velocity field at the interface of
(CVs) of a Delaunay mesh. It is locally conservative and suited for heterogeneous media may not be as relevant as the ones generated
unstructured grids. It has been widely employed for the simulation by homogenization of rock-fluid properties or the neglect of cap-
of multiphase flow in porous media (Monteagudo and Firoozabadi illary effects in multiphase flow in fractured media.
2004; Verma 1996; Helmig 1997; Helmig and Huber 1998; Bas- The second focus of this work is proper accounting for reser-
tian et al. 2000; Geiger et al. 2003) and the convergence of the voir wettability. In some reservoirs, wettability covers a wide
method for two-phase immiscible flow in porous medium has al- range from strong water-wetting to strong oil-wetting (Rao et al.
ready been proved (Michel 2003). 1992; Morrow 1990; Robin 2001). Even in a single reservoir,
Numerical treatment of heterogeneity in the framework of the wettability can change from strong water-wetting at the bottom to
CV method has been extensively studied in the past (Edwards oil-wetting at the top of the oil column.
2002; Edwards and Rogers 1998; Prevost 2000; Aavatsmark et al. The paper is outlined as follows: We first present the governing
1998a, b). Nevertheless, those works have focused on absolute equations for two-phase immiscible and incompressible flow for
permeability heterogeneity and anisotropy in single-phase flow. fractured media. Next, we briefly describe the CV discretization of
The main concern in those works is the use of full tensor perme- the governing equations with the numerical treatment for hetero-
ability and the accurate generation of streamlines (required by the geneity in rock-fluid properties and different wettability states.
streamline numerical method). It is well known that the standard Then, we provide numerical simulation of field-scale and lab-scale
CV method produces inaccurate velocity fields around the inter- examples and numerical tests for mixed-wettability conditions in
faces of heterogeneous media as the contrast in permeability is 2- and 3D fractured media. Finally, concluding remarks are made.
Governing Equations
Two-Phase Incompressible Flow in Porous Media. The two-
* Now with ConocoPhillips. phase incompressible, immiscible-flow displacement in porous
Copyright © 2007 Society of Petroleum Engineers
media can be described by two equations (Monteagudo and
Firoozabadi 2004; Aziz and Settari 1979):
Original SPE manuscript received for review 16 July 2005. Revised manuscript received 10
November 2006. Paper (SPE 98108) peer approved 17 November 2006. −ⵜ ⭈ 共␭w + ␭n兲ⵜ⌽w − ⵜ ⭈ ␭nⵜ⌽c − 共qw + qn兲 = 0 . . . . . . . . . . . (1)

September 2007 SPE Journal 355

 Fig. 1—Control volume cell for heterogeneous media. (a) Standard CV method; (b) locally homogeneous CV.

⭸␾Sw CV Discretization
− ⵜ ⭈ ␭wⵜ⌽w − qw = 0, . . . . . . . . . . . . . . . . . . . . . . . . . . . (2)
⭸t We use the CV methodology for the spatial discretization of both
the flow potential and saturation equations. In this method, the
where ␾ is the porosity, ␭i is the mobility of phase i defined by integration is performed over the CV cells, which are constructed
␭i=kkri /␮i , and k is the absolute permeability; kri , ␮i , and qi are the as the median dual of a Delaunay triangulation in 2D or tetrahe-
relative permeability, viscosity and source/sink term of phase i, drization in 3D. Fig. 1a shows the 2D CV cell formed by joining
respectively, ⌽w is the water flow potential, Sw is the water satu- the barycenters of Delaunay triangles with their edge midpoints for
ration, and ⌽c is the capillary flow potential which is defined unfractured media. In Fig. 2, we show the CV cell for fractured
(Karimi-Fard and Firoozabadi 2003) by media with homogeneous matrix (Fig. 2a) and with heterogeneous
⌽c = ⌽n − ⌽w = Pc + 共␳n − ␳w兲gz. . . . . . . . . . . . . . . . . . . . . . . . (3) matrix (Fig. 2b) media.
Eqs. 1 and 2 have three types of terms:
In Eq. 3, Pc is the capillary pressure, ␳i is the density of phase i, g 1. Time derivative, ⭸␾Sw /⭸t.
is the acceleration of gravity, and z is the vertical coordinate, 2. Source/sink, qi for i⳱{n, w}.
positive in the upward direction. 3. Divergent term ⵜ·F, where the flux vector F⳱F(Sw) may be
Eq. 1 is referred to as the flow-potential equation, which is either (␭w+␭n)ⵜ⌽w and ␭nⵜ⌽c in Eq. 1, or ␭wⵜ⌽w in Eq. 2.
elliptic in nature, while Eq. 2 is referred to as the saturation equa- All mobilities are computed at the boundaries of the CV cell
tion, which can be seen as a diffusion-convection equation (Peace- with the proper upwinding, based on the flow direction at the
man 1977). The main variables are the wetting-phase flow poten- boundary. Next, we provide the discretization of all these terms for
tial, ⌽w , and saturation, Sw. In most of our numerical simulations, 2D fractured media with homogeneous and heterogeneous matrix.
the boundary conditions are assumed to be impervious: Extension to 3D media is straightforward.
vi ≡−␭iⵜ⌽i⳱0 for both phases i⳱{n, w}, where vi is the velocity
vector of phase i. We have also incorporated the Dirichlet bound-
Discretization for Fractured Media With Homogeneous Ma-
ary condition to include an aquifer in a field-scale problem that
trix. In the following, we first add the matrix and fracture terms
will be simulated in this work. In that case, ⌽w⳱⌽D and Sw⳱1 at
and then approximate the integration in the finite-volume.
the bottom of the reservoir.
Time Derivative.

兰 ⭸t 共␾S 兲dA ≈ 冉 A ␾ 冊
Formulation for Fractured Media. Details of the formulation for
⭸ ⭸S wm ⭸S wf
fractured media are provided in Monteagudo and Firoozabadi m m
+ Af␾ f , . . . . . . . . . . . . . . (7)
(2004). The most remarkable feature of the formulation is that no
w
⭸t ⭸t
A
computation of matrix-fracture transfer term is required. This is
possible because of the integration in the CV cell and the use of the where superscripts m and f denote matrix and fracture variables.
crossflow equilibrium concept, which leads to capillary continuity Eq. 7 represents the change in the volume of the wetting phase in
(Firoozabadi and Hauge 1990; van Duijn et al. 1994). Using the the volume element per unit time. It consists of two terms, one
crossflow equilibrium, we write the equality of flow potentials at from the matrix contribution and the other from the fracture con-
the interface between the matrix and fractures: tribution. If there is no fracture in the volume elements, then the
contribution is from the matrix only. Note that in 2D, the product
⌽ im = ⌽ if for i = 兵n, w其. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (4) Am␾m (Am is the area of the matrix in the element; in 3D it should
be replaced by the volume of the matrix in the element) and Af␾ f
From Eq. 4, we readily obtain
(Af is the area of the fracture in the element) provide the volume of
⌽ cm共S wm 兲 = ⌽ cf 共S wf 兲. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (5) the matrix and fracture in the volume element, respectively. We
can express Eq. 7 in terms of S m w by using Eq. 6 and the chain rule:

冉 冊
Eq. 5 implies that water saturation may be discontinuous across the
⭸ dS wf ⭸S wm
matrix-fracture interface with different matrix and fracture capil-
lary pressure functions. A clear physical relationship is established 兰 ⭸t ␾S dA ≈
w Am␾m + Af␾ f
dS wm ⭸t
. . . . . . . . . . . . . . . . (8)
between S m f
w and S w :
A

Source/Sink Term.
S wf = 关P cf 兴−1P cm共S wm兲. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (6)
Therefore, the system of Eqs. 1 and 2 can be integrated with the
CV method in terms of matrix variables only.
兰q dA≈A q
A
w
m m
w + Afq wf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (9)

356 September 2007 SPE Journal

 Fig. 2—Control volume cell containing a fracture. (a) Homogeneous matrix; (b) heterogeneous matrix.

Divergence Term. Source/Sink Term. This term can be written as

nm

兰q dA ≈ 兺 A q
nb

兰ⵜ · FdA = 兰F · n d⌫ ≈ 兺 关F · ne兴 .
A ⌫A
A
j=1
j . . . . . . . . . . . . . . . (10)
A
w
k=1
k k
w. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (14)

Divergence Term. Integration of the divergence term is pre-

The approximation in Eq. 10 allows us to relate the rate of change
sented in Eq. 10. Each triangle contained in the CV cell is locally
in the volume of the wetting phase in the volume element in terms
homogeneous, and the saturation of different media can be related
of the normal fluxes across various edges (boundary elements) j of
through Eq. 6.
measure ej (in 3D we need to use the edge surface area). The
number of boundary elements is nb. The summation over all the Formulation Using Normalized Saturations. Normalized satu-
boundary elements provides the net flux of the wetting phase over ration is defined as
the control volume. We remark that for the CV cell containing a
fracture, the summation in Eq. 10 contains also the 1D flow along Sw − Swi
the fracture direction. Sw = , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (15)
1 − Swi − Sor
A detailed derivation of the system of equations for the IMPES
formulation can be found in Monteagudo and Firoozabadi (2004). where Swi is the irreducible-water saturation and Sor is the residual-
oil saturation.
Discretization for Fractured Media With Heterogeneous Ma- By using the chain rule, Eq. 2 can be expressed as:
trix. Fig. 2b shows a CV cell with a fracture and two different rock ⭸␾Sw
matrices in its interior. Similar to the matrix-fracture interface, we ␣ − ⵜ ⭈ ␭w共Sw兲ⵜ⌽w − qw = 0, . . . . . . . . . . . . . . . . . . . . . (16)
invoke capillary continuity at the interface between different rocks. ⭸t
Therefore, we can relate the saturation of all rocks at the interface where ␣⳱1−Swi−Sor
of a heterogeneous medium.
Time Derivative. The integration of the time-derivative term is Numerical Treatment for Large Contrast
approximated by extending the idea of Eqs. 7 and 8. Between Matrix And Fracture
Capillary Pressures
nm
⭸ ⭸
兰 ⭸t 共␾S 兲dA ≈ 兺 ⭸t S ␾ A ,
A
w
k=1
k
w
k k
. . . . . . . . . . . . . . . . . . . . . . . (11)
In our previous work (Monteagudo and Firoozabadi 2004), the
contrast between matrix and fracture capillary pressure was ex-
pressed by the ratio Bm /Bf using the capillary pressure model in
where nm is the number of different media inside the CV cell. Eq. Eq. 12. All the simulations in that work were carried out for a finite
11 may be expressed in terms of the saturation of the reference ratio, that is, Bf >0. When Bm⳱Bf⳱0, the model of Monteagudo
and Firoozabadi (2004) is applicable because S m w ⳱S w and dS w /
f f
medium S +w. The reference medium can be chosen on the basis of
parameter B in the capillary pressure model of the form dS m
w ⳱1. There is a numerical problem when B m >0 and B f⳱0 or
when Bm /Bf is large; with strong water wetting and relatively thick
P kc = −B k lnS wk , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (12) fractures, Bm /Bf would be large. Below we present a solution to
this problem.
where superscript k⳱1 . . . nm is an index for the different media Fig. 3 illustrates the imbibition process for the fracture and the
inside the control volume. Note that by using the capillary conti- matrix media. Fig. 3a shows the fracture and matrix capillary
nuity concept, a relationship can be established between the dif- pressure curves, for Bf⳱0.6 atm and Bm⳱1.0 atm, respectively
ferent media saturations. For example, in the case of fractured (Bm /Bf⳱1.67). Fig. 3b depicts the fracture and matrix capillary
media: S mw /S w⳱exp(−B /B ).
f m f pressure for Bf⳱0.1 atm and Bm⳱1.0 atm, respectively (Bm /
In this work, we select the reference medium with the highest Bf⳱10). Imbibition of a CV cell containing a matrix-fracture in-
value of Bk in Eq. 12. terface with these capillary pressure curves will follow the path
From Eq. 6 and the chain rule we derive: through the points A, B, and C as indicated by the arrows, provided

冋 册
the initial water saturation is zero. Because capillary continuity
nm
⭸ dS wk ⭸ + expressed by Eq. 6 must be satisfied at the interface during the
兰 ⭸t 共␾S 兲dA ≈ 兺 dS
A
w
k=1
+
w
␾kAk S . . . . . . . . . . . . . . . . (13)
⭸t w
imbibition, S fw⳱0 and dS fw /dS m
w ⳱0 along the path from A to B for
capillary pressures depicted in Figs 3a and 3b. On the other hand,

September 2007 SPE Journal 357

Fig. 3—Imbibition paths in fracture (left column) and in matrix (right column) at the interface. (a) Bf =0.6 atm, and Bm=1.0 atm; (b)
Bf =0.1 atm, and Bm=1.0 atm.

f m
along the path from B to C the derivative dS w /dS w will be If Pc>␧, then
much larger when Bf ⳱0.1 atm (Fig. 3b) than when Bf ⳱0.6 atm
⭸ ⭸S wm
(Fig. 3a). From this analysis, one can observe that in the event of
Bf⳱0 then Bm /Bf→⬁ and there is apparently no way to relate S m w
兰A
⭸t
共␾Sw兲dA ≈ 共Am␾m兲
⭸t
. . . . . . . . . . . . . . . . . . . . . . . . . . (18)
and S fw along the path B to C with the formulation outlined pre-
viously because dS w /dS w ⳱⬁. To avoid this numerical problem,
f m
Otherwise,
the left term in Eq. 8 may be expressed in an alternative way:
⭸ ⭸S wf
兰 ⭸t 共␾S 兲dA ≈ 共A ␾ 兲
冉 冊 冉 冊
f f
, . . . . . . . . . . . . . . . . . . . . . . . . . . (19)
dS wf ⭸S wm dS wm ⭸S wf
w
⭸t
A
Am␾m + Af␾ f ≡ Am␾m + Af␾ f . . . (17)
dS wm ⭸t dS wf ⭸t where ␧ is a tolerance.

Then we have two equivalent expressions for the left term in Eq. Numerical Examples
8 that can be distinctly used along different sections of the imbi- In this section, we present several numerical examples. First, we
bition path shown in Fig 3b. If chosen correctly, the numerical simulate oil production from data extracted from a large fractured
problem described here is avoided. Eq. 17 also provides a physical reservoir. The second example covers data from a laboratory test
interpretation of the imbibition process inside the CV cell; for the of water injection in a stack of matrix blocks at different injection
w ⳱0, which implies imbibition in the matrix
path A-B, dS fw /dS m rates. Next, we present the results of numerical simulation for
alone, and for the path B-C dS mw /dS w⳱0, mainly only the fracture
f
water injection in 2D and 3D for mixed-wettability state. The 2D
saturation changes. We have implemented this criterion in our and 3D meshes were generated with the packages triangle (Shew-
model to switch between the two equivalent expressions in Eq. 17. chuk 1996) and tetgen (Si 2002), respectively. All simulations
That is, were run with a Pentium PC 2.66 GHz.

358 September 2007 SPE Journal

Fig. 4—Cross-sectional study. (a) Horizontal view; (b) XZ cross section between Wells 13V and 6V showing horizontal Well 6H.

Field-Scale Example With Horizontal Production Well. A hori-

zontal view of the reservoir simulated is shown in Fig. 4a. A
region with an areal dimension of 1200×500m (shown as a rect-
angle with dashed lines) was selected for simulation. The region
contains two vertical wells, 13V and 6V, and one horizontal well,
6H. Only the horizontal well is producing. The region contains
several faults (the dark solid line). In Fig. 4a, only the three main
faults (the grey solid lines) that cross the region are shown. The
horizontal well is contained in one of the main faults. Fig. 4b
shows the fault plane XZ joining wells 13V and 6V with details of
the location of the horizontal well. We use a simple approach to
represent wells by a sink term and assign a triangular flow distri-
bution along the horizontal well (Joshi 1991) (see Fig. 5). The
reservoir has a thickness of 424 m with an aquifer at the bottom. Fig. 5—Triangular flow distribution in a horizontal well.

September 2007 SPE Journal 359

 Fig. 6—Conforming tetrahedrization mesh used in the simulation of field-scale problem. (a) Fracture mesh; (b) matrix mesh.

The aquifer is treated as the Dirichlet boundary, where water flow Lab-Scale Example for a Stack of Matrix Blocks. We simulate
potential and saturation are specified. the displacement of n-decane by water from a stack of chalk matrix
Fig. 6 shows the conforming tetrahedrization in the fractured blocks. Description of the equipment and experimental data are
region. The horizontal well 6H is shown as a dashed thick line. We provided by Pooladi-Darvish and Firoozabadi (2001). Previous
have incorporated three main and seven minor fractures as quad- attempts to model the waterflooding have not been successful
rilaterals (see Table 1 for vertex coordinates). Fracture thickness when one focuses on the water/oil front in the fractures (Pooladi-
was set to 0.1 cm. The mesh contains 4,000 nodes. Preliminary Darvish and Firoozabadi 2001; Terez and Firoozabadi 1999).
sensitivity study showed that with this refinement, there is very In the setup, 12 matrix blocks were assembled; fractures are the
little change in oil-recovery predictions. Fluid and rock properties space between blocks and between blocks and the container walls.
are shown in Tables 2 and 3, respectively. The reservoir is lay- The 3D conforming mesh for the porous media with 1,800 nodes
ered. Fig. 7 shows the permeability distribution. Rock-fluid inter- is presented in Fig. 11 (note that there are three matrix blocks per
actions are shown in Fig. 8. Note that the relative permeability and row). Tables 4 and 5 present fluid and rock properties, respec-
capillary pressure curves may not show a consistent wettability state. tively. Rock-fluid interactions parameters are shown in Table 6.
The initial distribution of fluids is determined using the vertical The capillary pressure model used is defined by Eq. 12, and the
equilibrium criterion (Pooladi-Darvish and Firoozabadi 2001): relative permeabilities are modeled with power laws:
Pc共z兲 = Pco − g共␳n − ␳w兲z , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (20) krw = k rw
0
S wnw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (21)
where Pco is the capillary pressure at the bottom of the reservoir krn = k0rn共1 − Sw兲nn, . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (22)
(z⳱0) and z is the vertical coordinate (positive in the upward
direction). Once Pc is known, Sw is computed from Eq. 12. where Sw is the normalized saturation defined by Eq. 15. All the
The simulation was carried out to 140 years of production, rock-fluid interaction parameters are reported in Pooladi-Darvish
keeping a constant total flow rate of 2400 B/D. The total CPU time and Firoozabadi (2001). Fig. 12 depicts the capillary pressure and
for this simulation was 12 hours. Water saturation contours at relative permeability. Note that capillary pressure in the fracture is
time⳱70 years are shown in Fig. 9. The corresponding simulated considered zero, as was evidenced by observation of a flat oil/
oil recovery curve is shown in Fig. 10. The initial reservoir pres- water interface in the fractures. The relative permeability for the
sure was 2,500 psia at z⳱0. The predicted average pressure across matrix shows a strong-water-wetting state.
the horizontal well after 17 years of production was approximately Using the previously described data (without any adjustment),
1.500 psia. The reported data for the same time period is 1,550 we simulated the experiments. There is good agreement between
psia. Our preliminary results are on the order of the reported data the simulation results and measured recovery data for all three

360 September 2007 SPE Journal

injection rates. Results are shown in Fig. 13. Our simulation results
also show the correct shape of the front in the fractures (see
Fig. 14). The front based on our simulation is in better agreement that may have very different rock-fluid properties with different
with measured data than in the previous works (Pooladi-Darvish endpoints. The idea is tested in the field-scale simulation of a
and Firoozabadi 2001; Terez and Firoozabadi 1999). One could layered reservoir.
obtain a better match by adjusting the matrix capillary pressure 3. A novel numerical procedure has been proposed when there is
curve slightly. The CPU time required for a displacement of 2 PV large contrast in capillary pressure at the matrix-fracture inter-
at each flow rate was 27 hours. Flow through fractures for Pc<␧ face (P fc≈0). The procedure was tested in the simulation of the
requires a small timestep to avoid instabilities. This results in the waterflooding in a stack of matrix blocks. The method not only
increased CPU time. provides proper prediction of oil recovery at different water
injection flow rates, but it also predicts a sharp water front in
Mixed-Wettability Examples. We present 2- and 3D numerical the fractures.
simulations to evaluate the performance of the CVDF method to 4. We also present the results of the numerical simulation for water
simulate waterflooding in porous media with mixed wettability injection in a mixed-wet medium using our CVDF model. In the
state. Because of the way the problem has been formulated, the past, the discrete fracture approach had not been used for water
mixed-wet state has not been implemented for two-phase immis- injection in fractured media with a mixed-wet state that has a
cible flow in fractured media (Bogdanov et al. 2003). higher nonlinearity than a water-wet state. Results show that our
The 2D example consists of a square of 25×25 m containing 6 formulation performs well in the mixed-wet state.
fractures. The 2D mesh with 2,700 nodes is shown in Fig. 15, with
the fractures marked as thick lines. Water is injected at a flow rate Nomenclature
of 3.96×10−7 m3/s (2.7×10−4 PV/day) at the lower left corner, and
the producing well is at the opposite corner. The initial water A ⳱ area, L2, m2
saturation is zero in the whole domain. B ⳱ parameter for Pc model Eq. 12, m/Lt2, Pa [atm]
The 3D example consists of a cube of 20×20×20 m with one e ⳱ measure of a boundary segment of a CV cell
fracture plane inside. A plot of the mesh for the fracture plane and F ⳱ generic flux crossing the boundaries of a CV cell
the matrix is shown in Fig. 16. Water is injected at a flow rate of k ⳱ absolute permeability, L2, m2
3.7×10−4 m3/s (0.02 PV/day) at coordinates [0 m, 0 m, 0 m], and kri ⳱ relative permeability of phase i⳱{n, w}
the producing well is placed at coordinates [20 m, 20 m, 20m]. nb ⳱number of boundary segments of a CV cell
Tables 7 and 8 present fluid and rock properties used for the nm ⳱number of different media inside a CV cell
examples. Rock-fluid interactions are shown in Fig. 17. Capillary
nn ⳱ exponent for the nonwetting-phase relative
pressure has two branches, positive and negative, indicating mixed
wettability state. These rock-fluid properties are used for both 2D permeability model Eq. 21
and 3D simulations. Because the capillary pressure model in Eq. nw ⳱ exponent for the wetting-phase relative permeability
12 is only suitable for a water-wet media, we compute dS fw /dS m w
model Eq. 22
numerically for the mixed-wetting state. Water saturation is ini- n ⳱ outward surface normal to the boundary of a CV cell
tialized using the positive Pc–Sw relationship shown in Fig. 17. Pc ⳱ capillary pressure, m/Lt2, Pa [atm]
The CPU time for the 2D geometry was 24 hours for a 1-PV qi ⳱ sink/source term of phase i, t−1, s−1
displacement. This represents a speed-up of three times when com- Si ⳱ saturation of phase i={n, w}
pared with a previous code (Karimi-Fard and Firoozabadi 2003) Sor ⳱ residual oil saturation
due to our formulation. For the 3D geometry, the CPU time was 20
minutes. Water saturation contours at 25 and 50% PV displace-
ment for the 2D configuration are shown in Fig. 18a and b, re-
spectively. For the 3D configuration, we show in Fig. 19a and b
the water saturation contours at 20% PV displacement for fracture
and matrix, respectively. It is observed that the water front in the
fracture moves fast in this problem when compared to the contours
obtained at the faults in the field-scale reservoir configuration.
These results reflect what we found in our previous work (Mon-
teagudo and Firoozabadi 2004); the contrast in capillary pressure
affects the performance of waterflooding in fractured media. Be-
cause of large nonlinearity of capillary pressure in mixed-wet me-
dia, CPU time is much higher than in water-wet media.

Conclusions
1. We have extended the control-volume discrete-fracture method
(Monteagudo and Firoozabadi 2004) to account for important
physical aspects in reservoir simulation, such as heterogeneity in
matrix rock and rock-fluid properties, different wettability con-
ditions, and large contrast in the matrix-fracture capillary pressures.
2. We propose the use of capillary continuity concept in order to
solve the problem of having a heterogeneous CV cell with rocks Fig. 7—Permeability field of the layered reservoir.

September 2007 SPE Journal 361

 Fig. 8—Field-scale simulation (rock-fluid properties). (a) Capillary pressure; (b) relative permeability.

Fig. 9—Water saturation contours after 70 years of production. (a) Fracture contours; (b) matrix contours.

Swi ⳱
irreducible water saturation
Si ⳱
normalized saturation of phase i⳱{n, w} (Eq. 15)
t ⳱
time, s
␣ ⳱
parameter in Eq. 16
⌫A ⳱
boundary of CV cell A
␧ ⳱
tolerance
␭i ⳱
mobility of phase i, L3t/m, m2/Pa-s
␭t ⳱
total mobility, L3t/m, m2/Pa-s
␮i ⳱
viscosity of phase i, m/Lt, Pa-s
␳i ⳱
density of phase i, m/L3, kg/m3
␾ ⳱
porosity
⌽c ⳱
capillary potential, m/Lt2, Pa
⌽D ⳱
flow potential specified at Dirichlet boundary, m/Lt2,
Pa
⌽i ⳱ flow potential of phase i, m/Lt2, Pa

Subscripts
n ⳱ nonwetting phase
Fig. 10—Oil recovery for the field-scale simulation. w ⳱ wetting phase

362 September 2007 SPE Journal

 merical Treatment of Multiphase Flows in Porous Media, ed. Chen,
Ewing, and Shi, 1–18. Berlin: Springer.
Bogdanov, I.I., Mourzenko, V.V., Thovert, J.-F., and Adler, P.M. 2003.
Two-phase flow through fractured porous media. Physical Review E 68
(2): 1–24. DOI: 10.1103/PhysRevE.68.026703.
Durlofsky, L.J. 1994. Accuracy of mixed and control volume finite element
approximations to Darcy velocity and related quantities. Water Re-
sources Research 30 (4): 965–973.
Edwards, M.G. 2002. Unstructured, control-volume distributed, full-tensor
finite-volume schemes with flow based grids. Computational Geo-
sciences 6: 433–452.
Edwards, M.G and Rogers, C.F. 1998. Finite volume discretization with
imposed flux continuity for the general tensor pressure equation. Com-
putational Geosciences 2: 259–290.
Firoozabadi, A. and Hauge, J. 1990. Capillary Pressure in Fractured Porous
Media. JPT 42 (6): 784–791; Trans., AIME, 289. SPE-18747-PA. DOI:
10.2118/18747-PA.
Geiger, S., Roberts, S., Matthai, S., and Zoppou, C. 2003. Combining finite
volume and finite element methods to simulate fluid flow in geological
media. ANZIAM Journal 44E: C180–C201.
Helmig, R. 1997. Multiphase Flow and Transport Processes in the Sub-
surface. First edition. Berlin: Springer.
Fig. 11—Conforming tetrahedrization mesh for lab-scale simu- Helmig, R. and Huber, R. 1998. Comparison of Galerkin-type discretiza-
lation. (a) Fracture mesh; (b) matrix mesh. tion techniques for two-phase flow in heterogeneous porous media.
Advances in Water Resources 21 (8): 697–711.
Joshi, S.D. 1991. Horizontal Well Technology. Tulsa: Penwell Publishing.
Superscripts
Karimi-Fard, M. and Firoozabadi, A. 2003. Numerical Simulation of Water
f ⳱ fracture Injection in Fractured Media Using the Discrete-Fracture Model and
m ⳱ matrix the Galerkin Method. SPERE 6 (2): 117–126. SPE-83633-PA. DOI:
10.2118/83633-PA.
Acknowledgments Michel, A. 2003. A finite volume scheme for two-phase immiscible flow
in porous media. SIAM J. Numer. Anal. 41 (4): 1301–1317.
This work was supported by the member companies of the Res-
Monteagudo, J.E.P. and Firoozabadi, A. 2004. Control-volume method for
ervoir Engineering Research Institute (RERI).
numerical simulation of two-phase immiscible flow in 2-D and 3-D
discrete-fracture media. Water Resources Research 40. DOI: 10.1029/
References 2003WR002996.
Aavatsmark, I., Barkve, T., Boe, O., and Mannseth, T. 1998a. Discretiza- Morrow, N.R. 1990. Wettability and its Effect on Oil Recovery. JPT 42
tion on unstructured grids for inhomogeneous, anisotropic media. Part (12): 1476–1484; Trans., AIME, 289. SPE-21621-PA. DOI: 10.2118/
I: Derivation of the methods. SIAM J. Sci. Comput. 19 (5): 1700–1716. 21621-PA.
Aavatsmark, I., Barkve, T., Boe, O., and Mannseth, T. 1998b. Discretiza- Peaceman, D. 1977. Fundamentals of Numerical Reservoir Simulation.
tion on unstructured grids for inhomogeneous anisotropic media. Part New York: Elsevier.
II: Discussion and numerical results. SIAM J. Sci. Comput. 19 (5):
1717–1736.
Aziz, K. and Settari, A. 1979. Petroleum Reservoir Simulation. London:
Applied Science Publishers.
Baliga, B. and Patankar, S. 1980. A new finite-element formulation for
convection-diffusion problems. Numerical Heat Transfer 3: 393–409.
Bastian, P., Helmig, R., Jakobs, H., and Reichenberger, V. 2000. Numeri-
cal simulation of multiphase flow in fractured porous media. In Nu-

September 2007 SPE Journal 363

 Fig. 12—Lab-scale waterflooding (rock-fluid properties). (a) Capillary pressure; (b) relative permeability.

Pooladi-Darvish, M. and Firoozabadi, A. 2001. Experiments and Model-

ling of Water Injection in Water-wet Fractured Porous Media. Journal
of Canadian Petroleum Technology 39 (3): 31–42.
Prevost, M. 2000. The streamline method for unstructured grids. MSc
thesis. Stanford, California: Stanford University.
Rao, D.N., Girard, M., and Sayegh, S. 1992. The influence of reservoir
wettability on waterflood and miscible flood performance. Journal of
Canadian Petroleum Technology 31 (6): 47–54.
Robin, M. 2001. Interfacila phenomena: Reservoir wettability in oil recov-
ery. Oil & Gas Science and Technology-Rev. IFP 56 (1): 55–62.
Shewchuk, J. 1996. Triangle: Engineering a 2D Quality Megenerator and
Delaunay Triangulator. In First Workshop on Applied Computational
Geometry, 124–133. Philadelphia: Assoc. for Comput. Mach.
Si, H. 2002. Tetgen. A 3-D Delaunay Tetrahedral Mesh Generator. v.1.2
Users Manual. Technical Report 4. Berlin: Weierstrass Institute for
Applied Analysis and Stochastics.
Terez, I.E. and Firoozabadi, A. 1999. Water Injection in Water-Wet Frac-
tured Porous Media: Experiments and a New Model With Modified
Buckley-Leverett Theory. SPEJ 4 (2): 134–141. SPE-56854-PA. DOI:
Fig. 13—Oil recovery for different water injection rates in the
stack of matrix blocks. 10.2118/56854-PA.

Fig. 14—Simulated and experimental fracture waterfront at different PV displacement in the stack of matrix blocks. (a) Front level
in the fracture; (b) water saturation in the fractures.

364 September 2007 SPE Journal

 van Duijn, C., Molenaar, J., and de Neef, M. 1994. The effect of capillary
forces on immiscible two-phase flow in heterogeneous porous media.
Tech. Report 94-103. Delft, The Netherlands: Delft University of Tech-
nology.
Verma, S. 1996. Flexible Grids for Reservoir Simulation. PhD dissertation.
Stanford, California: Stanford University.
Willhite, P. 1986. Waterflooding. Second edition. Richardson, Texas: So-
ciety of Petroleum Engineers.

Jorge Palomino Monteagudo is a staff reservoir engineer in the

Reservoir Simulation Development Division of ConocoPhillips,
Houston. He holds a BS degree from the Universidad Nacional
de Ingenieria in Lima, Peru, and MS and DSc degrees from the
Universidade Federal do Rio de Janeiro in Brazil, all in chemical
engineering. He is currently a Technical Editor for SPEREE. Ab-
bas Firoozabadi is a senior scientist and director at RERI. He also
teaches at Yale U. and at Imperial College, London. e-mail:
[email protected]. His main research activities center on thermo-
dynamics of hydrocarbon reservoirs and production and on
multiphase-multicomponent flow in fractured petroleum reser-
voirs. Firoozabadi holds a BS degree from Abadan Inst. of Tech-
nology, Abadan, Iran, and MS and PhD degrees from the Illi-
nois Inst. of Technology, Chicago, all in gas engineering.
Fig. 15—Conforming triangulation mesh for 2D mixed-wet frac- Firoozabadi is the recipient of both the 2002 SPE Anthony Lucas
tured media. Gold Medal and the 2004 SPE John Franklin Carll Award.

Fig. 16—Conforming tetrahedrization mesh for 3D mixed-wet fractured media. (a) Fracture; (b) matrix.

Fig. 17—Mixed-wet fractured system (rock-fluid properties). (a) Capillary pressure; (b) relative permeability.

September 2007 SPE Journal 365

Fig. 18—Water saturation contours for injection into a 2D mixed-wet fractured media with initial Sw=0. (a) 25% PV displacement;
(b) 50% PV displacement.

Fig. 19—Water saturation contours at 20% PV displacement for injection into a 3D mixed-wet system where initial Sw field is
computed with vertical equilibrium concept. (a) Fracture contours; (b) matrix contours.

366 September 2007 SPE Journal