Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
TwoPhase Flow Simulations in a Natural Rock Fracture using the VOF Method
D. Crandall1, G. Ahmadi2, D.H. Smith3, and G. Bromhal4
1URS Washington Division, National Energy Technology Laboratory
Morgantown, WV 26507 USA
Dustin. Crandall@ur.netl.doe.gov
2Mechanical and Aeronautical Engineering Department, Clarkson University
Potsdam, NY 13699, USA
gahmadi@tclarkson.edu
National Energy Technology Laboratory
Morgantown, WV 26507 USA
Duane. Smith @netl.doe.gov
4National Energy Technology Laboratory
Morgantown, WV 26507 USA
Grant.Bromhal@netl.doe.gov
Keywords: Geological, relative permeability, CO2 sequestration
Abstract
Geological sequestration of carbon dioxide (CO2) is the process of placing CO2 into subsurface formations in such a way that
it will remain permanently stored. In the United States deep brine aquifers potentially provide the largest storage capability for
geologically sequestered CO2 (NETL, 2008). Within the relatively lowpermeability rocks that form brine aquifers, fractures
exist that can act as natural fluid conduits. Understanding how CO2 is transported through initially liquid saturated rock
fractures is important for the prediction of storage capacity and monitoring, verification and accounting of CO2 in potential
sequestration locations.
Our computational study examined the motion of immiscible fluids as they were transported through models of a fracture in
Berea sandstone. The natural fracture geometry was initially scanned using microcomputerized tomography (CT) at a fine
volumepixel (voxel) resolution by Karpyn et al. (2007). This CT scanned fracture was converted into a numerical mesh for
twophase flow calculations using the finitevolume solver FLUE N T and the volumeoffluid method. Simulations of single
and twophase flows in the reconstructed fracture geometry were performed. Initial twophase simulations of air injection into
a water saturated fracture were performed and compared to experimental results from a bench scale model manufactured using
stereolithography. In both experiments and simulations the invading air moved intermittently, rapidly invading largeaperture
regions of the fracture (Crandall et al. 2009). Relative permeability curves were developed to describe the motion of the two
fluids. These permeability curves can be used in reservoirscale discrete fracture models for predictions of multicomponent
fluid motion within fractured geological formations. The numerical model was then upgraded to better mimic the subsurface
conditions at which CO2 will move into brine saturated fractures. The simulation results show that, with the predicted lower
interfacial tension between subsurface fluids, an increased volume of the lessviscous CO2 fills the rough fracture geometry.
Introduction from millimeters to miles. The high permeability (k) of open
fractures allows for enhanced fluid mobility and greater
Fluid movement within fractured geological media is of fluid flow rates. Thus, a thorough understanding of the
importance in subsurface activities including oil recovery behavior of fluid flow through fractured porous media will
(Selley 1998), CO2 sequestration (Klusman 2003), enable more accurate prediction of fluid motion on the
geothermal energy extraction (Hanano 2004), and nuclear reservoir scale.
waste disposal predictions (Selroos et al. 2002). Fractures The inclusion of fractures in reservoir scale simulators
exist in nearly all reservoir rock formations, ranging in size requires approximations of some sort. The continuum
Paper No
approach to reservoir simulation assumes that the fractures
and unfractured zones of porous media can be coupled
through some relationship, such as continuity of pressure
between the fracture walls and the porous media (Wu et al.
2004). Discrete fracture and channel network simulators,
which explicitly model flow through fracture networks,
often use a Darcy's Law relationship to model fluid flow
through the fractures (Seleroos et al. 2002).
Fractures may contribute to enhanced flow over regions
miles in length, especially when the flow in connected
fracture networks is considered (McKoy and Sams 1997).
The open apertures (b) of fractures are quite small; a
fracture with a mean b of 1 mm is considered large. Fracture
b are defined in this study as the perpendicular distance
between rock surfaces across the open fracture. The
longnarrow geometry of fractures has prompted the use of
"cubiclaw" relations to describe flows with discrete
fracture models. The cubic law is derived from the solution
of the NavierStokes equations for momentum in a wide,
narrow aperture. Simple descriptions of how fluids move
through fractures are a necessary requirement for efficient
reservoir scale discrete fracture models to function.
At the corescale (i.e. cm) fractures are highly
heterogeneous. Consisting of two rough surfaces in a
naturally heterogeneous material, the b between walls varies
throughout the fracture. Regions of zeroaperture (touching
fracture walls) are common within natural fractures
(Zimmerman et al. 1992, Karpyn et al. 2007). Apertures
have been shown to vary significantly within fractures
(Rcnsilul et al. 2000, Amadei and Illangasekare 1994,
Dougan et al. 2000). The relative motion of fluids through
individual natural fractures has been shown to be tortuous
(Tsang 1984). Tortuosity (0) of fractures has been defined as
the additional distance a fluid particle must travel to
circumnavigate through the heterogeneous fracture
geometry. Flow through fractures has also been shown to
primarily occur within 'channels' (PyrakNolte et al. 1990,
Brown et al. 1998), or regions of higher than average flow.
Twophase flow in fractures is affected by the tortuous fluid
path and high flow rates associated with constricted regions.
By identifying the parameters which cause these
irregularities in the fluid flow the present study hopes to aid
in the development of relationships which can be accurately
used in discrete fracture models of twophase reservoir
flows.
Experimental studies (Brown 1989, Hughes and Blunt 2001,
Zimmermann et al. 1992) have shown the channeling (or
tortuous nature of the flow) which occurs within single
natural fractures. CT scans of fractured sandstone have
shown the distribution of multiple fluids within fracture to
be disconnected and dependant on the fracture geometry
(RangelGerman et al. 1999, Karpyn et al. 2007). Other
studies have shown the relative permeability (k,) of
fractured porous media for a variety of fluid/fluid twophase
studies (Berkowitz 2002). The k, of flow in an idealized
fracture model was shown to be highly dependent on the
correlation of aperture heights within fractures (Pruess and
Tsang 1990). An experimental study through a single
fracture has been conducted and compared to a localized
cubic law model, which was shown to over predict the flow
rate experimentally observed (Konzuk and Kueper 2004).
Numerical studies have been reported, and for the most part
corroborate the experimental observations. Porelevel
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
models (Section 2.1.3) have been used to show the motion
of multiple fluids within idealized fracture models (Glass et
al. 1998, Glass et al. 2001, Hughes and Blunt 2001) under
various conditions. A percolation model was used to show
the channeling affect observed in experimental studies,
within a correlated aperture model (PyrakNolte et al. 1990).
The dispersion of tracers was shown to be affected by
localized asperities within fracture (Roux et al. 1998),
which was related to the properties of a macroscopic
(idealized fracture) with a localized cubiclaw model. Finite
difference simulations of flow through realistic fracture
models have shown that alteration of the fracture surface
(smoothing of asperities) results in a more homogenous
fluid invasion front (Lespinasse and Sausse 2000). All of
these works have shown behavior different than that
observed within parallelplates, but no unifying relationship
that accurately describes the flow of fluids within fractured
porous media has gained widespread acceptance.
We present full NavierStokes simulations of single and
multiphase flows through a reconstructed CT scanned
fracture. We show agreement with previously reported 0
and k values and report a novel k, relationship for airwater
flows.
Nomenclature
A Fracture crosssectional area (m2)
b Fracture aperture (m)
k Permeability (m2)
k, Relative permeability ()
L Fracture length (m)
L, Particle i path length (m)
P Pressure (Pa)
Q Volumetric flow rate (m3s')
W Fracture width (m)
Greek letters
0 Tortuosity ()
p Fluid density (kgm3)
P Absolute viscosity (Pas)
Numerical Model
The fracture used in this study was created within a 2.59 cm
diameter and 9.2 cm long core of Berea sandstone, using a
modified Brazilian technique and CT scanned at Penn State
University (Karypn et al. 2007). CT scanning of this
fracture was performed to determine the fracture properties
and to conduct flow experiments. A 240 [tm voxel
(volumepixel) data set was obtained from these scans.
The 240 [tm voxel data was 'cleaned' using an inhouse
code, which removed voids within the rock that were
unconnected to the rest of the fracture. The cleaned
threedimensional fracture geometry is shown in Figure 1. A
grayscale aperture map is used to show the variation in the
b throughout the domain. A profile of the fracture through
the midplane is shown in Figure 1 as well, with the aperture
height increased to illustrate the rough fracture walls. The
asscanned fracture volume was 1216.37 mm3 and the
cleaned fracture volume was 1214.45 mm3; a reduction of
less than one percent between the two models.
Paper No
Heigh[ J '0b 04 S h* .% 'I1#
Figurel: Fracture aperture (b) map and roughwalled
fracture profile.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
the extra distance that a fluid particle will travel as it
traverses the macroscopic length of a fracture. This
Lagrangian path is longer than the fracture length due to
restrictions within the fracture. Thus, 0 is a measurement of
the heterogeneity of a fracture. To measure this with the
numerical model a low injection velocity (0.1 "m/") was
applied to a singlephase flow of air through the fracture
and the steady state solution was obtained. Submicron
particles were released from the constant velocity inlet and
the particles were tracked using the stepbystep particle
tracker in FLUENTTM (2005). Brownian motion and
gravitational forces were not included, so that the particles
only followed the fluid flow. These stepbystep particle
paths were exported as a tabdelineated text file and read
into a spreadsheet, with the location of the particle (in
Cartesian x, y, and z coordinates) written at timesteps of
108. The distance each particle traveled between each time
step, L,, was calculated from,
Velocty
Inlet
SA
Figure 2: Computational fracture mesh.
A mesh of 688,000 unstructured tetrahedral cells was
applied to the full fracture model with the preprocessor
GAMBITTM. This level of mesh refinement was found to be
adequate for single phase flows; two phase flow models
were refined to approximately 1 million cells. The
unstructured single phase model surface mesh applied to the
entire fracture is shown in Figure 2, with a magnified insert
of a crosssectional profile to better illustrate the level of
refinement. The inlet and outlet boundary conditions are
shown as well. The surrounding fracture surfaces were
treated as nonslip walls.
A series of numerical simulations were performed using the
computational fluid dynamics solver, FLUENTTM.
Momentum was solved using a 2nd order upwind scheme,
pressure was solved using the pressure staggering option
scheme, pressurevelocity coupling was performed using
the embedded PISO algorithm, and the volume fraction
within each cell was determined with the geometric
reconstruction method. The solution gradients used by these
solvers were obtained by averaging node values along each
cell face to obtain solutions at the cell centers. The
geometries used for these analyses were threedimensional
models obtained from CT scans of a fracture in Berea
sandstone by Karpyn et al. (2007). Singlephase studies to
determine the fracture's 0, effective aperture bef, and
effective permeability kfe were conducted. Twophase
studies of the fracture were performed to determine the kr of
air and water in the fracture.
SinglePhase Flow Results
From the single phase simulations performed the 0, b, and
ke have been calculated. These values have been observed
to be similar to those expected for flow through a single
fracture. The methods used to calculate these properties and
the recorded values are discussed here.
As was previously mentioned, 0 is a quantity that describes
z1 )2
where the subscripts refer to the time steps, i being the
current time and i1 the previous time. L, was summed over
the entire distance the particles traveled and these particle
lengths were averaged for all the particles released; 200
particles were evaluated over the injection face. The 0 was
calculated with,
EL
L
where the fracture length L = 9.192 cm. The particles were
found to travel a distance 44% greater than the fracture
length, or 13.24 cm. Several representative particle paths
are shown in Figure 3. Note that all the particles traveled
through a narrow region near the bottom of this image,
regardless of their release location. This 0 value was
confirmed as being a reasonable value through discussions
with petroleum engineers familiar with this type of analysis
(Kadambi 2006).
Figure 3: Representative particle pathways through
fracture.
Darcy's Law for flow through porous media is often used as
an approximation of flow through fracture in reservoir scale
simulators (Wu et al. 2004). Darcy's Law predicts the
volumetric flow rate, Q, of a fluid with viscosity p through
a porous medium with crosssectional area A and length L
L, = (x, x,1)2 + y, y,)2 + (Z,
Paper No
when a pressure gradient, AP is imposed. The k can be
calculated by rearranging Equation 3.
APA
Q = k (3)
/L
This was accomplished with the numerical model by
changing the constant rate velocity inlet shown in Figure 2
to a constant pressure inlet, imposing pressure gradients of
0.1, 1, 10, and 100 Pa across the fracture, and measuring the
resultant mass flow rate of air and water through the
fracture. The properties of the working fluids for these
calculations are shown in Table 1. The mass flow rate was
used from the FLUENTTM simulations and converted to
volumetric flow rate with the densities (p) in Table 1. A plot
of the recorded flow rates for the case of water as the
working fluid are shown in Figure 4. From these results and
Equation 3, the effective permeability kef of the fracture was
calculated. kef refers to the fact that the fracture is not a
porous medium, but the narrow b resists flow in such a
manner that this relationship can be used with some
accuracy. The similarities of k and kef are shown visually in
Figure 5. The kef was calculated as 3.41(1010) m2 = 345
Darcies. This is order of magnitudes greater than the k of
the surrounding sandstone (Karypn et al. 2007) and matches
well to similar experimental and numerically measured
fracture permeabilities.
Table 1. Fluid properties for simulations.
p (kg/m3) (kg/m.s)
Air 1.225 1.8(105)
Water 998.2 1.003(10 3)
CO2 534 4(10 )
Brine 1052 103
10 07
c100 48
E
S10,
0
LL 10
10 11
0.1 1 10 1oc
Pressure Drop (Pa)
Figure 4: Volumetric flow rates of water through full
fracture model.
Sp
Figure 5: Permeability (k) through porous medium (left)
and effective permeability (keg) through a fracture (right).
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Using the same AP and Q results the bfe of the fracture was
calculated. The bfe assumes that the flow through the
fracture can be described by the cubiclaw for flow through
parallel plates. The cubic law is derived from the
NavierStokes equations with the assumptions of a narrow b
perpendicular to the flow through a wide and long set of
parallel plates. This relation has been used in discrete
fracture models and is given as
b3W AP
12= L
12/1 L
The bff of the fracture was calculated as 300.4 umr from
Equation (4). The similarity of this value to the smallest
aperture in the fracture model (240 gim) indicates that the
flow resistance within the fracture is dominated by the
smallest apertures, a fact that has been previously reported
in twodimensional fracture studies (Crandall et al. 2009)
and experimental studies of flow through fractures (Konzuk
and Kueper, 2004).
A representative image of pressure contours through the
fracture, with water as the working fluid and with an
imposed pressure gradient of 15 Pa is shown in Figure 6. As
can be seen the pressure drop across the fracture is
nonuniform, with constricted regions reducing the pressure
more rapidly. This is most apparent in the 12 Pa range of
the fracture flow shown in Figure 6.
PRM5.(Po
15.0
LTh 1.,
3.8
0.0
Figure 6: Pressure contours of water through fracture.
TwoPhase Flow Results
Using a refined fracture mesh simulations were performed
using the VolumeofFluid (VOF) model. Preliminary
simulations were performed using air and water as the
immiscible fluids, fluidfluid properties are listed in Table 2.
To approximate fluid behavior during sequestration brine
and CO2 properties were modeled as well, with the values
reported by Pashin and McIntyre (2003) for the Black
Warrior Basin used to approximate the downhole
conditions. Specifically these values represent the
approximate values of CO2 and brine with a salinity of
0.085 (85000 ppm NaC1) at a pressure of 11.7 MPa and
temperature of 69.6 C. Under these conditions a y of
35 "n/m was modeled, based upon the reported CO2brine y
at reservoir conditions by Chalbaud et al. (2008).
Table 2. Fluidfluid properties for twophase simulations.
Fluids Vn Pout Y 0
(cm/s) (kPa) (m/m) ()
(a) CO2 Brine 0.1 11700 35 90.0
(b) Air H20 0.1 1.01 72 80.0
(c) Air H20 1.0 1.01 72 90.0
Paper No
Three images of fluid invasion into the fracture, initially
filled with either water or brine are shown in Figure 7. As can
be seen the saturation varies significantly with the change in
the fluidfluid properties. The brineCO2 system has a much
lower y, thus the capillary resistances are much lower and the
invading CO2 easily moves into small regions of the fracture.
The saturations of air shown in Figure 7(b) and 7(c) show
that there is a dependency of the air saturation on the
direction of flow.
The relative permeabilities k, of the fluids moving through
this fracture were calculated as the permeability of each
individual fluid within four equal sized segments normalized
by the absolute permeability in that segment, e.g. k, I /.
The individual fluid permeabilities were calculated with,
, = Q,'
V_
where the subscripts indicate individual fluid properties. The
pressure gradient and flow of each fluid was measured across
the four individual segments of the fracture. Only when the
fluid was present at both edges of a segment could the k, be
determined. This segmentation of the fracture increased the
number of k, values that could be calculated, plotted as a
function of the water saturation (within each segment) in
Figure 8. Several important properties of these curves are
listed here. The irreducible water saturation in the first three
segments was approximately 5%, 8% and 11.5%, and thus a
wide range of air k, values were recorded at these values. The
k, of air reduced sharply from 1 to nearly 105 as the water
saturation increased from 0% to 30%. At water saturations of
greater than 30% the invading air did not form a continuous
path between segment faces and no k, values were recorded.
Conversely, the water k, was observed to reduce dramatically
at very low water saturation values, but remained fairly
constant from saturations of approximately 15% to 95%. A
sharp increase in the water k, was observed as the water
saturation in the segments approached 100%.
Figure 7: Twophase flow saturations of (a) CO2 (blue)
and brine (b) air (blue) and water injected in Side 1 and (c)
air (blue) and water injected in Side 2
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The wide range of values obtained from this unsteady
simulation of air into the waterfilled fracture is shown by
the all the 'calculated' points in Figure 8. Due to the sporadic
motion of the invading air, Q varied dramatically for nearly
identical saturation values. In a similar fashion, the pressure
within the invading air changed substantially as regions of
the fracture were invaded. These k, values may be best
described in terms of the probability of values at a distinct
saturation, or as the range of expected values. Work is
ongoing to explore these ideas, as well as to identify the
fluid properties that most significantly change the observed
flow saturations. The averaged values shown in Figure 8
were averaged over ranges of 2% and 5% for the air and
water, respectively.
0.01
2 0.001
0.0001
0.00001
0.000001 
0%
A
A
L^A
h^ ^o
AwSragd Wow kr
MjsurokgAir J
olAvo goed~torkr
10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
Water Saturation
Figure 8: Relative permeability (kr) curves of air and
water within the fracture.
Conclusions
We have performed simulations of single phase and
multiple phase flows in a reconstructed CTscanned fracture.
From the single phase simulations the tortuosity, effective
permeability and effective aperture were all determined.
Multiphase simulations of CO2 invading an initially brine
filled fracture revealed that, with low interfacial tension
values, the CO2 tends to fill the majority of the fracture,
including the small aperture locations. Similar simulations
of air invading a water filled fracture, with a higher
interfacial tension value, showed more fingering and less air
saturation. Relative permeability curves of air invading the
water filled fracture revealed a set of curves dramatically
different from the typical curves for flow in porous media,
due primarily to fluidfluid interactions along the length of
the fracture.
References
Amadei, B. and T. Illangasekare. 1994 A mathematical
model for flow and solute transport in nonhomogeneous
rock fractures. Int. J. Rock Mech. Min Geomech. Abstr.
36:719731.
Berkowitz, B. 2002. Characterizing flow and transport in
fractured geological media: A review. Adv. Water Res,
25:861884.
Brown, S.R., A. Caprihan, and R. Hardy. 1998.
Experimental observation of fluid flow channels in a single
fracture. J. Geophys. Res., 103:51255132.
Brown, S.R. 1989. Transport of fluid and electric current
through a single fracture. J. Geophys. Res., 94:94299438.
Paper No
Chalbaud, C., M. Robin, J.M. Lombard, F. Martin, P
Egermann, and H. Bertin, 2008. Interfacial tension
measurements and wettability evaluation for geological CO2
storage, Adv. Water Resour., 32(1):98109.
Crandall, D., G Ahmadi, and D.H. Smith. 2010.
Computational modeling of fluid flow through a fracture in
permeable Rock. Transport Porous Media, In Press.
Crandall, D., G Ahmadi, and D.H. Smith. 2009. Modeling
of immiscible, twophase flows in a natural rock fracture.
Proceedings of the ASME 2009 Fluids Engineering Division
Summer Meeting, FEDSM78138.
Crandall, D., G Ahmadi, D. Leonard, M. Ferer, and D.H.
Smith. 2008. A new sterolithography experimental porous
flow device. Rev. Sci. Instrum. 79(4):044501.
Dougan, L.T., PS. Addison, and W.M.C. McKenzie. 2000.
Fractal analysis of fracture: A comparison of dimension
estimates. Mech. Reas. Comm., 27:383392.
FLUENT. 2005. FLUENT 6.2 User's Guide, FLUENT Inc.
Lebanon, NH.
Glass, R.J., H. Rajaram, M.J. Nicholl, and R.L. Detwiler.
2001. The interaction of two fluid phases in fractured media.
Cur. Opin. Col. Int. Sci., 6:223235.
Glass, R.J., M.J. Nicholl, and L. Yarrington. 1998. A
modified invasion percolation model for lowcapillary
number immiscible displacements in horizontal rough
walled fractures: Influence of local inplane curvature. Water
Resource. Res., 34:32153234.
Hanano, M. 2004. Contribution of fractures to formation
and production of geothermal resources. Renew. Sustain.
Energy Rev., 8:223236.
Hughes, R.G and M.J. Blunt. 2001. Network modeling of
multiphase flow in fractures. Adv. Water Res., 24:409421.
Kadambi, J.R. 2006 Personal communication.
Karpyn, Z.T., A.S. Grader, and PM. Halleck. 2007.
Visualization of fluid occupancy in a rough fracture using
microtomography, J. Colloid Interface Sci.,
307(1):181187.
Konzuk, J.S. and B.H. Kueper. 2004. Evaluation of cubic
law based models describing singlephase flow through a
roughwalled fracture. Water Res. Reas., 40:W02402.
Klusman, R.W. 2003. Rate measurements and detection of
as microseepage to the atmosphere from an enhanced oil
recovery/sequestration project, Rangely, Colorado, USA.
App. Geochem., 18:18251838.
Lespinasse, M. and J. Sausse. 2000. Quantification of fluid
flow: hydromechanical behavior of different natural rough
fractures. J. Geochem. Explor., 69:483486.
National Energy Technology Laboratory, United States
Department of Energy. 2008. Carbon Sequestration Atlas of
the United States and Canada, www.netl.doe.gov, 2008.
Selley, R.C. 1998. Elements of Petroleum Geology, 2nd Ed.
Academic Press. Toronto, Canada.
McKoy, M.L. and Sams, W.N. 1997. Tight gas simulation:
Modeling discrete irregular stratabound fracture network
flow including dynamic recharge from the matrix. Number
P17, Presented at US Dept. of Energy's Natural Gas
Conference, Houston TX.
Pashin, J.C. and M.R. McIntyre. 2003. Temperaturepressure
conditions of the Black Warrior basin: implications for
carbon sequestration and enhanced coalbed methane
recovery, Int. J. Coal Geol. 54(3):167183.
Pruess, K. and Y.W. Tsang. 1990. On twophase relative
permeability and capillary pressure of roughwalled
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
fractures. Water Res. Reas., 26:19151926.
PyrakNolte, L.J., D.D Nolte, L.R. Myer, and N.GW. Cook.
1990. Fluid flow through single fractures. In Rock Joints, ed.
Barton and Stephansson, Balkema, Rotterdam, 405412.
RangelGerman, E., S. Akin, and L. Castanier. 1999.
Multiphaseflow properties of fractured porous media. SPE
54591, Presented at 1999 SPE Western Regional Meeting,
Anchorage, Alaska, May 2628.
Renshaw, C.E., J.S. Dadakis, and S.R. Brown. 2000.
Measuring fracture apertures: A comparison of methods.
Geophys. Res. Lett., 27:289292.
Roux, S., F. Plouraboud, and J.P Hulin. 1998. Tracer
dispersion in rough open cracks. Trans. Porous Media,
32:97116.
Tsang, Y.W 1984. The effect of tortuosity on fluid flow
through a single fracture. Water Resour. Res.,
4W0827:12091215.
Wu, Y.S., L. Pan, and K. Preuss. 2004 A physically based
approach for modeling multiphase fracturematrix
interaction in fractured porous media. Adv. Water Res.,
27:875887.
Zimmerman, R.W., D.W. Chen, and N.G.W. Cook. 1992.
The effect of contact area on the permeability of fractures. J.
Hydrology, 139:7996.
