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 low-permeability 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 micro-computerized tomography (CT) at a fine volume-pixel (voxel) resolution by Karpyn et al. (2007). This CT scanned fracture was converted into a numerical mesh for two-phase flow calculations using the finite-volume solver FLUENT® and the volume-of-fluid method. Simulations of single and two-phase flows in the reconstructed fracture geometry were performed. Initial two-phase 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 large-aperture 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 reservoir-scale discrete fracture models for predictions of multi-component 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 less-viscous CO2 fills the rough fracture geometry.

General Note:

The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows

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 low-permeability 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 micro-computerized tomography (CT) at a fine
volume-pixel (voxel) resolution by Karpyn et al. (2007). This CT scanned fracture was converted into a numerical mesh for
two-phase flow calculations using the finite-volume solver FLUE N T and the volume-of-fluid method. Simulations of single
and two-phase flows in the reconstructed fracture geometry were performed. Initial two-phase 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 large-aperture
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 reservoir-scale discrete fracture models for predictions of multi-component
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 less-viscous 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
long-narrow geometry of fractures has prompted the use of
"cubic-law" relations to describe flows with discrete
fracture models. The cubic law is derived from the solution
of the Navier-Stokes 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 core-scale (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 zero-aperture (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' (Pyrak-Nolte et al. 1990,
Brown et al. 1998), or regions of higher than average flow.
Two-phase 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 two-phase 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
(Rangel-German 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 two-phase
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. Pore-level

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 (Pyrak-Nolte 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 cubic-law 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 parallel-plates, but no unifying relationship
that accurately describes the flow of fluids within fractured
porous media has gained widespread acceptance.
We present full Navier-Stokes simulations of single and
multi-phase flows through a reconstructed CT scanned
fracture. We show agreement with previously reported 0
and k values and report a novel k, relationship for air-water
flows.

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
(volume-pixel) data set was obtained from these scans.
The 240 [tm voxel data was 'cleaned' using an in-house
code, which removed voids within the rock that were
unconnected to the rest of the fracture. The cleaned
three-dimensional 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
as-scanned 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.

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 single-phase flow of air through the fracture
and the steady state solution was obtained. Sub-micron
particles were released from the constant velocity inlet and
the particles were tracked using the step-by-step particle
tracker in FLUENTTM (2005). Brownian motion and
gravitational forces were not included, so that the particles
only followed the fluid flow. These step-by-step particle
paths were exported as a tab-delineated text file and read
into a spreadsheet, with the location of the particle (in
Cartesian x, y, and z coordinates) written at time-steps of
10-8. 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 pre-processor
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 cross-sectional 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 non-slip 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, pressure-velocity 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 three-dimensional
models obtained from CT scans of a fracture in Berea
sandstone by Karpyn et al. (2007). Single-phase studies to
determine the fracture's 0, effective aperture bef, and
effective permeability kfe were conducted. Two-phase
studies of the fracture were performed to determine the kr of
air and water in the fracture.

Single-Phase 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

z-1 )2

where the subscripts refer to the time steps, i being the
current time and i-1 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 cross-sectional 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.

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 cubic-law for flow through
parallel plates. The cubic law is derived from the
Navier-Stokes 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 two-dimensional 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
non-uniform, 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.

Two-Phase Flow Results

Using a refined fracture mesh simulations were performed
using the Volume-of-Fluid (VOF) model. Preliminary
simulations were performed using air and water as the
immiscible fluids, fluid-fluid 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 down-hole
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 CO2-brine y
at reservoir conditions by Chalbaud et al. (2008).

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 fluid-fluid properties. The brine-CO2 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 10-5 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: Two-phase 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

The wide range of values obtained from this unsteady
simulation of air into the water-filled 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
on-going 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 CT-scanned fracture.
From the single phase simulations the tortuosity, effective
permeability and effective aperture were all determined.
Multi-phase 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 fluid-fluid interactions along the length of
the fracture.