WATER RESOURCES RESEARCH, VOL. 47, W05534, doi:10.1029/2010WR009707, 2011
A study of conceptual model uncertainty in large‐scale CO2
storage simulation
Shuiquan Li,1 Ye Zhang,1 and Xu Zhang2
Received 30 June 2010; revised 4 February 2011; accepted 1 March 2011; published 25 May 2011.
[1] In this study, multiscale permeability upscaling is combined with a sensitivity study
of model boundary condition to identify an optimal heterogeneity resolution in developing
a reservoir model to represent a deep saline aquifer in CO2 storage simulation. A three‐
dimensional, fully heterogeneous reservoir model is built for a deep saline aquifer in western
Wyoming, where each grid cell is identified by multiple material tags. On the basis of
these tags, permeability upscaling is conducted to create three increasingly simpler site
models, a facies model, a layered model, and a formation model. Accuracy of upscaling
is evaluated first, before CO2 simulation is conducted in all models. Since at the injection
site, uncertainty exists in the nature of the reservoir compartment, end‐member boundary
conditions are evaluated, whereby brine production is introduced to control formation fluid
pressure. The effect of conceptual model uncertainty on model prediction is then assessed
for each boundary condition. Results suggest that for the spatial and temporal scales
considered, without brine production, optimal complexity of the upscaled model depends
on the prediction metric of interest; the facies model is the most accurate for capturing
plume shape, fluid pressure, and CO2 mass profiles, while the formation model is adequate
for pressure prediction. The layered model is not accurate for predicting most of the
performance metrics. Moreover, boundary condition impacts fluid pressure and the amount
of CO2 that can be injected. For the boundary conditions tested, brine production can
modulate fluid pressure, affect the direction of mobile gas flow, and influence the accuracy
of the upscaled models. In particular, the importance of detailed geologic resolution is
weakened when viscous force is strengthened in relation to gravity force. When brine
production is active, variability of the predictions by the upscaled models becomes smaller
and the predictions are more accurate, suggesting a subtle but important interplay between
heterogeneity resolution, fluid driving forces, and model predictions.
Citation: Li, S., Y. Zhang, and X. Zhang (2011), A study of conceptual model uncertainty in large‐scale CO2 storage
simulation, Water Resour. Res., 47, W05534, doi:10.1029/2010WR009707.
1. Introduction both onshore and offshore settings [Orr, 2009], particularly
for commercial‐scale operations.
[2] Carbon dioxide (CO2) is believed to be a driving force [3] In deep aquifers, CO exists in a supercritical “gas”
behind recent observations of global climate change. To 2phase that is less dense than formation brine. Under the
reduce the amount of CO2 entering the atmosphere, a variety imposed pressure gradient from injection, it will spread out
of actions are proposed, including CO2 capture from indus- laterally as well as rise and flow toward the caprock under
trial sources and subsequent storage into deep, permeable buoyancy, dissolving into brine as it migrates. Along the
geologic formations (carbon capture and storage, or CCS) trailing plume, a portion of the CO can become trapped due
[Intergovernmental Panel on Climate Change (IPCC), 2005]. 2to the imbibition of brine and associated CO relative per-
Candidates for such formations include unminable coal 2meability hysteresis [Flett et al., 2004]. Thus during and after
seams, depleted oil and gas reservoirs, and deep saline injection, dissolved, trapped, and mobile CO can all exist
aquifers. The last category, ubiquitous in the world’s sedi- 2in a storage formation. Over longer times, dissolved CO
mentary basins and easily accessible from CO2 point sources,
2
can react with the formation rock matrix, creating mineral
constitutes up to 90% of the total geostorage volume. Deep precipitates, or dissolution, or both, though fluid‐rock reac-
saline aquifers are thus considered ideal for CO2 storage in tions are considered slow in quartz‐dominated sandstones
[IPCC, 2005]. Currently, a variety of ongoing and planned
saline aquifer storage projects exist [Michael et al., 2009a,
1Department of Geology and Geophysics, University of Wyoming, 2009b], including both pilot‐scale and larger operations. At
Laramie, Wyoming, USA.
2 these sites, while CO2 storage is demonstrated to be tech-Schlumberger Information Solutions, Schlumberger, Houston,
Texas, USA. nically feasible, several issues must be considered before
commercial‐scale injection can become a reality.
Copyright 2011 by the American Geophysical Union. [4] First, commercial‐scale injection requires a storage
0043‐1397/11/2010WR009707 capacity that is orders of magnitude larger than those cur-
W05534 1 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
rently assessed by the existing operations [Michael et al., identify which one might be optimal for CO2modeling. Here,
2010]. At such scales, to quantify CO2 flow and storage the word “optimal” is used rather than “accurate” since it is
and to predict its footprint following injection and closure, widely acknowledged that a better resolution of subsurface
deep saline aquifers need to be characterized at the regional heterogeneity will lead to more accurate predictions of fluid
scale and, accordingly, large‐scale geologic site models must flow and transport. For example, in the work of Zhang et al.
be built. However, unlike hydrocarbon reservoirs which can [2006], we demonstrated that a facies model of higher geo-
be characterized with site laboratory or field data, deep saline logic resolution is more accurate than a formation model in
aquifers often suffer from extreme data scarcity. Using predicting groundwater flow. In thework of Zhang andGable
standard geophysical and borehole techniques to characterize [2008], we further demonstrated that although some aspects
them will make the storage project prohibitively expensive. of transport modeling (e.g., plume centroid and covariance)
For example, the cost of a seismic survey increases greatly can be captured by the facies model, other aspects (e.g.,
with lateral extent of the survey, depth of investigation, and plume breakthrough and tailing) require detailed heteroge-
resolution desired to map the storage system. Borehole neity. In CO2 modeling, however, particularly in assessing
sampling, which can provide downhole petrophysical prop- the problem of large‐scale storage in regional aquifers,
erties of the aquifer, e.g., porosity () and permeability (k), detailed description of heterogeneity in three dimensions
needs to be minimized to reduce cost and leakage risk [Carey down to the smallest resolvable continuum scale (e.g., core
et al., 2007]. Owing to these limitations, most saline aquifers measurements) is likely unobtainable. Thus a relevant ques-
are not extensively characterized, and accordingly, CO2 tion is not whether finer heterogeneity resolution should lead
modeling carried out as part of a site assessment study to more accurate predictions (we believe that the answer to
employs a variety of assumptions to simplify the construct of this question is yes) but whether we can identify a level of
the site model. In particular, depending on the quality and heterogeneity resolution in a geological storage model that is
accessibility of site‐specific data, subsurface environment is sufficiently accurate in predicting a variety of CO2 perfor-
commonly parameterized assuming permeability homoge- mance metrics.
neity at various scales. [7] The issue of optimal heterogeneity resolution is being
[5] For example, to predict basin‐scale CO2 storage explored by petroleum researchers concerned with reser-
capacity and plume footprint in response to a hypothetical voir production in extremely data‐poor environments [e.g.,
injection scenario in the Illinois Basin, the target saline Castellini et al., 2003; Friedmann et al., 2003]. Using a set
aquifer is modeled at the scale of a geological formation of characterization data for a deepwater reservoir, Milliken
[Birkholzer and Zhou, 2009; Zhou et al., 2010]. At such et al. [2007] built two families of geologic models at dif-
scales, field or laboratory data are often limited, and homo- ferent complexities. The same oil production scenario was
geneous properties are assigned to the model units [e.g., Law simulated to determine if the simpler models (e.g., a forma-
and Bachu, 1996; McPherson and Cole, 2000; Maldal and tion model consisting of a series of stacked channels) could
Tappel, 2004; Pruess et al., 2004; Nordbotten et al., 2005; be used to predict the behaviors of the more complex
Hesse et al., 2006; Sasaki et al., 2008; Stauffer et al., 2009]. models (e.g., each channel in the formation incorporates
On the other hand, where more data are available, within‐ five facies, where each facies is characterized with a dif-
aquifer heterogeneity has been modeled as distinct facies ferent and k distribution). The goal is to identify an
[Johnson et al., 2001; Hollowway et al., 2004; Obi and optimally simple model for flow predictions, since the cost
Blunt, 2006; Primera et al., 2009]. These units can include of building the complex models is extremely high. However,
interconnected high permeability (high‐k) channels embed- though the simpler models were identified to provide aspects
ded in low permeability (low‐k) clay, or reversely, low‐k of predictions similar to those of the complex models, the
clay or shale lenses embedded in sand or, in the case of two families of models are not strictly comparable, i.e., they
Doughty and Pruess [2004], multiple facies with variable do not provide the same bulk flow predictions under the
mean and k. By explicitly modeling facies distribution, same production scenarios.
effect of facies‐controlled heterogeneity on CO2 flow and [8] Besides the uncertainty in developing a site conceptual
storage can be quantified. Furthermore, approaches other model for CO2 simulation, another issue that is relevant for
than facies modeling have also been adopted. For example, the success of large‐scale storage is injectivity [Stauffer et al.,
Kumar et al. [2004] divided an aquifer into many layers, with 2009; Michael et al., 2010]. As commercial‐scale injection
each layer assuming a different permeability. Geostatistics‐ proceeds, a large volume of brine can be displaced from the
based heterogeneity has been investigated, which describes injection site. In a deeply buried aquifer, without an outlet to
either random or semistructural distribution of permeability dispose of the brine, formation fluid pressure can build up,
within aquifers [Flett et al., 2007; Ide et al., 2007; Qi et al., creating a footprint that can extend to regional scales [Zhou
2007]. Simulation results further suggest that models incor- et al., 2008; Birkholzer et al., 2009]. Though bottomhole
porating high‐resolution heterogeneity are needed to accu- pressure (BHP) constraint can be set to shut down the injector,
rately assess the various storage schemes [Juanes et al., 2006; little CO2 can in turn be injected, severely limiting injectivity.
Qi et al., 2009]. This view is echoed by petroleum researchers Without this constraint, however, damaging fluid pressure
who stated that the capture of small‐scale reservoir geometry buildup can potentially fracture the formation and its cap-
down to the scales of bedding planes can improve model rocks, creating leakage pathways. Neither option is accept-
estimates of hydrocarbon recovery [Elfenbein et al., 2005]. able for commercial‐scale storage of CO2 in saline aquifers.
[6] In CO2 modeling, a variety of approaches have thus During injection, a set of measures must be taken to address
been used to capture heterogeneity in a storage aquifer, some pressure buildup and brine displacement while ensuring
resolving heterogeneity down to fine scales, while others injectivity. These measures should also take into account
use simple configurations assuming homogeneity. These site boundary characteristics, i.e., whether the reservoir is
approaches have rarely been compared among one another to compartmentalized or not. Without appropriate well tests, this
2 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
information also suffers uncertainty. To date, most proposed [10] In the remainder of the article, the upscaling method is
strategies to address pressure buildup are based on brine introduced first, followed by a description of the reservoir
production [Flett et al., 2008; Qi et al., 2009; Surdam et al., models. A sensitivity study is presented, including simulation
2009; Wolery et al., 2009]. However, they have not been technique, parameters used, and injection design. Results
evaluated within the context of alternative conceptual models illustrate the impact of both conceptual model uncertainty and
and uncertain site boundary condition. Will a brine produc- boundary condition on CO2 flow and storage predictions.
tion strategy developed assuming formation homogeneity Findings of the study are discussed and summarized before
still be effective if the underlying heterogeneity is revealed? directions for future research are indicated.
And, is the boundary condition and its uncertainty important
to flow prediction? These questions are difficult to assess 2. Method
without prior information on reservoir heterogeneity.
[9] In this study, multiscale permeability upscaling is 2.1. Permeability Upscaling
combined with a sensitivity study of model boundary con- [11] To assess the effect of heterogeneity resolution on
dition to identify an optimal model complexity in developing CO2 prediction, besides a FHRM, a set of simple (or upscaled)a reservoir model for CO2 storage simulation. Unlike past models can be created under different simplifying assump-
research where models of various complexities were devel- tions, e.g., a formation that is homogeneous, one exhibit-
oped at different sites, a suite of reservoir models, each of a ing multiple facies, or one consisting of multiple layers. To
different complexity, is developed for the same site. The sim- facilitate model comparison, an upscaling analysis is used so
pler models are created based on an underlying and (assumed) that the simple models become equivalent to the FHRM in
fully known heterogeneous reservoir model (FHRM), which flow prediction. In the following sections, both intrinsic and
serves as a reference for predictions. To create a consistent relative permeability upscaling are discussed.
framework for model comparison, permeability upscaling is [12] In intrinsic permeability upscaling, equivalent perme-
used to link the simple models to the FHRM. CO2 injection is ability (k*) can be computed based on results of single‐phase
then simulated in all models. By comparing a suite of predic- steady state flow modeling in the FHRM. k* is assigned to
tion metrics against that of the FHRM, our goal is to identify a each unit of the upscaled models to represent the effect of
simple model that is sufficiently accurate for predicting gas unresolved, subunit permeability heterogeneity on single‐
migration, CO2 storage, and pressure buildup, thus contribut- phase predictions. Since in developing a site model for CO
ing to the development of a cost‐effective strategy in devel- 2storage, the model units can be irregular, this upscaling is
oping site models for CO2 simulations. The methodology accomplished using a numerical technique described by
presented is general and can be applied to any storage site. It is Zhang et al. [2006], where a 2‐D system was first analyzed.
tested here on a reservoir model developed for a proposed Herein, the earlier formulation is extended to 3‐D, with a few
commercial‐scale injection site in western Wyoming. modifications:
2 D E 3 8 9
@F @F @F
6 @x @y @z 0 0 0 0 0 0 > hqxi1 >6 1 1 1 7 > >6 D E 7 >> >>6 0 0 0 @F @F @F 70 0 0 7 > >66 @x 1 @y 1 @z 1 D E 7
> qy 1 >
6 776 8 9 >
>> >>>
6 0 0 0 0 0 0
@F @F @F
x y 7 kxx > hqzi >6 @ @ @z
> > 1
6 D E
1 1 7
1 7 >> >>
>> >> > > >6 >6 @F @F @F
7
@x @y @z 0 0 0 0 0 0 77 >> kxy6 2 2 2 D E 7 >> >
>> >>> hqxi2 >>>
66
> > >
6 0 0 0 @F @F @F
7 > > > >
@x @y @z 0 0 06 2 2 2 D E 7
7 >> kxz >> >> qy 2 >>
66
7
77 >
> >> >> >>
6 0 0 0 0 0 0 @F @F @F6 @x 2 @y 2 @z 2 7
>> kyx >7 <> >=> >
> hq
> z
i >2
6 < >
>=
6 76 . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 > kyy6 D E 7 > >
¼ ðÞ . . .
> >> >>
ð1Þ
66 7 > > > >6 @F @F @F 0 0 0 0 0 0 77 >>> kyz >>> >>> hqxi >6 @x m @y m @z m
m >
66 D E 7
7 >> > >
>>
6 0 0 0 @F @F @F 0 0 0 7 > k
>> >zx > qy >>
66
@x m @y mm @z m 77 >> >> >> >>
66
D E 7 > > > >
0 0 0 0 0 0 @F @F @F 77 >> k> zy
> > hqzi >
@x m @y @z m > > m >
66 m6 7
7 > > >
6 0 1 0 1 0 0 0 0 0 77
: ; > >
kzz >> 0 >>
6 > >6 77 >> >>6
46
0 0 1 0 0 0 1 0 0 77 >> 0 >5 :>
>>
;0 1 0 0 0 1 0 1 0 0
3 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Table 1. Symbols and Units of This Studya [15] Third, the numerical technique implemented in
equation (1) can be considered as a global method within the
Parameter Unit Symbol
context of single‐phase upscaling (see reviews [e.g.,Wen and
Porosity fraction Gomez‐Hernandez, 1996; Renard and de Marsily, 1997;
Permeability mdarcy k Sanchez‐Vila et al., 2006]). That is, k* of any unit of the
Natural log permeability mdarcy lnk
Log permeability mdarcy logk upscaled models is calculated using averages of potential
Horizontal permeability mdarcy k gradients and fluxes which are computed from global flowH
Vertical permeability mdarcy kV simulations that extend beyond the particular unit. It differs
Pressure Psia P from the local‐based methods in that upscaling is concerned
Temperature °F T with a model region (e.g., a facies unit) instead of seeking
Phase saturation fraction Si
Concentration of CO in brine ppm C coarse‐grid parameters from fine‐grid k distributions. The2 i
Residual water saturation fraction Sresw later is used in coarse‐graining research for the specific goal
Critical gas saturation fraction Scg of achieving higher computation efficiency in fluid flow
Residual gas saturation fraction Sresg modeling. Though an important topic on its own (e.g., the
Water endpoint relative permeability fraction kendrw
Gas endpoint relative permeability fraction kend global method is invariably more expensive than the localrg
Injection rate Mscf/d r methods, whether the upscaling concerns model regions ori
Production rate Mscf/d rp coarsened grids), this topic is not pursued here. The FHRM
Pore volume cf Vp and the upscaled models share the same grid and are
Component mass lb mol or Mt ‐ considered alternative conceptual models; thus numerical
Phase Density lb/ft3 ri
Phase Viscosity cP m discretization error does not come into play in model com-i
Fluid Potential Psi F parison. Nevertheless, it is useful to point out that global
Fluid potential gradient Psi/ft rF methods can be more accurate than local methods, since flow
Aquifer influx rate stb/d QR channeling due to permeability correlation outside the up-
aMscf/d is 1000 standard cubic feet per day; lb mol is pound mole; Mt is scaling domain can be accounted for in calculating k*. This is
one million metric ton; stb/d is standard barrel per day. confirmed by observations that the accuracy of a local method
increases when the size of the upscaling block surrounding
where h i represents spatial averaging among grid cells that the grid cell increases [Wen et al., 2003]. Given advances
belong to the same upscaling unit; qx, qy, qz are components of in computing techniques, limitation of the global method
the single‐phase Darcy flux; F is fluid potential (water is the may be overcome in time and its strength better realized in
fluid chosen for upscaling, thus F = rgh, h is hydraulic head); the future. In this study, the technique is used to help
r and m are fluid density and dynamic viscosity, respectively; establish single‐phase flow equivalence among various
subscripts 1, 2, …, m denote flow experiments conducted model representations.
using the FHRM, driven by distinct boundary conditions [16] The implementation of equation (1), however, is not
(m is ideally a large number, subjecting the model to straightforward with reservoir models employing irregular
many flow stimulations); kxx, …, kzz are components of the grids (i.e., grid lines nonorthogonal to one another nor
upscaled k* (symbols and units are listed in Table 1). For aligned with the simulation axes). The FHRM is built with a
each unit of the upscaled models, equation (1) must be commercial package (Petrel, 2009, available at http://www.
assembled and solved once. More details on how equation (1) slb.com/services/software/geo/petrel.aspx), resulting in an
is assembled from multiple flow experiments can be found in irregular corner point grid. The flow simulation is done
the work of Zhang et al. [2006]. For 3‐D upscaling, a few using Eclipse 300 [Schlumberger, 2009], an integrated finite
points need to be elucidated. difference (FD) simulator that interfaces with this grid. CO2
[13] First, the total number of flow experiments m must be storage is modeled with the GASWAT module of Eclipse
≥3 to obtain a unique solution. As discussed by Zhang et al. 300. A single simulation platform is thus used in permeability
[2006], when the upscaling domain size is larger than the upscaling, its verification, and CO2 modeling. However,
underlying ln(k) correlation ranges, k* becomes insensitive to Eclipse does not export fluid potential gradient, only potential
the number of m used and the associated flow patterns. To at the centroid of the grid cells. Owing to the irregular grid,
obtain k*, fewer flow experiments can thus be theoretically the gradient cannot be computed directly with finite differ-
carried out (in that 2‐D study, both linear and nonlinear flow ence. A postprocessing program is thus written which con-
patterns were tested for m = 2 and m = 4, respectively). verts the corner point grid to a finite element (FE) grid. This
Whether this is still true in 3‐D will be tested here as well. program is verified by conducting permeability upscaling for
[14] Second, the last three equations of equation (1) are a 3‐D synthetic aquifer with an orthogonal grid [Zhang et al.,
constraints added to ensure symmetry in k*. Depending on 2011]. Finite difference is used to compute the gradient,
the magnitude of the gradient terms in the coefficient matrix, which is compared to that computed with the FE program.
these equations may need to be rescaled to reduce numerical Results are consistent when multiple variances are tested. In
truncation error and to ensure a stable inversion. However, this study, each grid cell contains eight corner points; thus the
the enforcement of symmetry does not guarantee the creation FE formulation is based on eight‐node cells and linear shape
of a positive definite k*, which is required for modeling functions. Since Eclipse is a general‐purpose transient simu-
physically correct flow by the upscaled models. (In 2‐D lator, a special module is required to drive steady state flow
experiments with images exhibiting truncated heterogeneity, through the model, which is needed for k* upscaling and
1∼2% of the upscaled k* was found to be non‐positive‐ its verification. Detailed information on how this module
definite.) In this study, equivalent k* is tested for positive‐ is set up is described elsewhere [Zhang et al., 2011].
definiteness before being assigned to the upscaled models.
4 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
[17] Besides intrinsic permeability upscaling, effective twice as thick as the lower zone, is divided into two layers.
relative permeability functions (one for CO2 and one for The lower zone, despite its variable heterogeneity, is not
brine) can be obtained using results of two‐phase flow subdivided because (1) heterogeneity in this zone is irregular
modeling in the FHRM [Ertekin et al., 2001]. This approach and is not structured as layers and (2) we wish to examine the
is referred to as dynamic upscaling, and if the grid is simul- consequences of using simple layers while ignoring under-
taneously coarsened, the upscaled relative permeability lying variability. Finally, a formation model is created rep-
functions are referred to as pseudofunctions. Here, since all resenting the reservoir with one unit (Figure 1).
models share the same grid, only the effective relative per- [20] The creation of simple models using the above divi-
meability functions are of relevance. The effective functions sions is nonunique, as alternative approaches for label-
can be adjusted from those of the local functions (or “rock ing stratifications exist, e.g., based on percolation cluster or
curves”) assigned to the FHRM, i.e., those assigned to indi- connectivity analysis. The divisions adopted here reflect ap-
vidual grid cells that belong to different relative permeability proaches in CO2 modeling in data‐poor settings, as summa-
families, which is typically controlled by and measured for rized in section 1. The number of units in the upscaled models
individual facies in the reservoir. In upscaling for relative is used as a proxy for model complexity, since data require-
permeabilities, the goal is to ensure that the upscaled models ment and associated cost of building the model will likely
predict the same flow rate for each fluid phase, similar increase with the level of heterogeneity resolution. For
to single‐phase flow equivalence in intrinsic permeability example, the FHRM is considered the most complex and
upscaling. The effective functions are generally a tensor expensive to build, followed by the facies model (less
property that can vary with flow rate (e.g., viscous‐ versus expensive, since potentially petrophysical properties of the
capillary‐dominated, viscous‐ versus gravity‐dominated, facies can be borrowed from analog data [Milliken et al.,
etc.), flow direction, and boundary condition [Ataie‐Ashtiani 2007]), which is followed by the layered model (e.g.,
et al., 2002; Pickup et al., 2005; Braun et al., 2005]. Studies thickness‐based division), and the formation model. Note,
have shown that if the upscaling domain is large compared however, to develop insights into model complexity, this
to ln(k) correlation range, and if the rock curves are uniform study is conceptual, whereby the FHRM provides a basis to
in the reservoir, and if capillary pressure effect is negligible, create the simpler models. For a real reservoir, if data
the effective relative permeability of upscaled models may availability is such that a FHRM can be developed, simpler
be approximated by the local (and uniform) rock curves of models would not likely be created from its homogenization
the heterogeneous model [Sáez et al., 1989; Amaziane et al., because one would generally prefer to build as detailed a
1991; Das et al., 2004]. The first criterion is satisfied in this model as the underlying data allow. The simple models here
study (discussed later). To simplify the analysis, the same are hypothetical; only if no data exist to create detailed
relative permeability functions are assigned to the grid cells heterogeneity in the first place, these models would be used.
of the FHRM, disregarding potential differences among grid They are meant to represent alternative models created due
cells belonging to different facies (i.e., intrinsic permeabil- to lack of data.
ity groups). In CO2 modeling, capillarity is also ignored. [21] For each simple model, frequency distribution of the
Therefore the upscaled models employ the same relative underlying, within‐unit horizontal permeability (kx) is shown
permeability functions as those of the FHRM. (Figure 2). For the formation model, kx varies over 5 orders of
magnitude. Its histogram is multimodal with a ln(k) variance
2.2. Model Creation of 3.2, reflecting heterogeneity of the entire reservoir. A
significant volume fraction occurs near the median perme-
[18] The FHRM and its equivalent, upscaled models are ability of 21 mdarcy, corresponding to the thick and more
created based on a reservoir model developed for a deep uniformly distributed upper reservoir (Figure 1). For the
saline aquifer in western Wyoming (Nugget Sandstone), a facies model, since was separately modeled for each facies
proposed commercial‐scale CO2 storage formation. The unit using Sequential Gaussian Simulation and log(kx) was
reservoir model is heterogeneous with 134,064 grid cells mapped from [Li et al., 2011], permeability histograms of
(Figure 1). It was first upscaled using a local method from a the facies units are approximately unimodal with lognormal
3‐million‐cell geostatistical grid which has a laterally iso- distributions. These histograms correspond to the individual
tropic permeability with a vertical‐to‐lateral anisotropy ratio modes in the histogram of the formation model, as expected.
of 0.2 [Li et al., 2011]. Diagonal upscaling was used [Li et al., For the layered model, the upper two units have unimodal
2011]; thus local permeability of the FHRM is a diagonal distributions occurring near the median permeability, while
tensor. The geostatistical model was built with a hierarchical the bottom unit has a bimodel distribution spreading over
approach, integrating structural, facies, and petrophysical the full range of variability. In the lower reservoir, visual
data. Upon grid coarsening to create the FHRM, the (fine‐ inspection reveals a sheetlike high‐k structure in the southern
grid) facies types were conformably mapped onto the coarse region (Figure 1), with preferential continuity along the N‐S
model. Thus each cell of the FHRM is labeled with a facies direction. This continuity, embedded in low‐k deposit, is
tag. Using these tags, a four‐unit facies model is created for captured by the facies model (unit 4) but not by the layered
this study (Figure 1). model. The frequency distribution of the within‐unit vertical
[19] In building reservoir models with sparse data, facies permeability has similar characteristics as those described for
can be “picked” from wireline logs. Using these picks, kx and will thus not be presented.
layered models are created assuming lateral continuity away
from well control. At other times, layers may be created by
simple division of the reservoir thickness. In this study, a 2.3. CO2 Simulation
three‐unit layered model is created based on the observed [22] CO2 simulation is conducted with GASWAT (2009
variability in the FHRM (Figure 1). The upper zone, being version) [Schlumberger, 2009], a multiphase compositional
5 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Figure 1. Reservoir models evaluated in this study. From top to bottom are shown FHRM (horizontal per-
meability in mdarcy), four‐unit facies model (facies IDs are shown), three‐unit layered model (layer IDs are
shown), and a formation model. All models employ 15 times vertical exaggeration. Arrow points north.
6 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
simulator applicable to modeling CO2 flow in deep saline
aquifers. Two phases are considered: a CO2‐rich supercritical
phase (“gas”) and a H2O‐rich liquid phase. CO2 density is
computed using a cubic equation of state tuned to experi-
mental measurements, while liquid density is corrected for
total dissolved solids (TDS). Between the two phases, two
components (CO2, H2O) are modeled; CO2 andH2O can exist
in both phases. GASWAT first solves the pressure and molar
density of each component. Mole fractions of the components
in the phases are then computed through a flash process,
where mutual solubilities of CO2 and H2O are calculated to
match experimental data. Though another module (CO2ST-
ORE) has been cited, it is not used here due to the limitations
in its applicable pressure (up to 8700 psi) and temperature
(12∼100°C) ranges. If 8700 psi is set as the injector BHP,
CO2 cannot be injected into the deeply buried reservoir at a
sufficiently high rate. Maximum Nugget sandstone temper-
ature near the injection site is ∼110°C, which can be
accommodated by GASWAT. With this option, temperature
of the reservoir can vary with depth, which provides static
values to determine fluid properties (no energy balance
equation is solved). A temperature field is assigned to the
model by interpolating and extrapolating data from temper-
ature logs [Li et al., 2011].
[23] Brine TDS is set at a uniform value of 77,487 ppm,
reflecting a measured value near the injection site. To eval-
uate residual trapping, a hysteretic relative permeability
function is assigned to the gas phase (Figure 3). A non-
hysteretic function is assigned to the liquid phase, assuming
that the rock matrix is water wet. Both functions are based
on those measured in CO2/brine experiments on the Viking
Sandstone [Bennion and Bachu, 2005, 2006a, 2006b]. To
model flow reversal in a grid cell before maximum gas sat-
uration is reached, scanning curve interpolation between the
bounding relative permeability curves is used [Carlson,
1981]. During simulations, a single pressure is calculated,
assuming zero capillarity. Besides heterogeneity andmobility
effects, CO2 flow is typically dominated by viscous force
during injection and gravity override after injection ceases.
Figure 2. Histograms of within‐unit horizontal permeabil-
ity (kx, x is oriented east‐west): (a) formationmodel; (b) facies
model; (c) layered model. Permeability is shown in log scale;
tick marks within one log cycle (e.g., 10 to 100) represent 20,
40, 60, and 80mdarcy. For the facies and layered models, his-
tograms of the individual units are superimposed into the Figure 3. Relative permeability functions used in CO2
same diagram, thus the associated high frequencies. simulations. Krg is relative permeability for the gas phase;
Krw is relative permeability for the liquid phase.
7 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Figure 4. Location of the CO2 injector and brine producers. Model shown is the FHRM. The injector is
placed at Shute Creek gas plant, the proposed injection site.
Capillarity can locally enhance trapping, though it is not considered). When turned on (during the injection phase
modeled here so as to simplify the upscaling analysis. only), the producers are placed on a total fluid production
[24] To initialize the model and to compute fluid proper- rate constraint set equal to the CO2 injection rate adjusted
ties, an initial hydrostatic pressure distribution is assigned to by reservoir volume. The producers act to create pressure
the model which is fully saturated with the liquid phase. Rock drawdowns in the reservoir, simultaneously extracting brine
compressibility is assigned using a typical value for sand- and reducing pressure. This “voidage replacement” strategy
stones (9.81 × 10−7 1/psi assigned at a reference pressure of is often used in reservoir modeling to simulate secondary
5801 psi). A single injector is placed at the Shute Creek gas recovery [Kumar et al., 2005]. In this study, all rates were
plant (proposed injection site), fully perforating the reservoir adjusted during an initial set of experiments, whereby the
(Figure 4). An injection rate is set at 83,515 Mscf/d (1.70 Mt/ simulation progress was monitored for convergence issues.
year), for a duration of 600 years. The injection phase is The final set of parameters give consistent and stable solu-
followed by a 2500‐year monitoring period, for a total sim- tions for all models, e.g., during injection, both the CO2
ulation time of 3100 years. The injection time is chosen to be injector and brine producers (if turned on) are continuously
artificially large so as to maximize the filling of CO2 in the operating without violating the BHP constraint.
reservoir when a single injector is used. The injection rate is [27] Boundary conditions (BC) of the model can pro-
also a target rate; the actual rate is modified by the simulator foundly influence the evolution of reservoir pressure, CO2
to satisfy a set of well constraints (discussed below). The plume, and formation sweep. Typical for deep aquifers (at the
actual rate can vary with time. injection site, Nugget Sandstone lies at depth of ∼13,000 ft),
[25] To prevent geomechanical damage to the storage significant uncertainty exists in the nature of the reservoir
formation and overlying caprocks, an injector BHP constraint compartment. Though well log correlation near the injection
is set using 1.8 × hydrostatic pressure at the reservoir depth site suggests lateral continuity of the sandstone [Li et al.,
(10,260 psi). This threshold is selected based on results of 2011], extrapolation to regional scale is difficult due to lack
pressure leak‐off tests (LOT) from two wells located near of well control away from the site. Far west of Shute Creek,
Shute Creek. LOT can be used to determine in situ fracture in an area called Wyoming‐Utah Thrust Belt, the Nugget
gradient of a formation. In the Nugget Sandstone, LOT yield Sandstone is thrusted toward the surface. Here, the formation
an equivalent mud weight of about 1.8 g/cm3, which trans- was observed to be compartmentalized by bounding faults
lates to the above fracture gradient if fresh water density is [Lindquist, 1988]. In this study, end‐member boundary
used for formation fluid (1.0 g/cm3). This fracture gradient is conditions are tested to assess their impact on CO2 storage.
consistent with those estimated for many hydrocarbon fields Four BC are investigated, with and without brine production.
[Eaton, 1969]. BC 1 is all‐closed‐no‐brine, where all boundaries are no‐flow
[26] During simulations, the injection rate is continuously (closed) and brine production is not activated. BC 2 is all‐
adjusted by the simulator so that the maximum formation closed‐with‐brine, where all boundaries are closed, but brine
fluid pressure does not exceed the BHP constraint. However, production is turned on during CO2 injection. BC 3 is open‐
the injector can eventually be shut down if no suitably small no‐brine, where all boundaries are closed except the north
rate can be found to satisfy the constraint. To control pressure and south faces, where an open boundary is assigned. No
buildup, six brine producers are completed near the southern brine is produced. BC 4 is open‐with‐brine, where north
end of the model. The producers are downdip from the face is open and south face is closed, while brine produc-
injector, each fully perforating the formation. (An opposite tion is activated during CO2 injection. In all the simulations
placement in the northern end will enhance updip migration with brine production, the producers are shut down when
and reduce the overall formation sweep; this option is not CO2 injection ceases, but the conditions at the boundaries
8 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Table 2. Equivalent Permeability of the Upscaled Models Computed Using m = 3, 4, and 6a
Upscaled Permeability (mdarcy)
Stratigraphic Model Unit m kxx kxy kxz kyy kyz kzz
Facies unit 1 3 22.61 0.22 −0.06 25.61 −0.11 1.36
4 22.61 0.22 −0.05 25.62 −0.12 1.32
6 22.60 0.22 −0.04 25.62 −0.13 1.29
unit 2 3 17.31 0.12 0.06 19.41 0.13 0.98
4 17.31 0.12 0.06 19.40 0.16 0.94
6 17.31 0.12 0.07 19.40 0.15 0.90
unit 3 3 7.14 1.55 0.38 23.08 −0.80 0.41
4 7.34 1.83 0.20 22.24 −0.39 0.21
6 7.40 1.91 0.14 21.99 −0.25 0.14
unit 4 3 1014.70 137.12 −0.12 1082.16 −15.35 0.42
4 1020.37 135.69 0.45 1077.41 −16.28 0.36
6 1027.24 136.93 0.44 1065.18 −16.70 0.34
Layered unit 1 3 20.99 0.21 0.15 23.89 −0.20 1.54
4 20.99 0.21 0.15 23.89 −0.19 1.54
6 20.99 0.21 0.15 23.89 −0.19 1.54
unit 2 3 19.36 0.29 −0.07 21.52 0.12 0.98
4 19.36 0.27 −0.05 21.58 0.14 0.91
6 19.35 0.27 −0.04 21.58 0.13 0.86
unit 3 3 302.86 72.42 −0.54 334.19 −4.92 0.36
4 303.31 72.36 −0.67 332.66 −5.22 0.33
6 304.09 72.34 −0.917 329.90 −5.703 0.32
Formation unit 1 3 92.35 18.59 −0.02 96.89 −1.03 0.69
4 92.36 18.58 −0.01 96.93 −1.02 0.64
6 92.44 18.59 −0.08 96.69 −1.18 0.61
aSee text for details. Location of the units is shown in Figure 1. When m = 4, the boundary conditions are x flow, y flow, z flow, and injection x flow. Two
additional sets of experiments are conducted for m = 4, using x flow, y flow, z flow, and injection y flow and injection z flow, respectively. The upscaled
permeability is very close to what is obtained with injection x flow. Results of these tests are not shown.
remain unchanged. The open boundary is modeled with the along the x axis or injection x flow); (5) water is injected at
Fetkovich Aquifer Facility [Schlumberger, 2009], which the same location and at the same rate and all boundaries
allows the resident brine, and later the injected CO2, to flow are closed except the north and south faces (injection y
out of the model into an external aquifer where pressure flow); (6) water is injected at the same location and at the
response is assumed to be felt uniformly throughout. This in same rate and all boundaries are closed except the top and
effect represents a reservoir that is semi‐infinite in the bottom faces (injection z flow). Note that experiments 1–3
direction where an open boundary is prescribed. The oppo- create linear flood patterns. The specified heads are selected
site situation of the closed system reflects the existence of such that in each experiment, flow is driven toward the
sealing faults that act to compartmentalize the reservoir. positive axis. Experiments 4–6 create diverging flow patterns
[28] To evaluate multiple conceptual models, all models where the water injector is fully perforating the reservoir.
are simulated under identical conditions (e.g., boundary The upscaled permeability for each unit is calculated based
condition, injection rate, brine production, well constraints), on the first three experiments, then based on four
the only difference being that in the simple models, equiva- experiments (x flow, y flow, z flow, and injection x flow
lent permeability and average bulk porosity are assigned to or injection y flow or injection z flow; three sets of k* are
each model unit. Thus no calibration is attempted by varying computed), and finally, based on all six experiments.
the parameters of these models to fit the predictions of the [30] Results suggest that for all units, k* are diagonally
FHRM. The bulk porosity is estimated by a volume‐weighted dominant full tensors (Table 2). All k* have also passed the
average of the subunit cell porosities. positive‐definiteness test; no physically incorrect values
were found. The diagonal dominance reflects the nearly flat
3. Results stratification in the FHRM and the fact that heterogeneity
principal correlation axes are approximately parallel to the
3.1. Permeability Upscaling: Sensitivity Analysis simulation axes (Figure 1). Because of stratification, the lat-
[29] To analyze the intrinsic permeability of the upscaled eral components of k* (e.g., kxx, kyy) are also greater than the
models, sensitivity analysis on the number of flow experi- vertical components (e.g., kzz), as expected. Unit 4 of the
ments (m) is conducted to understand whether the computed facies model has much higher equivalent permeability than
k* is sensitive tom and the associated flow configuration. Up unit 3 of the same model (both lie in the lower reservoir),
to six global boundary conditions are tested (m = 6): (1) x flow attesting to the fact that the facies model has captured the
(specified heads along the west and east faces; no‐flow on all bimodal distribution in this region by separating the deposit
other faces); (2) y flow (specified heads along the north and into two units. In the layered model, the lower reservoir is
south faces; no‐flow on all other faces); (3) z flow (specified represented by its own unit 3, whose k* lies between those of
heads along the top and bottom faces; no‐flow on all other units 3 and 4 of the facies model. Further, k* of units 1 and 2
faces); (4) water is injected at a constant rate into the model of the facies model are similar to those of units 1 and 2 of the
center and all boundaries are closed except the east and west layered model, despite the different division scheme. These
faces (i.e., injection‐induced radial pattern dominated by flow four units lie in the thicker upper reservoir, where the level of
9 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Table 3. Aquifer Inflow Rates Computed By All Modelsa [34] Histogram of the prediction error in fluid potential
(DF = Ffw − Fref) is also shown for each model, one for
Permeability X Flow (W‐E) Y Flow (N‐S) Z Flow (Vertical) each boundary condition tested. Under the same boundary
Model QR Error(%) QR Error(%) QR Error(%) condition, the potential drop used to drive flow through the
FHRM 56161.25 ‐ 44248.63 ‐ 19413.37 ‐ models is the same. But, for the same model (Figure 5), the
Facies 49368.26 −12.87 36021.43 −20.50 19408.29 −0.03 potential drop varies with boundary condition. The histogram
Layered 57219.8 1.87 47622.66 7.35 19254.10 −0.82 ranges from unimodal to bimodal, with a similar distribution
Formation 63305.82 11.96 57323.73 25.75 19322.87 −0.47
when models are simulated under the same boundary con-
aQR is aquifer inflow rate. For each upscaled model, a percent error is dition. Flow direction thus controls the shapes of the error
computed based on the results of the FHRM. histograms rather than difference in conceptual models.
Under the same boundary condition,DF range is consistently
variability is low. The lateral components of k* of these units the smallest for the facies model. For the layered and for-
range from 17.3 to 25.6 mdarcy, encompassing the median mation models, the spread of DF is larger, although they are
value of the full distribution (21.0 mdarcy). In the formation similar in magnitude.
model, only a single k* is computed. [35] The above results suggest that accuracy in flow rate
[31] Table 2 also lists the components of k* computed prediction appears to be controlled by boundary condition
using each of the flow configurations. By comparing their (i.e., whether flow is parallel or perpendicular to stratifica-
values, changes in k* are observed to be small. For example, tion), while accuracy in fluid potential prediction appears to
for the layered model the largest change in k* components is be controlled by themodel (i.e., finer heterogeneity resolution
0.4% (betweenm = 3 and 4) and 1.3% (betweenm = 3 and 6); has better accuracy). The magnitude of flow rate prediction
for the facies model the largest change is 0.6% and 1.6%, error does not correspond to the magnitude of fluid potential
respectively; for the formation model, the largest change is prediction error, e.g., the facies model is not necessarily the
6.7% and 0.2%, respectively. This suggests that k* is mostly most accurate in flow rate prediction. This suggests that when
insensitive to the number of flow experiments used in up- multiple boundary conditions are considered, no model is
scaling and the associated flow patterns. The additional best in predicting both flow rate and fluid potential. The
nonlinear flow experiments do not significantly influence the magnitude of the errors, however, is considered acceptable
upscaled permeability compared to those obtained with linear compared to those encountered in previous studies. For
floods. The upscaled permeability obtained from linear flow example, 2‐D results of Zhang et al. [2006] gave rise to a
experiments (i.e., m = 3) are used in subsequent single‐phase maximum flow rate prediction error of −31% and a maximum
verification tests and CO2 storage modeling. head prediction MRE of 9%. In that study, errors were also
found to be sensitive to boundary condition and homogeni-
zation level. Given the overall compatibility of the results
3.2. Permeability Upscaling: Verification Tests with the expected error margins, the upscaled models are
[32] Equivalent permeability computed for the upscaled considered adequate to represent the FHRM as equivalent
models is verified by conducting single‐phase (water), models for single‐phase processes.
steady‐state, incompressible flow simulations in all models,
with the FHRM providing reference predictions. The same 3.3. CO2 Flow Modeling
boundary conditions used in upscaling (i.e., x flow, y flow, z
[36] Commercial‐scale CO injection is simulated in all
flow) are tested, following Zhang et al. [2006]. Three sets of 2
models until the reservoir is sufficiently filled, followed by
Constant‐Head Darcy tests are conducted, with each set
postinjection simulation to monitor plume migration. To
simulating all models driven by the same boundary condition.
understand pressure buildup in the reservoir, two boundary
Each upscaled model computes an aquifer inflow rate, i.e.,
conditions are tested to represent a reservoir that is com-
flow rate between an external aquifer of higher head and the
partmentalized and one that is open on the sides (in this case,
reservoir model (Table 3). A percent relative error is com-
brine will be displaced out of the reservoir into an external
puted based on the flow rate of the FHRM. This error ranges
aquifer). These two conditions define end‐member scenarios,
from less than 1% to 26%, with the z flow simulations gen-
given the uncertainty of the boundary condition at the injec-
erally the most accurate. (An outflow rate, between the res-
tion site. To also address injectivity, for half of the runs, six
ervoir model and an external aquifer of lower head, is also
brine producers will operate during CO2 injection. Finally, tocomputed. At steady‐state, its value is verified to be nearly
evaluate performance of the upscaled models against that of
identical to the inflow rate and is thus not shown.)
the FHRM, four prediction metrics are used.
[33] To assess the accuracy in predicting fluid potential, a
[37] 1. The first is total predicted gas‐in‐place (GIP) and
mean error (ME) is defined for each upscaled model:
X gas storage ratio (GSR) at the end of the simulation (Table 5).n
¼ 1 ðiÞ ðiÞ ð Þ GSR is mass fraction of the total dissolved and trapped gas inME Ffw Fref 2n ¼ GIP. It is less than 1.0 due to the existence of mobile gas.i 1
where F(i)fw represents fluid potential of the ith grid cell com- Table 4. MRE in Fluid Potential Computed for Each Upscaled
puted by an upscaled model, F(i)ref represents the same fluid Model
potential of the FHRM, and n is the total number of grid cells
in the model. A dimensionless mean relative error (MRE) is Permeability X Flow (W‐E) Y Flow (N‐S) Z Flow (Vertical)
Model MRE (%) MRE (%) MRE (%)
also computed by normalizing ME with absolute potential
drop across the model which drives flow (Table 4). TheMRE Facies 2.24 4.28 0.37
ranges from less than 1% to 11%, with the facies model Layered 4.18 10.57 0.64
consistently being the most accurate. Formation 4.16 9.84 0.94
10 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Figure 5. Histogram of the prediction error in fluid potential (DF = Ffw − Fref; in psi) for each upscaled
model, under (left) x flow, (middle) y flow, and (right) z flow conditions. (top) Facies model; (middle) lay-
ered model; (bottom) formation model.
[38] 2. The second is average reservoir fluid pressure and [40] 4. The fourth is total predicted gas (in mass unit) over
its evolution over time. It is lower than the injector BHP but time or gas profiles.
greater than the producer BHP (if activated). [41] In Table 5, results pertaining to the first performance
[39] 3. The third is gas saturation map and its evolution metric are presented for all boundary conditions. In all gas
over time. categories, facies and formation models are more accurate
11 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Table 5. Gas‐in‐Place in the Reservoir Predicted By All Models at the End of Simulationa
Gas‐in‐Place (lb mol) (×1010) Storage Ratio Total Injected Gas Injection Time Injection Rate
Model BC Dissolved Trapped Mobile Total (%) Pound Mole (×1010) Megaton (years) (Mt/yr)
FHRM 1 0.03 0.16 0.04 0.23 0.84 0.23 44.47 50 0.89
2 0.20 1.97 1.72 3.89 0.56 4.23 817.87 600 1.36
3 0.32 2.92 2.21 5.45 0.59 5.49 1061.49 600 1.77
4 0.24 2.14 1.62 4.00 0.60 5.10 986.09 600 1.64
Facies 1 0.03 0.17 0.03 0.23 0.87 0.23 44.47 50 0.89
2 0.20 2.09 1.70 3.99 0.57 4.12 796.60 600 1.33
3 0.32 2.81 2.31 5.44 0.58 5.48 1059.56 600 1.77
4 0.25 2.28 1.59 4.12 0.61 5.03 972.55 600 1.62
Layered 1 0.03 0.17 0.03 0.23 0.87 0.23 44.47 50 0.89
2 0.21 1.78 1.26 3.25 0.61 3.52 680.59 600 1.13
3 0.21 2.10 1.97 4.28 0.54 4.36 843.01 600 1.41
4 0.23 2.00 1.36 3.59 0.62 5.04 974.48 600 1.62
Formation 1 0.03 0.17 0.03 0.23 0.86 0.23 44.47 40 1.11
2 0.21 2.22 1.76 4.19 0.58 4.35 841.07 600 1.41
3 0.18 2.29 2.88 5.35 0.46 5.37 1038.29 600 1.73
4 0.23 2.49 2.38 5.10 0.53 5.13 991.89 600 1.65
aGIP is gas‐in‐place. Four boundary conditions are tested (see text for details). The total injected gas is listed, as well as the actual injection rate. For some
cases, CO2 has flowed out of the reservoir; thus total GIP < total injected gas. A CO2 storage ratio is defined as (Dissolved GIP + Trapped GIP)/Total GIP ×
100%. For BC 1, the injection time is defined by the time when the actual injection rate is larger than 1000 Mscf/day.
than the layered model. For all models, BC 1 results in much observed. This postinjection gravity segregation is signifi-
less CO2 being injected, but the highest GSR. In a sealed cant, owing to the fairly high vertical permeability of the
system, the injector is shut down when the BHP constraint is models (i.e., cell permeability of the FHRM is weakly
reached. The higher ratio may have been due to the fact that anisotropic while permeability of the upscaled models is
formation fluid pressure in BC 1 simulations is the highest derived from those of the FHRM). Beneath this mobile high‐
among all boundary conditions tested. Higher pressure leads saturation zone, reservoir cells are dominated by trapped
to more CO2 dissolution in brine as well as to higher gas CO2, often at or below the residual gas saturation (0.48). This
saturation at the end of injection. This higher gas saturation, fraction of the gas plume is immobile and will not migrate
in turn, leads to more residual trapping, following a scanning further.
curve that is closer to the bounding imbibition curve. Simu- [44] In sections 3.3.1–3.3.5, for each boundary condition
lations under BC 2, 3, and 4, however, lead to similar GSR, tested, time profiles of average reservoir pressure, saturation,
but the amount of total gas injected and total gas stored vary and gas predictions are presented. Gas saturation is visualized
by ∼27%. This suggests that, given the current well con- along a north‐south transect through the model, crossing the
straint, boundary conditions can impact the amount of gas injector at Shute Creek. The saturation shown is the total
that can be injected and subsequently stored in the reservoir. combined mobile and trapped CO2. Dissolved CO2 in for-
[42] For each simulation case, the actual injection rate mation brine is not shown. The dissolved plume has a shape
achieved is computed by dividing the total injected similar to the gas‐phase plume, but its (aqueous) concentra-
gas by the duration of injection (Table 5). It ranges from tion is more uniformly distributed in space.
0.89∼1.77 Mt/year, with BC 3 and 4 predicting higher rates. 3.3.1. BC 1: All‐Closed‐No‐Brine
Considering that these cases simulate a reservoir linked to an [45] Under this BC, fluid pressure rises quickly in all
infinitely large external aquifer, this result is expected. With models from the initial hydrostatic pressure (∼5700 psi), and
one injector, the actual operation rate will likely lie within the the injector is shut down when the reservoir pressure reaches
above range, between a completely sealed and a completely the BHP constraint after ∼50 years of simulation (Figure 6).
open system. However, the simulations use an artificially The gas‐phase plume is small at the end of 600 years
long injection time to create a single plume to fill the res- (reservoir is still monitored after the injector is shut down),
ervoir. In actuality, the injection time will be shorter, corre- filling a small part of the formation near the injector. During
sponding to the lifetime of power plants. If a single injector the actual injection time, the average rate is ∼0.89 Mt/year,
is used, the extent of the actual plume at the decommission too small to sustain a commercial operation. All models
time will be much smaller than those simulated here. Clearly, predict similar pressure as well as shape and size of the gas‐
realistic simulations will require more injectors operating and dissolved‐phase plumes (not shown). Under this extreme
over shorter times. Such scenarios will be presented in condition of low injectivity, gas migration is limited and all
section 4. The focus of this section is conceptual and does models are approximately of equal accuracy.
not reflect actual practices. 3.3.2. BC 2: All‐Closed‐With‐Brine
[43] Under boundary conditions 2, 3, and 4, which are of [46] With brine production, significantly more CO2 can
greater interest, the dominant flow of CO2 in the formation is be injected into the reservoir and formation pressure never
lateral and upward during injection, driven by the pressure exceeds the BHP constraint (Figure 7). During injection, the
gradient in the reservoir and buoyancy of the CO2. After average reservoir pressure drops gradually until reaching a
injection ceases, gravity override dominates. The gas‐phase steady, subhydrostatic value (∼3800 psi) at the end of injec-
plume continues to rise, spreading out beneath the formation tion (600 years; 219,000 days). This is due to brine produc-
top. At the end of simulation, a laterally extensive high‐ tion. The level of pressure drop is controlled by the water rate,
saturation zone, consisting mostly of mobile CO2, can be e.g., larger drop can result from higher production rate (when
12 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Figure 6. BC 1, all‐closed‐no‐brine. (top) Gas‐phase CO2 saturation along a north‐south transect after
600 years of simulations. Shown is the prediction by the FHRM. (bottom) Average reservoir pressure over
time predicted by all models.
the voidage constraint is removed). All models simulate a create pressure drawdown in the south which counteracts
large regional‐scale plume, with the facies model being the updip flow.Among all the upscaledmodels, the faciesmodel is
most accurate in reproducing the plume shape of the again the most accurate in both plume shape and pressure
FHRM, at both time scales (i.e., end‐of‐injection and end‐of‐ predictions, and the layered model is the least accurate. The
monitoring). It is also the most accurate in predicting the formation model predicts the high‐saturation mobile gas zone
pressure profile. The layered model appears to be the least reasonably well, though it underestimates the extent of the
accurate in predicting both pressure and plume shape, e.g., a trapped gas in the downdip reservoir. In this model, despite
tongue of gas plume is seen in the lower reservoir which is its reasonable pressure prediction, lack of heterogeneity results
absent in the FHRM. The predictions of the formation in less plume spreading and, as a result, less trapping. In BC 2,
model are not significantly different from those of the facies however, formation model is comparatively more accurate;
model. Despite a slightly larger deviation in predicting the brine production in that case enhances downdip sweep and thus
pressure profile, it is sufficiently accurate to predict the trapping, despite the lack of heterogeneity.
extent of the plume footprint at both time scales. Overall, [49] Gas plume predicted by the layered model is signifi-
predictions of the upscaled models are not significantly cantly different from those of the other models, with most of
different from those of the FHRM. the injected CO2 channeling through the lower unit to reach
3.3.3. BC 3: Open‐No‐Brine and escape through the open boundaries on the north and
[47] Under this BC, again, significant amounts of CO2 can south sides. This feature persists during repeated simulations
be injected into the reservoir and formation pressure never when results are inspected at earlier times; it is an authentic
exceeds the BHP constraint (Figure 8). However, compared outcome, which develops soon after the injection starts. The
to the pressure profiles of BC 2, pressure rises in all models lower unit of this model, which contains a bimodal perme-
during injection and then drops to the preinjection level later ability distribution (Figure 2), is homogenized by a single
on. This rise (with the exception of the layered model) is equivalent k*. Though single‐phase tests under linear floods
gradual and pressure buildup in the reservoir is small. This is do not indicate significant errors in flow rate and fluid
because, during injection, formation brine is continuously potential predictions (they are similar in magnitude to those
being displaced into the external aquifer through the open estimated for the facies and formation models), such an
boundaries. Pressure buildup near the injector is thus contin- equivalent k* appears inadequate in two‐phase modeling.
uously attenuated. Compared to the facies model that conforms to heterogeneity,
[48] With the exception of the layered model, all models the layered model averages permeabilities from grid cells
predict a large, regional‐scale plume, with significant belonging to facies with different mean k. Though this does
updip (mobile) gas migration occurring during monitoring. not cause problems in single‐phase prediction, it is suffi-
Accordingly, in these three models (FHRM, facies, forma- ciently inaccurate for CO2 modeling, which is additionally
tion), sweep efficiency is higher in the updip region. This affected by mobility, gravity, and flow rate (i.e., strength of
updip migration is more significant compared to that pre- the viscous driving force). The pressure profile predicted by
dicted by the samemodels under BC 2, where brine producers this model is also the least accurate. It rises higher during
13 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Figure 7. BC 2, all‐closed‐with‐brine. (top) Gas‐phase CO2 saturation along the same transect of
Figure 6. (left) End of injection (600 years or 219,000 days) and (right) end of monitoring. First row
is FHRM; second row is facies model; third row is layered model; fourth row is formation model.
(bottom) Average reservoir pressure over time predicted by all models.
14 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Figure 8. BC 3, open‐no‐brine. (top) Gas‐phase CO2 saturation. (left) End of injection and (right) end of
monitoring. First row is FHRM; second row is facies model; third row is layered model; fourth row is for-
mation model. (bottom) Average reservoir pressure over time predicted by all models.
15 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Figure 9. BC4, open‐with‐brine. (top) Gas‐phase CO2 saturation. (left) End of injection and (right) end of
monitoring. First row is FHRM; second row is facies model; third row is layered model; fourth row is for-
mation model. (bottom) Average reservoir pressure over time predicted by all models.
16 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
17 of 23
Figure 10. Gas‐in‐place over time predicted by all models under BC 2, 3, and 4. (left) BC 2; (middle) BC 3; (right) BC 4.
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
injection compared to those predicted by the other models. At [53] The mobile gas grows linearly during injection, since
∼160,000 days (438 years), it reaches the BHP constraint and the total amount of the injected gas is increasing. After the
the injection rate is adjusted downward to maintain a lower injection ceases, the amount of the mobile gas drops gradu-
pressure until the end of injection (the simulator is able to ally, corresponding to the simultaneous increase of both the
find lower enough rates so that the injector is not shut down). trapped and dissolved gas.
This is not the behavior of the FHRM. [54] Similar to what is observed with the plume shape, for
3.3.4. BC 4: Open‐With‐Brine all the gas categories, throughout the simulation time, when
[50] Under this BC, the main difference from BC 3 is that BC 2 and 4 are used (brine producers are active), variability
brine producers are active during CO2 injection. Several of the predictions by the upscaled models is smaller and the
important differences are noted from the results of BC 3. First, predictions themselves are more accurate. Without brine
average pressure profiles predicted by all models are nearly production, as in BC 3, the facies model is again significantly
equally accurate. Owing to the producer rate constraint, all more accurate than the layered and formation models.
pressures are maintained at approximately the initial value.
Clearly, brine production exerts an important control on the 4. Discussions
evolution of formation pressure. Second, during the injec-
tion phase, in all models, the dominantly updip migration [55] Feasibility studies on CO2 storage in deep saline
(Figure 8) is replaced by downdip migration (Figure 9). aquifers use reservoir simulation with a geologic site model.
Accordingly, all models predict more sweep efficiency and However, multiple conceptual models can be developed
trapping in the downdip region. Clearly, brine production depending on data support. Since increased costs can be
exerts an important control on the direction of the mobile gas incurred from building more complex models, it is important
flow. Third, the facies model, though again the most accurate, to identify an optimal heterogeneity resolution in suchmodels
is not significantly more accurate than the other upscaled which can provide adequate predictions in CO2 modeling.
models. The layeredmodel, though again the least accurate, is However, in simulating CO2 flow, gas migration and trapping
not as incorrect as it is simulated in BC 3. This is similar to are affected by multiple processes, e.g., viscous flow from
what is observed in simulations under BC 2. These results active engineering (i.e., injection and production), mobility
suggest that brine production overrides the effect of hetero- effect (i.e., different relative permeabilities and fluid viscos-
geneity, making the reservoir appear more homogeneous. In ities), gravity segregation, formation heterogeneity, chemical
other words, the importance of representing heterogeneity reactions (only dissolution is modeled here), and site bound-
with higher geologic resolution and realism (i.e., facies ary condition. Since these processes operate and interact over
model), which is shown to be important in BC 3 simulations, varying spatial and temporal scales, the assessment of model
is significantly dampened. In BC 3, without brine production, complexity (i.e., identification of an appropriately simple site
how heterogeneity is represented is clearly important in model) is not straightforward. In this study, multiscale per-
determining the relative accuracy of the upscaled models. meability upscaling is combined with a sensitivity study on
3.3.5. CO Mass Profiles model boundary condition to identify an optimal model2
[ ] For the last performance metric, dissolved, mobile, complexity in developing a reservoir model for CO2 simula-51
and trapped CO computed by all models are shown over tion. Several insights are gained and are discussed as follows:2
time, under BC 2, 3, and 4 (Figure 10). In all models, with [56] In permeability upscaling, equivalent permeability
increasing time, more gas dissolves into brine, as expected. (one obtained for each unit of the upscaledmodels) is found to
During injection, more gas is dissolved per unit time than be largely insensitive to the number of flow experiments used
during monitoring (observed change in slope at the end of and the associated flow patterns. Simple linear floods appear
injection), likely due to the fact that per unit time, more grid to provide sufficient information for its calculation without
cells are swept during injection. More cells are contacted by the need to introduce complex flow patterns. This is consis-
gas, thus more dissolution. Once injection ceases, fluid tent with the results of Zhang et al. [2006]. This could be due
velocity becomes smaller (dominated by the slow and upward to the fact that when the upscaling domain is large compared
gravity segregation), and fewer cells are contacted by gas per to the characteristic size of heterogeneity (i.e., lnk correlation
unit time. Though gas dissolution is also affected by reservoir ranges), equivalent permeability can approach an effective
pressure, the level of pressure variation does not significantly value that is independent of the boundary condition used in
impact the rate of gas dissolution, e.g., during injection, upscaling [Renard and de Marsily, 1997]. Thus the upscaling
average pressure is dropping in BC 2, rising in BC 3, and domain in this study is suitably large for the calculation of a
maintaining a stable value in BC 4. Flow effect appears to stable k*. In Nugget sandstone variogrammodeling of facies‐
dominate gas dissolution rate, compared to the effect of specific generally yields correlation ranges that are
pressure variation. smaller than the size of the encompassing facies unit [Li et al.,
[52] The trapped gas also grows over time, as CO migrates 2011]. Correlation ranges of lnk are also smaller, due to the2
continuously in lateral and vertical directions, resulting in linear − log(k) transforms used to populate k from (one for
residual trapping. Additional trapping also occurs due to the each facies). For the layered and formation models, the size of
nonnegligible critical gas saturation assigned to the CO their units is even larger. In addition, unlike past research2
relative permeability drainage curve. Growth rate of the where model domains for upscaling were contiguous [Zhang
trapped gas is higher during injection. At this stage, higher et al., 2006], in this study, the shape of the facies units is
fluid velocity in the reservoir again results in more cells being extremely discontinuous, reflecting the geostatistical algo-
swept by gas per unit time. After the injection ceases, most rithms used in facies modeling (Sequential Indicator Simu-
trapping likely occurs at the trailing edge of the plume during lation). Still, k* of the facies units are stable and physically
the slow upward migration. reasonable, i.e., error characteristics in single‐phase verifi-
cation tests are consistent with those observed in past work.
18 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Figure 11. BC 3, open‐no‐brine, with two injectors (producers are inactivated). Gas‐phase CO2 satura-
tion. (left) End of injection (50 years) and (right) end of monitoring. First row is FHRM; second row is facies
model; third row is layered model; fourth row is formation model.
[57] In CO2 modeling, for the spatial and temporal scales is affected by brine production; the importance of detailed
considered, without brine production, optimal complexity of heterogeneity resolution is weakened when the viscous force
the upscaled model is found to depend on the prediction is strengthened in relation to the gravity force, i.e., increased
metric of interest. The facies model is the most accurate at lateral pressure gradient established by brine production in
capturing plume shape, fluid pressure, and CO2mass profiles, addition to that created by CO2 injection. Thus when flow is
while the formation model is adequate for pressure predic- more viscous dominated, variability of the predictions by the
tion. The layered model is not accurate for predicting most of upscaled models becomes smaller and the predictions are
the performance metrics, suggesting that both heterogeneity more accurate, suggesting a subtle but important interplay
resolution and homogenization strategy are important to between heterogeneity resolution, fluid driving forces, and
obtaining accurate predictions by the simple models. In terms model predictions. Brine production, besides being useful
of resolution, facies rather than the formation model is for controlling fluid pressure, has the potential benefit of
required if detailed plume shape is of interest. However, the offsetting conceptual model uncertainty. This observation has
layered configuration, despite its higher resolution than the implications for modeling CO2 storage in data‐poor settings,
formation model, is less accurate. Heterogeneity in the facies where an efficient and cost‐effective strategy is needed to build
model is characterized by an underlying unimodal lognormal a site model. For example, if brine production is employed, a
k. In the layered and formations models, for some units, simple model will perhaps suffice for the predictions of both
heterogeneity is homogenized over multimodal distributions. plume shape and reservoir pressure.
This suggests that k* obtained for model units with unimodal [59] The above observations, however, cannot be gener-
distributions is more accurate for multiphase prediction. alized easily, since they are based on one FHRM with a fixed
[58] Boundary condition is found to impact not only for- heterogeneity pattern and variance. With a single injector,
mation fluid pressure but also howmuch CO2 can be injected. very long injection time was used, but realistic cases will
If the reservoir is compartmentalized, brine production can likely use multiple injectors operating over shorter times. To
control and modulate pressure buildup as well as enhance test the results under more realistic conditions, additional
CO2 injectivity. If the reservoir is not compartmentalized, simulations are conducted with BC 3 (no brine production)
brine production does not significantly impact pressure but and BC 4 (brine production). Using the same constraints (e.g.,
can affect the direction of mobile gas flow. Importantly, the injector BHP, voidage replacement if applicable), following a
adequacy of an upscaled model in predicting CO2 storage streamline analysis in the FHRM to identify main reservoir
19 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Figure 12. BC 4, open‐with‐brine, with two injectors and two producers. Gas‐phase CO2 saturation. (left)
End of injection (50 years) and (right) end of monitoring. First row is FHRM; second row is facies model;
third row is layered model; fourth row is formation model.
connectivity, two injectors and two producers were sited in to the enhanced viscous flow. Again, brine production acts to
the model (Figure 11 and Figure 12). The injectors are located reduce the variability in the predictions made by the simple
in the central region, while the producers are located near the models. This result appears to be independent of the well
north and south boundaries. The producers are inactivated design and the time scale of injection.
in BC 3 runs. Over a period of 50 years, CO2 is injected [61] In this study, pure CO2 is injected into the reservoir.
at each well at 5 times the previous target rate, followed by Given the formation temperature and pressure condition,
3050 years of monitoring, for the same total simulation time fluid mobility ratio is fixed. No sensitivity study is conducted
of 3100 years. Over the 50‐year injection period, BC 3 runs by varying the fluid composition. However, mobility effect
achieved an average injection rate of 11 Mt/year, while BC 4 can also affect model complexity. For example, by improving
runs achieved an average rate of 17Mt/year. The largest coal‐ the mobility between the displacing and displaced phases
fired power plant near Shute Creek, the Jim Bridger plant, (e.g., injecting CO2 mixed with brine [Qi et al., 2009]), the
produces CO2 at approximately 18.5 Mt/year [Allis et al., importance of detailed heterogeneity on model predictions
2003]. Given the current boundary condition uncertainty could be reduced, i.e., the upscaled models become more
and BHP constraint, sequestration of this much CO2 will accurate. If, however, the mobility ratio is made to be unfa-
require two injection wells, with brine production needed if vorable (e.g., injecting CO2 to displace heavy oil), the effect
the reservoir is compartmentalized. of detailed heterogeneity may become more important, thus
[60] In BC 3, without brine production, gas plume appears performance of the upscaled models may degrade [Kumar
to be viscous dominated during injection, while gravity et al., 2005]. This is another reason why results of the pres-
becomes important postinjection (Figure 11). Similar to the ent study are difficult to generalize.
previous observation (Figure 8), the facies model is the most [62] Although brine production is shown to be unimportant
accurate while the layered model predicts incorrect physics for pressure control if the reservoir is “open,” that this result is
during postinjection when viscous force is reduced. In BC 4, obtained from modeling an infinitely large external aquifer,
the producers extract brine from both ends of the model which may not be realistic even though the Nugget formation
(Figure 12), inducing a different flow pattern from those occurs at the regional scale. In reality, pressure pulse from
previously simulated under the same BC (Figure 9). How- injection will eventually reach an aquifer boundary or sealing
ever, despite the changing flow pattern, the layered model has fault. Thus boundary conditions simulated here represent
predicted the correct postinjection gravity override, likely due end‐member scenarios. Before injecting CO2 at Shute Creek,
20 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
well tests of adequate durations are needed to detect and important prediction metrics of the heterogeneous model.
analyze the boundary, i.e., to determine the extent of reservoir This points to the possibility that simple models can be useful
compartment or distance to faults. These data, along with for CO2 storage simulation, rather than striving for the most
pressure and plume monitoring during the injection and complex and costly model. However, simulation results here
closure phases, can be used to help refine the boundary are only relevant to cases where a uniform relative perme-
condition assigned to the current models. ability prevails in the storage reservoir, which may not be
[63] Finally, in the storage schemes utilizing brine pro- realistic given the magnitude of variation of the intrinsic
duction, a large volume of brine can be produced and must permeability. Future work will address upscaling of relative
be disposed of. Potential options for brine disposal include permeabilities where multiple, facies‐based relative perme-
desalinization, beneficial use of the treated water, and dis- ability families can be assigned to the FHRM and suitable
posal of residual high‐density brine via reinjection [Surdam effective functions sought for the upscaled models. Further,
et al., 2008]. However, brine treatment and reinjection are though the facies model is identified as an optimal model
expensive, future work will seek means to optimize the (under brine production, perhaps the formation model), in
entire operation, thereby both enhancing storage security practice, equivalent permeability is not easily obtainable.
and reducing cost. Future work will develop calibration strategy to obtain
inverse permeabilities for the upscaled models that are con-
5. Conclusions sistent with the equivalent values.
[67] FHRM of this study was upscaled from a fine‐grid
[64] In this study, model complexity is evaluated by geostatistical model for computation efficiency. Future work
developing multiple conceptual models for the Nugget will directly model the fine gird using high‐performance
Sandstone, a deep saline aquifer in westernWyoming. A fully computing. With higher speed, a comprehensive uncertainty
heterogeneous reservoir model (FHRM) is first built, where analysis can be conducted. For example, in creating the
each grid cell is identified by multiple material tags (e.g., geostatistical model, uncertainties that have not been evalu-
facies identification). On the basis of these tags, permeability ated include kV /kH, geostatistical parameters in facies and upscaling is conducted to create three increasingly simpler modeling, and − log(k) transforms. These uncertainties
models, a facies model (four units), a layered model (three pertain to how accurately the FHRM represents the site
units), and a formation model (one unit). The accuracy of condition, which will be relevant when the goal of the sim-
upscaling is evaluated by conducting single‐phase verifica- ulation is to predict the actual performance during CO
tion tests in the upscaled models, with the FHRM providing 2injection. However, site‐specific analysis can only capture
reference predictions. Then, large‐scale CO2 storage simu- uncertainty at a single location. To further test insights of
lation is conducted in all models to assess the effect of this study, future work should address reservoirs at different
conceptual model uncertainty on predicting CO2 flow and depth, T/P regime, fluid types, permeability structure, het-
storage. At the injection site, since uncertainty exists in the erogeneity variance, etc. Experience with these systems will
characteristics of the reservoir compartment, end‐member lead to insights that can be more readily generalized.
boundary conditions are evaluated, whereby brine production
is introduced to control formation fluid pressure. The effect of
conceptual model uncertainty onmodel prediction is assessed [68] Acknowledgments. This material is based on work supported by
for each boundary condition, with results elucidating the the U.S. Department of Energy National Energy Technology Laboratory
interplays between fluid driving force, heterogeneity resolu- (award DE‐NT0004730) and the School of Energy Resources, University
tion, and boundary condition. of Wyoming. Partial funding for this study was provided by NSF grantEAR‐0838250. The views of the authors do not reflect those of the agencies.
[65] Results suggest that equivalent permeabilities We acknowledge software donation from Schlumberger, Inc. We acknowl-
obtained for the upscaled models are stable values indepen- edge the reviewers’ comments which helped improve the content and orga-
dent of the flow configurations used in upscaling. In CO nization of the manuscript.2
modeling, for the spatial and temporal scales considered,
without brine production, optimal complexity of the upscaled References
model depends on the prediction metric of interest; the facies Allis, R., T. Chidsey, C. Morgan, J. Moore, and S. White (2003), CO2
model is the most accurate at capturing plume shape, fluid sequestration potential beneath large power plants in Colorado plateau‐
pressure, and CO mass profiles, while the formation model southern RockyMountain region, USA, paper presented at SecondAnnual2
is adequate for pressure prediction. The layered model is not Conference on Carbon Sequestration, Dep. of Energy, Alexandria, Va.
Amaziane, B., A. Bourgeat, and J. Koebbe (1991), Numerical simulation
accurate for predicting most of the performance metrics. and homogenization of two‐phase flow in heterogeneous porous media,
Moreover, boundary condition impacts fluid pressure and the Transp. Porous Media, 6, 519–547.
amount of CO2 that can be injected. For the boundary con- Ataie‐Ashtiani, B., S. M. Hassanizadeh, and M. A. Celia (2002), Effects of
ditions tested, brine production can modulate fluid pressure, heterogeneities on capillary pressure‐saturation‐relative permeability
affect the direction of mobile gas flow, and influence the relationships, J. Contam. Hydrol., 56, 175–192.Bennion, B., and S. Bachu (2005), Relative permeability characteristics for
accuracy of the upscaled models. In particular, the impor- supercritical CO2 displacing water in a variety of potential sequestration
tance of detailed geologic resolution is weakened when vis- zones in the western Canada sedimentary basin, Rep. SPE 95547, Soc. of
cous force is strengthened in relation to gravity force. When Petrol. Eng., Richardson, Tex.
brine production is active, variability of the predictions by Bennion, B., and S. Bachu (2006a), The impact of interfacial tension and
pore‐size distribution/capillary pressure character on CO2 relative perme-the upscaled models becomes smaller and the predictions are ability at reservoir conditions in CO2 ‐brine systems, Rep. SPE 99325,
more accurate. Soc. of Petrol. Eng., Richardson, Tex.
[66] The method of this study is applicable to the study of Bennion, B., and S. Bachu (2006b), Supercritical CO2 and H2S‐brine drain-
other systems. Here, under the scenarios tested, an upscaled age and imbibitions relative permeability relationships for intergranular
model with adequate geological resolution is found to capture sandstone and carbonate formations, Rep. SPE 99326, Soc. of Petrol.Eng., Richardson, Tex.
21 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Birkholzer, J. T., and Q. Zhou (2009), Basin‐scale hydrogeologic impacts Juanes, R., E. Spiteri, F. J. Orr, and M. Blunt (2006), Impact of relative per-
of CO2 storage: Capacity and regulatory implications, Int. J. Greenhouse meability hysteresis on geological CO2 storage, Water Resour. Res., 42
Gas Control, 3, 745–756. W12418, doi:10.1029/2005WR004806.
Birkholzer, J. T., Q. Zhou, and C.‐F. Tsang (2009), Large‐scale impact of Kumar, A., M. Noh, G. Pope, K. Sepehrnoori, S. Bryant, and L. Lake
CO2 storage in deep saline aquifers: A sensitivity study on pressure (2004), Reservior simulation of CO2 storage in deep saline aquifers,
response in stratified systems, Int. J. Greenhouse Gas Control, 3, Rep. SPE 89343, Soc. of Petrol. Eng., Richardson, Tex.
181–194. Kumar, M., V. Hoang, C. Satik, and D. H. Rojas (2005), High mobility
Braun, C., R. Helmig, and S. Manthey (2005), Macro‐scale effective con- ratio water flood performance prediction: Challenges and new insights,
stitutive relationships for two‐phase flow processes in heterogeneous Rep. SPE 97671, Soc. of Petrol. Eng., Richardson, Tex.
porous media with emphasis on the relative permeability‐saturation rela- Law, D., and S. Bachu (1996), Hydrogeological and numerical analysis of
tionships, J. Contam. Hydrol., 76, 47–85. CO2 disposal in deep aquifers in the Alberta sedimentary basin, Energy
Carey, B., C. Christopher, and J. Gale (2007), Overview of the IEA Green- Convers. Manage., 37, 1167–1174.
house Gas R&D Programme’s Wellbore Integrity Network, paper pre- Li, S. Q., Y. Zhang, X. Zhang, and C. G. Du (2011), Geological modeling
sented at Technical Workshop on Geosequestration: Well Construction and fluid flow simulation of acid gas disposal in western Wyoming,
and Mechanical Integrity Testing, Environ. Prot. Agency, Los Alamos, AAPG Bull., in press.
N. M. (Available at http://www.epa.gov/ogwdw000/uic/pdfs/page_uic_ Lindquist, S. (1988), Practical characterization of eolian reservoirs for
albuquerque_presentations.pdf.) development: Nugget sandstone, Utah‐Wyoming thrust belt, J. Sediment.
Carlson, F. (1981), Simulation of relative pearmeability hysteresis to the Geol., 56, 315–339.
non‐wetting phase, Rep. SPE 10157, Soc. of Petrol. Eng., Richardson, Maldal, T., and I. Tappel (2004), CO2 underground storage for snohvit gas
Tex. field development, Energy, 29, 1403–1411, doi:10.1016/j.energy.2004.
Castellini, A., A. Chawathe, D. Larue, J. Landa, F. Jian, J. Toldi, and 03.074.
M. Chien (2003), What is relevant to flow? A comprehensive study McPherson, B., and B. Cole (2000), Multiphase CO2 flow, transport and
using a shallow marine reservoir, Rep. SPE 79669, Soc. of Petrol. sequestration in the Powder River Basin, Wyoming, U.S.A, J. Geochem.
Eng., Richardson, Tex. Explor., 69, 65–69.
Coburn, T. C., P. Corbett, A. Datta‐Gupta, J. Jensen, M. Kelkar, D. Oliver, Michael, K., G. Allinson, A. Golab, S. Sharma, and V. Shulakova (2009a),
and C. D. White (2006), A virtual conversation on the impact of geosta- CO2 storage in saline aquifers, II, Experience from existing storage
tistics in petroleum, production, and reservoir engineering, AAPG operations, Energy Proc., 1, 1973–1980.
Comput. Appl. Geol., 5, 23–33. Michael, K., M. Arnot, P. Cook, J. Ennis‐King, R. Funnell, J. Kaldi,
Das, D. B., S. M. Hassanizadeh, and B. E. Rotter (2004), A numerical study D. Kirste, and L. Paterson (2009b), CO2 storage in saline aquifers, I,
of micor‐heterogeneity effects on upscaled properties of two‐phase flow Current state of scientific knowledge, Energy Proc., 1, 3197–3204.
in porous media, Transp. Porous Media, 56, 329–350. Michael, K., A. Golab, V. Shulakova, J. Ennis‐King, G. Allinson, S. Sharma,
Doughty, C., and K. Pruess (2004), Modeling supercritical carbon dioxide and T. Aiken (2010), Geological storage of CO2 in saline aquifers—A
injection in heterogeneous porous media, Vadoes Zone J., 3, 837–847. review of the experience from existing storage operations, Int. J. Green-
Eaton, B. (1969), Fracture gradient prediction and its application in oilfield house Gas Control, 4, 659–667.
operations., J. Petrol. Technol., 42, 1353–1360. Milliken, W., M. Levy, S. Strebelle, and Y. Zhang (2007), The effect of
Elfenbein, C., P. Ringrose, and M. Christie (2005), Small‐scale reservoir geologic parameters and uncertainties on subsurface flow: Deepwater
modeling tool optimizes recovery offshore norway, J. World Oil, 226, depositional systems, Rep. SPE 109950, Soc. of Petrol. Eng., Richardson,
45–59. Tex.
Ertekin, T., J. Abou‐Kassem, and G. R. King (2001), Basic Applied Reser- Nordbotten, J., M. Celia, and S. Bachu (2005), Injection and storage of CO2
voir Simulation, Soc. of Petrol. Eng., Richardson, Tex. in deep saline aquifers: Analytical solution for CO2 plume evolution dur-
Flett, M., R. Gurton, and I. Taggart (2004), The function of gas‐water rel- ing injection, Trans. Porous Media, 58, 339–360, doi:10.1007/s11242-
ative permeability hysteresis in the sequestration of carbon dioxide in 004-0670-9.
saline formations, Rep. SPE 88485, Soc. of Petrol. Eng., Richardson, Obi, E., and M. Blunt (2006), Streamline‐based simulation of carbon diox-
Tex. ide storage in a north sea aquifer, Water Resour. Res., 42, W03414,
Flett, M., R. Gurton, and G. Weir (2007), Heterogeneous saline formations doi:10.1029/2004WR003347.
for carbon dioxide disposal: Impact of varying heterogeneity on contain- Orr, F. M. (2009), Onshore geological storage of CO2, Science, 325,
ment and trapping, J. Petrol. Sci. Eng., 57, 106–118, doi:10.1016/j. 1656–1658.
petrol.2006.08.016. Pickup, G., K. Stephen, J. Ma, P. Zhang, and J. Clark (2005), Multi‐stage
Flett, M., et al. (2008), Gorgon project: Subsurface evaluation of carbon upscaling: Selection of suitable methods, Trans. Porous Media, 58,
dioxide disposal under barrows island, Rep. SPE 116372, Soc. of Petrol. 191–216.
Eng., Richardson, Tex. Primera, A., W. Sifuentes, and N. Rodriguez (2009), CO2 injection and
Friedmann, F., A. Chawathe, and D. Larue (2003), Assessing uncertainty in storage: A new approach using integrated asset modeling, Rep. SPE
channelized reservoirs using experimental designs, Rep. SPE 85117, 121970, Soc. of Petrol. Eng., Richardson, Tex.
Soc. of Petrol. Eng., Richardson, Tex. Pruess, K., J. Garcia, T. Kovscek, C. Oldenburg, J. Rutqvist, C. Steefel,
Hesse, M., H. Tchelepi, and F. J. Orr (2006), Scaling analysis of the migra- and T. Xu (2004), Code intercomparison builds confidence in numerical
tion of CO2 in saline aquifers, Rep. SPE 102796, Soc. of Petrol. Eng., simulation models for geologic disposal of CO2, Energy, 29, 1431–1444,
Richardson, Tex. doi:10.1016/j.energy.2004.03.077.
Hollowway, S., A. Chadwick, E. Lindeberg, I. Czernichowski‐Lauriol, and Qi, R., T. LaForce, and M. Blunt (2007), Design of carbon dioxide storage
R. Arts (2004), Best practice manual from Saline Aquifer CO2 Storage in a North Sea aquifer using streamline‐based simulation, Rep. SPE
Project, Statoil, Strondheim, Norway. (Available at http://www.co2store. 109905, Soc. of Petrol. Eng., Richardson, Tex.
org/TEK/FOT/SVG03178.nsf/Attachments/SACSBestPractiseManual. Qi, R., T. LaForce, and M. Blunt (2009), Design of carbon dioxide storage
pdf/FILE/SACSBestPractiseManual.pdf.) in aquifers, Int. J. Greenhouse Gas Control, 3, 195–205.
Ide, S., K. Jessen, and F. J. Orr (2007), Storage of CO2 in saline aquifers: Renard, P., and G. de Marsily (1997), Calculating equivalent permeability:
Effects of gravity, viscous, and capillary forces on amount and timing of A review, Adv. Water Resour., 20, 253–278.
trapping, Int. J. Greenhouse Gas Control, 1, 481–491. Sáez, A. E., C. J. Otero, and I. Rusinek (1989), The effective homogeneous
Intergovernmental Panel onClimate Change (2005), Special report on carbon behavior of heterogeneous porous media, Trans. Porous Media, 4,
dioxide capture and storage, Geneva. (Available at http://arch.rivm.nl/env/ 213–238.
int/ipcc/pagesmedia/SRCCSfinal/IPCCSpecialReportonCarbondioxide- Sanchez‐Vila, X., A. Guadagnini, and J. Carrera (2006), Representative
CaptureandStorage.htm.) hydraulic conductivities in saturated groundwater flow, Rev. Geophys.,
Johnson, J., J. Nitao, C. I. Steefel, and K. Knauss (2001), Reactive transport 44, RG3002, doi:10.1029/2005RG000169.
modeling of geologic CO2 sequestration in saline aquifers: The influence Sasaki, K., T. Fujii, Y. Niibori, T. Ito, and T. Hashida (2008), Numerical
of intra‐aquifer shales and the relative effectiveness of structural, solubil- simulation of supercritical CO2 injection into subsurface rock masses,
ity, and mineral trapping during prograde and retrograde sequestration, Energy Conversion Manage., 49, 54–61, doi:10.1016/j.enconman.2007.
Lawrence Livermore Natl. Lab., Livermore, Calif. (Available at http:// 05.015.
www.netl.doe.gov/publications/proceedings/01/carbonseq/P28.pdf.) Schlumberger (2009), Eclipse technical description manual, Abingdon,
U. K.
22 of 23
W05534 LI ET AL.: A STUDY OF CONCEPTUAL MODEL UNCERTAINTY W05534
Stauffer, P., H. Viswanathan, R. Pawar, and G. Guthrie (2009), A system Zhang, Y., and C. W. Gable (2008), Two‐scale modeling of solute
model for geologic sequestration of carbon dioxide, Environ. Sci. transport in an experimental stratigraphy, J. Hydrol., 348, 395–411,
Technol., 43, 565–570. doi:10.1016/j.jhydrol.2007.10.017.
Surdam, R. C., K. Clarey, R. D. Bentley, S. J. E., and Z. Jiao (2008), Desa- Zhang, Y., C. W. Gable, and M. Person (2006), Equivalent hydraulic con-
linization project feasibility, in Water Associated with Coal Beds in ductivity of an experimental stratigraphy ‐ Implications for basin‐scale
Wyoming’s Powder River Basin, Explor. Mem., vol. 2, edited by D. A. flow simulations, Water Resour. Res., 42, W05404, doi:10.1029/
Copeland and M. L. Ewald, pp. 295–311, Wyo. State Geol. Surv., 2005WR004720.
Laramie, Wyo. Zhang, Y., B. Liu, and C. Gable (2011), Three‐dimensional hydraulic
Surdam, R. C., Z. Jiao, P. Stauffer, and T. Miller (2009), An Integrated conductivity upscaling for hierarchical sedimentary deposits at multi-
Strategy for Carbon Management Combining Geological CO2 Seques- ple scales based on a synthetic aquifer, Transp. Porous Media, 87(3),
tration, Displacement Fluid Production, and Water Treatment, Chal- 717–737, doi:10.1007/s11242-010-9711-8.
lenges in Geol. Resour. Dev., vol. 8, Wyo. State Geol. Surv., Laramie, Zhou, Q., J. Birkholzer, T. C. F., and J. Rutqvist (2008), A method for
Wyo. quick assessment of CO2 storage capacity in closed and semi‐closed
Wen, X.‐H., and J. Gomez‐Hernandez (1996), Upscaling hydraulic saline formations, Int. J. Greenhouse Gas Control, 2, 626–639.
conductivities in heterogeneous media: An overview, J. Hydrol., 183, Zhou, Q., J. T. Birkholzer, E. Mehnert, Y.‐F. Lin, and K. Zhang (2010),
ix–xxxii. Modeling basin‐ and plume‐scale processes of CO2 storage for full‐scale
Wen, X.‐H., L. Durlofsky, and M. Edwards (2003), Use of border regions deployment, Ground Water, 48, 494–514.
for improved permeability upscaling, Math. Geol., 35, 521–547.
Wolery, T. J., R. Aines, Y. Hao, W. Bourcier, T. Wolfe, and C. Haussman S. Li and Y. Zhang, Department of Geology and Geophysics, University of
(2009), Fresh water generation from aquifer pressured carbon storage: Wyoming, 1000 E. University Ave., Laramie, WY 82071, USA. (sli2@uwyo.
Annual report FY09, Rep. lLNL‐TR‐420857, Lawrence Livermore Natl. edu; yzhang9@uwyo.edu)
Lab., Livermore, Calif. X. Zhang, Schlumberger Information Solutions, Schlumberger, 5599
San Felipe, Suite 1700, Houston, TX 77056, USA.
23 of 23
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.