Dissolution Trapping of Carbon Dioxide in ... - ACS Publications

Jun 9, 2017 - School of Earth Sciences, The Ohio State University, Columbus, Ohio ... of Earth and Environmental Sciences, Wright State University, Da...
1 downloads 0 Views 2MB Size
Subscriber access provided by UNIV OF ARIZONA

Article

Dissolution Trapping of Carbon Dioxide in Heterogeneous Aquifers Mohamad Reza Soltanian, Mohammad Amin Amooie, Naum Gershenzon, Zhenxue Dai, Robert Ritzi, Fengyang Xiong, David R. Cole, and Joachim Moortgat Environ. Sci. Technol., Just Accepted Manuscript • Publication Date (Web): 09 Jun 2017 Downloaded from http://pubs.acs.org on June 12, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Environmental Science & Technology is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33

Environmental Science & Technology

Dissolution Trapping of Carbon Dioxide in Heterogeneous Aquifers Mohamad Reza Soltanian,†,‡ Mohammad Amin Amooie,† Naum Gershenzon,¶ Zhenxue Dai,§,k Robert Ritzi,¶ Fengyang Xiong,† David Cole,† and Joachim Moortgat∗,†,⊥ †School of Earth Sciences, The Ohio State University, Columbus, OH 1

[email protected], Exponent, 1055 East Colorado Boulevard, Suite 500 Pasadena, CA 91106 ¶Department of Earth and Environmental Sciences, Wright State University, Dayton, OH §Los Alamos National Laboratory, Los Alamos, New Mexico, USA kCollege of Construction Engineering, Jilin University, Changchun, China ⊥Correspondence: 303 Mendenhall Laboratory, 125 South Oval Mall, Columbus, OH 43210, Tel: +1-614-688-2140, Fax: +1-614-292-7688 E-mail: [email protected]

2

Abstract

3

The geologic architecture in sedimentary reservoirs affects the behavior of density-

4

driven flow and the dispersion of CO2 -rich brine. The spatial organization and con-

5

nectivity of facies types play an important role. Low-permeability facies may sup-

6

press fingering and reduce vertical spreading, but may also increase transverse mix-

7

ing. This is more pronounced when geologic structures create preferential flow path-

8

ways through connected facies types. We perform high-resolution simulations of three-

9

dimensional (3D) heterogeneous formations whose connectivity cannot be represented in

1

ACS Paragon Plus Environment

Environmental Science & Technology

10

two-dimensional models consistent with percolation theory. This work focuses on the

11

importance of 3D facies-based heterogeneity and connectivity on advection-diffusion

12

transport of dissolved CO2 . Because the dissolution of CO2 and the subsequent density

13

increase of brine are the driving force for gravitational instabilities, we model the phase

14

behavior with the accurate cubic-plus-association equation-of-state, which accounts for

15

the self-association of polar water molecules and the cross-association between CO2 and

16

water. Our results elucidate how the spatial organization of facies affects the dynamics

17

of CO2 convective mixing. Scaling relations for the evolution of a global dispersion-

18

width provide insights that can be universally applied. The results suggest that the

19

long-term evolution and scaling of dispersion are surprisingly similar for homogeneous

20

and (binary and multiscale) heterogeneous porous media.

21

Introduction

22

Carbon dioxide (CO2 ) injection in geologic formations is a promising short- and long-term

23

strategy to slow the increase in atmospheric CO2 concentration. 1–5 To assure decision makers

24

and the public of the feasibility, economics, and storage permanence of Carbon Capture

25

and Storage (CCS), a thorough understanding of the basic science behind CCS is critical.

26

However, comprehensive assessment of CCS is challenging due to the complex physical and

27

chemical processes controlling the transport of CO2 in geologic reservoirs. This includes

28

CO2 phase behavior and different types of trapping mechanisms (e.g., residual, dissolution,

29

structural, and mineral). 6 In this work, we focus on dissolution or solubility trapping of CO2

30

in brine, which is believed to be an intermediate-term (0–1000 years) mechanism. The focus

31

is on the dissolution mechanisms taking place in the aquifer after a plume reaches the cap

32

rock. Residual, 7 structural, 8 and mineral trapping are not considered.

33

Solubility trapping is enhanced when gravitational instabilities (fingering) are triggered

34

by a local increase in brine density as CO2 dissolves into brine in the top of an aquifer. 9–13

35

Whether the interface between CO2 -rich brine and fresh brine becomes gravitationally un2

ACS Paragon Plus Environment

Page 2 of 33

Page 3 of 33

Environmental Science & Technology

36

stable depends on the ratio of advection to molecular diffusion. When a fingering instability

37

is triggered, dissolved CO2 is mixed throughout the aquifer at advective time-scales, which

38

can be much faster than diffusive transport. This fast spreading towards higher depths im-

39

proves the storage capacity of a given aquifer and decreases the leakage risk in case of cap

40

rock failure.

41

Heterogeneity in petrophysical properties (e.g., permeability and porosity) has a signif-

42

icant impact on density-driven flow of CO2 . 11,13–15 Both rock and fluid properties directly

43

affect the characteristics of fingers including the critical time for the onset of the instability

44

(tc ) and a critical wavelength of fingers (λc ). The results of linear stability analyses, as well

45

as numerical simulations show that: 16–21

46

47

φµ , k∆ρg 2  φµ = c2 D , k∆ρg

λc = c 1 D

(1)

tc

(2)

48

where µ is brine viscosity, φ is porosity, D is the diffusion coefficient for CO2 in brine, k is

49

permeability, ∆ρ is the maximum density increase of the aqueous phase upon CO2 disso-

50

lution, and g is the gravitational acceleration. The c1 and c2 are proportionality constants

51

that vary by orders of magnitudes in the cited references.

52

Convective mixing and solubility trapping can be enhanced by subsurface heterogene-

53

ity. 11,14 Geological heterogeneity is controlled by the spatial organization of facies types which

54

are three-dimensional (3D) bodies of rocks, whose differentiation provides a useful concep-

55

tual framework to characterize heterogeneity in attributes of interest (e.g., permeability). 22

56

The architecture of facies types may exhibit sharp discontinuities across boundaries between

57

depositional features, for instance across the sandstone-shale contacts that are common in

58

fluvial deposits. In fact, some potential CO2 reservoirs have sedimentary architecture re-

59

flecting fluvial deposition with complex architecture created by depositional processes. Such

60

processes often result in variable connectivity and highly tortuous flow pathways. 6,23–25 Well-

3

ACS Paragon Plus Environment

Environmental Science & Technology

61

known examples include the Victor interval of the Ivishak Formation, 8,26 lower Mt. Simon

62

Sandstone, 25 the lower Tuscaloosa Formation, 23,27 the lower Paaratte Formation, 27,28 and

63

the Morrow Sandstone. 29

64

Facies connectivity significantly affects CO2 migration by creating preferential flow path-

65

ways. 8,13,30–32 Connected zones of high-permeability facies exhibit shorter critical times and

66

wavelengths, and thus experience more pronounced fingering. In this work we perform 3D

67

simulations to (1) represent 3D geologic structure and connectivity, (2) understand the fluid

68

dynamics of gravitational fingering when 3D heterogeneity is considered, and (3) provide

69

scaling relations to predict the long-term dissolution trapping of CO2 . Earlier 3D simulation

70

studies of CO2 convective mixing have only considered homogeneous media. 19

71

The volume proportion of each facies is the most important factor controlling the con-

72

nectivity. 33 Percolation theory suggests that randomly placed cells of one type have at least

73

one connected tortuous pathway through a cluster that spans opposing boundaries if the

74

volume proportion of that type is above the 3D percolation threshold of 31.16%. 33,34 If the

75

spatial organization of a facies type is represented within a finite domain, a tortuous fully

76

connected pathway is formed at even lower volume proportions. 35 In this work we quantify

77

connectivity by studying connected clusters of high-permeability cells within two different

78

types of geologic models using information on how geologic structure affects the percolation

79

principle. 35–38 In the first model we create synthetic binary distributions of facies types using

80

a Markov Chain approach. Facies architectures are represented by physically quantifiable

81

attributes such as mean lengths and volume proportions. The second geologic model is a

82

multiscale complex representation of a fluvial architecture, which consists of several facies

83

types spanning from the cm to the hundred meter scale. In both models, we create different

84

realizations, each with a different volume proportion of high-permeability facies.

85

Our main objective is to study the impact of three-dimensional facies connectivity on

86

the density-driven convective mixing of dissolved CO2 . By using models for which the con-

87

nectivity of facies types are known, we aim to quantify the role of connectivity on solubility

4

ACS Paragon Plus Environment

Page 4 of 33

Page 5 of 33

Environmental Science & Technology

88

trapping in general, and the dynamics of CO2 mixing in particular.

89

Methodology

90

Heterogeneity structures and connectivity

91

We consider (1) a synthetic binary heterogeneity model, and (2) a heterogeneity model that

92

reflects a multiscale fluvial architecture. In both models we vary the volume proportion

93

of high-permeability facies to study connectivity. There are several approaches to study

94

connectivity including the use of (1) flow- and transport-based measures, 39 (2) critical path

95

analysis, 40 and (3) facies-based approaches linked with percolation theory. 8,35,36 Here we use

96

the facies-based approach and define the cell connectivity in three dimensions.

97

Binary heterogeneity

98

We use generic heterogeneity models in binary media with spatial distribution of high-

99

permeability facies within a background of lower permeability facies. We model a 50 m ×

100

50 m × 50 m domain, discretized by a fine 80 × 80 × 80 grid with 512,000 grid cells.

101

We show the fundamental importance of 3D heterogeneities arising from facies archi-

102

tecture. We choose a parsimonious set of two facies types. 41,42 The multifacies model is

103

described in the next section. Different structures are created using a Markov Chain model

104

and indicator simulation with quenching in T–PROGS. 43 Rather than focusing on any sin-

105

gle site or data set, we represent general permeability distributions consisting of high- and

106

low-permeability facies. The volume proportions of high-permeability facies are increased

107

from 5 to 16, 28, and 40% (Figure S1 in the Supporting Information shows realizations with

108

16 and 28%). In all binary cases, we use an isotropic integral scale (λ) of 10 m. The effect of

109

small-scale anisotropy within grid cells is presented as a special case in the Results and Dis-

110

cussion. A high and a low permeability of 5, 000 md and 500 md are assigned to the two facies

111

types, respectively. The porosity of high- and low-permeability facies are calculated as 0.34 5

ACS Paragon Plus Environment

Environmental Science & Technology

112

and 0.13, respectively, using k = aφb where a and b are constants. 44,45 We do not consider

113

permeability variations within each facies but assume that the spatial organization of facies

114

types is the major source of heterogeneity. We choose domain size, cell size, and permeability

115

values based on prior experience and Equation (1), which shows that λc is proportional to

116

1/k. Thus, the cell size should be small enough to capture the smallest fingers within high

117

permeability facies. At the same time, the domain size should be large enough compared to

118

λc so that the domain encompasses the largest fingers within low-permeability facies. For

119

comparison, we also create homogeneous domains with permeabilities and porosities based

120

on the means of the heterogeneous domains (Table S1).

121

We use a search algorithm to determine how high-permeability facies are connected in

122

the 3D heterogeneous binary model. As an example, for the volume proportion of 28% this

123

algorithm determined that 96% of high-permeability facies (27% of all cells) are connected

124

in a single cluster that spans between any pair of opposing domain boundaries, creating

125

a percolating cluster (see Figure S1d). It is obvious that the volume proportion of 28%

126

is below the volume proportion of 31.16% known as percolation threshold for a random,

127

infinite cubic lattice. However, the structure added to the finite-domain model reduced the

128

threshold, as is consistent with previous studies. 35,38,46 The percentage of connected cells

129

for volume proportions of 5, 16, and 40% are ≪ 1, 1.4, and 39.6%, respectively. There

130

is no connected percolating cluster for volume proportions of 5 and 16%. However, there

131

is a connected percolating cluster for volume proportions of both 28 and 40%. Note that

132

two-dimensional (2D) slices from a 3D model do not have the connectivity of the 3D domain

133

they are extracted from, as per percolation theory. 46

134

Multiscale heterogeneity

135

Many candidate reservoirs for geologic CO2 sequestration have fluvial depositional structures.

136

Recent studies have led to new conceptual and quantitative models for such architectures

137

over a range of spatial scales that are relevant to fluid flow in oil and deep brine reservoirs. 47

6

ACS Paragon Plus Environment

Page 6 of 33

Page 7 of 33

Environmental Science & Technology

138

This approach to characterizing fluvial architectures, including scales of stratification, volume

139

proportions, relative sizes, and 3D geometry for facies types at each scale, was incorporated

140

into a multiscale heterogeneity model by Guin et al. 38 and Ramanathan et al. 37

141

The multiscale heterogeneity model 8,25,41,48 represents fluvial deposits dominated by sandy

142

gravel and the associated conglomeratic reservoirs. At the largest scale, the model comprise

143

active channels and compound bars. The compound bars, in turn, comprise unit bars.

144

Within the unit bars there are the smallest scale cross-strata with different textures. Im-

145

portantly, there are high-permeability sets of cross-stratified open-framework conglomerate

146

(OFC), which are known to create preferential flow pathways. The individual OFC cross-sets

147

are represented in the model, and clusters of continuously connected OFC cells create pref-

148

erential flow pathways. Guin et al. 38 found that tortuous spanning pathways exist through

149

connected OFC cells if OFC volume proportions were at or above 20%, independent of other

150

input parameters. The number, size, and orientation of such OFC clusters change with the

151

volume proportions of each facies. Even for a given volume proportion, these properties

152

change across hierarchical levels. 37,38 For instance, connected OFC cells at the smallest scale

153

(level I) form flow pathways that vertically span (i.e., connect across) a single unit bar at the

154

next scale (level II). Connections across unit bars at level II enhance lateral branching and

155

the many clusters within unit bars connect into a smaller number of larger clusters at the

156

scale of multiple unit bars. At the scale of a compound bar deposit (level III), these clusters

157

are typically connected into one or two large clusters that span all the domain boundaries.

158

We use this complex heterogeneity model to investigate the effect of OFC connectivity on

159

the dynamics of density-driven flow of CO2 .

160

We generate three realizations in which the volume proportions of OFC are 19, 24, and

161

28%. The volume proportion controls the percentage of OFC cells that are connected in

162

clusters to create preferential flow pathways. 8 Such spanning clusters may develop at OFC

163

proportions below the well-known percolation threshold of 31.16%. The search algorithm

164

described in the previous section is used to identify connected clusters of OFC cells, and

7

ACS Paragon Plus Environment

Environmental Science & Technology

165

to quantify their attributes. The proportion of all connected cells are 2.5, 17.4, and 25.6%,

166

respectively, for volume proportions of 19, 24, and 28% OFC cells. There are spanning

167

clusters in the cases with 24 and 28% OFC, but not for 19%. Therefore, the multiscale

168

heterogeneity models range from one with no spanning cluster (19%) to one with about 90%

169

connected OFC in one spanning cluster (28%).

170

The permeability fields of cases with 19 and 28% volume proportions are illustrated in

171

Figure S2. The domain size is 100 m × 100 m × 15 m, discretized by a 2 m × 2 m × 0.05 m

172

grid-cell size with 750,000 grid cells. For each facies type we use lognormal permeability

173

distributions as found in natural deposits. The permeability of each cell is assigned from

174

the probability density functions (PDFs) defined for each facies type (see Figure 3, panels

175

on the right, in Ramanathan et al. 37 ). The geometric mean of permeabilities for 28, 24, and

176

19% cases are 720, 891, and 1051 md, respectively.

177

The thickness-to-length ratio in the geometry of facies is significantly less than unity at

178

all spatial scales. Therefore, the heterogeneity gives rise to an effective anisotropy at the

179

scale of the sedimentary facies types. Anisotropy within the 5 cm thickness of grid cells

180

is not considered (as in Gershenzon et al. 8 ). Note that the grid-cell size is defined by the

181

characteristic size of the smallest scale facies types in the model architecture. Larger grid-

182

cell sizes cannot capture the small-scale heterogeneities. Therefore, to make the problem

183

computationally feasible, there is a trade-off between grid-cell size and domain size.

184

Numerical methods and boundary conditions

185

We use numerical methods presented in earlier works. 21,49–56 The governing equations are

186

presented in Supporting Information, Text S1. We consider a fluid mixture with two compo-

187

nents (water and CO2 ). Note that we consider fluid compression due to the pressure build-up

188

during CO2 injection into an aquifer with impermeable boundaries and volume swelling of

189

brine due to CO2 dissolution. This is in contrast to previous studies on gravitational fingering

190

in which an incompressible aqueous phase is considered. 8

ACS Paragon Plus Environment

Page 8 of 33

Page 9 of 33

Environmental Science & Technology

191

To solve the governing equations we use higher-order finite element methods. Specifi-

192

cally, we adopt a second-order discontinuous Galerkin (DG) method to solve the transport

193

equation. 57 Using DG is beneficial in studying flow instabilities since it is able to capture

194

the small-scale onset of viscous and gravitational fingers on relatively coarse grid cells by

195

reducing numerical dispersion. This helps us to resolve λc with only few grid cells. 21 We

196

use the mixed hybrid finite element (MHFE) method to simultaneously solve for continuous

197

pressure and velocity fields. MHFE provides an accurate and globally continuous velocity

198

field for highly heterogeneous domains. This feature is critical given the sharp contrasts in

199

the permeability fields arising from the spatial organization of sedimentary facies.

200

Phase behavior is based on the cubic-plus-association (CPA) equation-of-state (EOS).

201

The CPA-EOS was developed to model the phase behavior of fluids that contain polar

202

molecules, such as water 58 or asphaltenes, 59 by using thermodynamic perturbation theory.

203

The CPA-EOS takes into account the self-association between water molecules, as well as

204

the (polarity-induced) cross-association between water and CO2 . 51,60,61 This approach is an

205

improvement over previous studies that relied on empirical correlations for the brine density

206

and Henry’s law for CO2 solubilities. We strictly enforce local thermodynamic equilibrium.

207

Also, the CPA-EOS has been shown to accurately reproduce measured densities of CO2 -

208

brine mixtures and the volume swelling due to CO2 dissolution. Volume swelling and the

209

associated movement of the interface between free CO2 in a gas cap and the underlying brine,

210

cannot be modeled by considering only the aqueous phase 17 and replacing the gas cap with

211

a constant-CO2 -composition Dirichlet boundary condition. 62

212

The temperature and the pressure at the bottom boundary are, respectively, 77 ◦ C and

213

100 bar. At these conditions, the aqueous phase density is ρw = 981 kg/m3 , the maximum

214

CO2 solubility is 1.7 mol%, and the CO2 -saturated aqueous phase density is ρ = 987 kg/m3

215

with ∆ρ = (ρ − ρw )/ρw = 0.6% (note that CO2 behaves as a supercritical fluid in this

216

pressure and temperature condition). The diffusion coefficient is D = 1.33 × 10−8 m2 s−1 .

217

We consider 3D bounded domains with impermeable no-flow Neumann conditions for all

9

ACS Paragon Plus Environment

Environmental Science & Technology

Page 10 of 33

218

boundaries. The domain is initially saturated with brine. CO2 is injected uniformly into

219

the top grid cells at a low rate (as a source term of 0.1% pore volume (PV) per year). The

220

domain remains in a single aqueous phase without a capillary transition zone. 63,64 Because

221

of the impermeable boundaries the pressure increases upon injection. We also simulated

222

a case with an open bottom boundary in which the pressure remains roughly constant. 64

223

However, the evolution of the dispersion width (and the mixing dynamics more generally) do

224

not appear to be sensitive to this boundary condition. The reason is that, while the pressure

225

evolution is quite different for the two different boundary conditions, the density difference

226

within the aqueous phase upon CO2 dissolution is nearly identical, and this is the driver for

227

convective mixing. Additional testing shows that the computational domain is sufficiently

228

wide that the lateral boundaries do not impact the temporal evolution of CO2 mixing.

229

The density-driven flow of CO2 is quantified by the global dispersion-width σz (the square

230

root of the spatial variance σz2 ) of the molar density of CO2 following Soltanian et al. 13 (see

231

also Agartan et al. 14 ). The choice of CO2 molar density, rather than composition, accounts

232

for the density change due to both the CO2 dissolution and the compressibility of aqueous

233

phase. The first three central moments of the distribution define the spatial variance σz2 of

234

the molar density of CO2 (C ≡ czCO2 ): σz2 (t)

N002 = − N000



N001 N000

2

(3)

235

where N001 /N000 expresses the coordinate location of the center of mole in the vertical (z)

236

direction. The Nijk is the ijk-th moment of the CO2 molar density distribution in space,

237

analogous to the moment of concentration distribution, 65 defined here as:

Nijk (t) =

Z

Lx 0

Z

Ly 0

Z

Lz

xi y j z k φC(x, y, z, t)dxdydz

(4)

0

238

with Lx , Ly , and Lz the domain size in x, y, and z directions, respectively. The total number

239

of moles of CO2 in the domain is physically represented by the zeroth spatial moment, N000 . 10

ACS Paragon Plus Environment

Page 11 of 33

Environmental Science & Technology

240

The N000 does not remain constant over time since CO2 is continuously injected into the

241

domain. The spatial variance is not sensitive to N000 increase since the first and second

242

moments are normalized by N000 . The variances in the x-direction and y-direction are

243

negligible due to the dominant downward migration of CO2 -rich brine. Fundamentally σz is

244

a measure of spreading and dispersion in porous media, however, it may also represent the

245

mixing state since more spreading usually leads to more mixing. 66

246

Results and Discussion

247

Binary media

248

We first study the dynamics of fingering in the synthetic binary media and their homogeneous

249

equivalents.

250

We provide Figure 1a that visually shows the profiles for mole fraction of CO2 after 3, 4.5,

251

6.5, 8.5, 11, and 13.5 years for homogeneous domains corresponding to binary media with

252

16 and 28% volume proportions of high-permeability facie types. The mean (homogeneous)

253

permeability and porosity increase for larger volume proportions of high-permeability facies

254

(see Table S1). It is clear that the onset of fingering occurs later for the homogeneous

255

case of 16% than that of 28% and there are more CO2 fingers initiated for the 28% volume

256

proportion. Using the values in Table S1 we find that λc (28%)/λc (16%) is about 0.8 and

257

tc (28%)/tc (16%) is about 0.65.

258

Phase behavior is a critical part of our simulations. A thermodynamically consistent

259

description of phase behavior results in a driving force for gravitational instabilities (∝ ∆ρ)

260

that evolves non-trivially upon CO2 dissolution (see Figure 5 in Soltanian et al. 13 ). This

261

evolution of ∆ρ results in four distinct flow regimes in our scaling analyses of homogeneous

262

simulations presented in Figure 2a, which is consistent with our novel findings for 2D simu-

263

lations in homogeneous media. 13 The alternating flow regimes are interpreted as follows.

264

The first flow regime corresponds to the well-known diffusive transport of CO2 before the 11

ACS Paragon Plus Environment

Environmental Science & Technology

ACS Paragon Plus Environment

Page 12 of 33

Page 13 of 33

Environmental Science & Technology

265

onset of gravitational fingers (see Figure 1a, panel 3 years for 16%). As shown in Figure 2a,

266

the σz exhibits a classical Fickian scaling of σz ∝ t0.5 . The front in this regime travels to a

267

diffusion penetration depth, which obeys a (Dt)0.5 Einstein-type scaling.

268

The second flow regime corresponds to the onset of fingering and can be recognized as

269

the transition to a convection-dominated regime, according to the sharp increase in σz in

270

Figure 2a. The dependence of tc and λc on k and φ imply that higher permeabilities result in

271

an earlier onset time with a larger number of small-scale fingers, as illustrated in Figure 1a.

272

The initial scaling of σz in Figure 2a for the fingering regime is σz ∝ t3.5 . This is also

273

consistent with our 2D simulation results. 13 Note that the onset time is earlier for larger

274

values of homogeneous permeability and porosity.

275

The second regime is associated with reduction in ∆ρ within the fast sinking fingers. 13

276

This is because the high rate of downward advection leads to stretching of CO2 -rich plume,

277

hence lowering the peak CO2 concentration and subsequently the maximum ρ. Meanwhile,

278

the stretching-induced sharp concentration gradients across the finger boundaries cause the

279

diffusion of CO2 into the ambient upwelling brine. This results in depletion of CO2 in the

280

sinking fingers, which further contributes to ∆ρ reduction. In a negative feedback loop, the

281

resulting reduced ∆ρ inhibits further stretching of fingers (note that ∆ρ is the driving force).

282

The second regime is therefore transient, and is followed by a third restoration regime.

283

The third regime is attributed to a fingering stagnation period, in which the depleted

284

fingers become replenished with new CO2 that is injected from the top boundary and dis-

285

solved in brine at a diffusive rate. We find diffusive scaling in this regime, mainly because of

286

persistent coalescence and merging of growing fingers resulting in coarsened fingers, as well

287

as the transverse diffusion from the CO2 -rich fingers into the brine as the CO2 concentra-

288

tions inside the fingers become higher again during the restoration. This regime transitions

289

into the final fully-developed convection once the ∆ρ, through constant CO2 injection, again

290

exceeds a critical limit that is required to drive the gravitational instabilities. 13

291

Finally, in the fourth regime the flow is again advection-dominated and shows a sharp

13

ACS Paragon Plus Environment

Environmental Science & Technology

292

increase in the growth rates of σz . This regime and its related scaling for σz are the most

293

important ones since they control the long-term evolution of σz , while the first three regimes

294

occur at shorter time-scales (∼10 years in our results). What defines the long-term fate of the

295

injected CO2 is the evolution of this fourth flow regime, the extent of which is determined by

296

the vertical size of the aquifer. Thick aquifers may be more effective in trapping CO2 through

297

solubility. It has been also shown that the capillary sealing efficiency of cap rocks –another

298

storage mechanism– is enhanced by CO2 sequestration at higher depths. 67 This suggests

299

that CCS projects may benefit from deeper formations, in spite of the higher costs involved.

300

Off-shore sedimentary aquifers, which can have thicknesses of the order of a kilometer, 71 are

301

attractive for the same reasons.

302

Note that such distinct flow regimes, especially the second and third regimes, are not

303

observed in models that assume incompressible flow and a constant pressure and CO2 -

304

composition at the top boundary of a brine-saturated domain. Under such conditions, the

305

maximum value for the density contrast (∆ρmax ) between pure and CO2 -saturated brine

306

remains constant. A constant value of ∆ρmax is also assumed in all linear stability analyses

307

that lead to expressions of the form Equations (1) and (2). Allowing for a variable ∆ρmax

308

results in the more complex dynamics that are captured in our simulations.

309

We also show in Figure S3 the mole fraction of CO2 after 1.5, 2.5, 4.5, 8.5, 11, and

310

13.5 years for binary heterogeneous domains with 16 and 28% volume proportion of high-

311

permeability facies. As a reminder, the ratio between porosities of the two facies is 2.6.

312

This means that grid-cells in the high-permeability facies have 2.6 times more pore volume

313

(PV) than the low-permeability facies grid-cells. Equations (1) and (2) suggest that lower

314

porosity facies have a lower λc and tc . Porosity affects the CO2 transport and the diffusive

315

flux (respectively, Equations (1) and (2) in Text S1), and its dispersion-width (Equation (3)).

316

As shown in Figure S3, the first fingers develop in high-permeability facies (see results of

317

1.5 and 2.5 years suggesting separated island of fingers through high-permeability facies).

318

However, after 4.5 years the lower permeability-porosity facies become saturated with CO2

14

ACS Paragon Plus Environment

Page 14 of 33

Page 15 of 33

Environmental Science & Technology

319

as well. This increases the density contrast and triggers instabilities. As a result, while the

320

advective Darcy flux (∝ k) is lower in low-permeability facies, the finger growth rates (∝ k/φ)

321

can be comparable to that in the higher permeability regions for both cases with 16 and 28%

322

volume proportions of high-permeability facies. This effect becomes more pronounced as the

323

volume proportion of low-permeability facies increases.

324

Figure 2b compares the temporal evolution of σz for different volume proportions of

325

high-permeability facies. The heterogeneous domains with fully connected clusters (28 and

326

40% cases with high connectivity of 27 and 39.6%, respectively) show a very similar evo-

327

lution of σz . The evolution for the 5 and 16% cases exhibit a greater difference because

328

the gravitational flow dynamics are more affected by the spatial organizations of both low-

329

and high-permeability facies. An interesting and novel finding is that the heterogeneous

330

simulations show fairly similar scaling relations as compared to the four flow regimes in

331

their homogeneous counterparts. Specifically, we observe both the first and third diffusion-

332

dominated flow regimes, and that these regimes last the longest for the 5 and 16% cases. The

333

similarities in behavior are clear in Figure 2d, which offers a direct quantitative comparison

334

between the homogeneous and heterogeneous simulations.

335

It is clear that the temporal evolution and the scaling of the third and the fourth regimes

336

are the same for homogeneous and heterogeneous simulations. In the first two flow regimes

337

deviations from the homogeneous case are related to the tortuous pathways arising from

338

the facies architecture. The former has significant practical implications in determining

339

the long-term evolution of CO2 convective mixing. Regardless of the complex fingering

340

and mixing dynamics at early times, the long-term storage permanence is determined by the

341

convective mixing in the fourth flow regime. Based on our findings, if the volume proportions

342

of facies types can be determined and their properties measured (e.g., from well-log data

343

or outcrops), then suitably averaged petrophysical properties (permeability and porosity)

344

should be sufficient to predict the long-term fate and storage performance of CO2 .

345

More generally, consider the variables that determine the mixing dynamics in Equa-

15

ACS Paragon Plus Environment

Environmental Science & Technology

346

tion (1) (and similarly in the Rayleigh number). Because of the low solubility of CO2 in

347

water (∼ 1–3 mol%), the maximum density difference ∆ρmax is only around 1% with modest

348

variation. Brine viscosity depends on temperature, but not significantly on CO2 composi-

349

tion. The diffusion coefficient for CO2 in brine depends on temperature and pressure but

350

also likely varies by less than an order of magnitude. Permeability can vary by many orders

351

of magnitude across different types of formations. However, 1) high permeability aquifers

352

are most favorable for CO2 storage, because they have the highest injectivity, and 2) if the

353

porosity is correlated with permeability at all, the combination φ/k in Equation (1) shows

354

less variation. For these reasons, the on-set time of instability, typical size (wavelength) of

355

fingers, and propagation speed may not vary as much as a dimensionless model may suggest

356

in terms of Rayleigh numbers that vary by many orders of magnitude. Typical Rayleigh

357

numbers for real CO2 storage formations may be more constrained.

358

Below we perform additional sensitivity analyses to better understand the potential effects

359

of heterogeneity and connectivity on the density-driven flow of CO2 . We perform simulations

360

1) with constant average porosities, 2) with within-cell scale anisotropy, 3) with shorter

361

correlation length for facies types, and 4) with permeability variation within each facies.

362

Homogeneous porosity

363

To separate the effects of heterogeneous porosities versus heterogeneous permeabilities, we

364

consider the same heterogeneous permeability distributions as before, but with a homoge-

365

neous (mean, as in Table S1) porosity.

366

Figure 1b shows the mole fraction of CO2 after 1.5, 2.5, 4.5, 8.5, 11, and 13.5 years

367

for a binary domain with 28% high-permeability facies. Media with 1) homogeneous k and

368

homogeneous φ (left column), 2) heterogeneous k and homogeneous φ (middle column), and

369

3) heterogeneous k and heterogeneous φ (right column) are shown. As a linear combination

370

of the porosities of the high- and low-permeability facies (0.34 and 0.13, respectively), the

371

mean porosity values always fall between the two (see Table S1). Therefore a k/φ ratio for

16

ACS Paragon Plus Environment

Page 16 of 33

Page 17 of 33

Environmental Science & Technology

ACS Paragon Plus Environment

Environmental Science & Technology

381

and φ as shown in the left panels. In high-permeability facies, tc and λc are shorter and

382

fingering is more pronounced. In general tc is smaller for all heterogeneous media (regardless

383

of porosity spatial variability) compared to homogeneous ones (Figure 2d).

384

The σz scaling is different for simulations with heterogeneous k but homogeneous φ

385

(Figure 2c). Similar to other cases, σz grows diffusively with σz ∝ t0.5 . The other three

386

distinct flow regimes observed in homogeneous and fully heterogeneous domains are smoothed

387

out in heterogeneous formations with uniform φ. After fingers are triggered the σz grows

388

faster in fully connected clusters, leading to a higher σz for both volume proportions of 28

389

and 40% with fully connected clusters. The volume proportion of 5% shows smaller σz .

390

For volume proportions with fully connected clusters, we see predominantly ballistic flow

391

regime (σz ∝ t1 ). The σz reaches the maximum σz at comparable times as the equivalent

392

homogeneous cases. For 5 and 16% volume proportions (below the percolation threshold), we

393

observe the ballistic regime only in part of the σz scaling. For example, volume proportion of

394

16% has a drop in σz after 10 years. The 5% case also shows some deviations from the ballistic

395

regime. This is related to the spatial organization and the location of high-permeability facies

396

and it is a realization-dependent phenomenon. The ensemble results of many realizations

397

should be considered (the number of realizations are not known a priori and could be as

398

many as hundreds) for cases that are below the percolation threshold (0–20%). However,

399

such analyses are computationally expensive and are outside the scope of this work.

400

Anisotropy ratio

401

As one other potentially important formation property, we consider small-scale (within grid

402

cells) anisotropy. Horizontal permeability, kh , is typically larger than the vertical kv . We

403

define the vertical to horizontal permeability as anisotropy ratio γ = kv /kh , which is typically

404

< 1. We create γ < 1 by increasing the horizontal permeability kh while keeping kv constant.

405

This way we are able to study anisotropy while comparing the results for the same buoyancy

406

scale (i.e., with equal driving force for vertical convective flux) but different γ. There are only

18

ACS Paragon Plus Environment

Page 18 of 33

Page 19 of 33

Environmental Science & Technology

407

a few works that have considered how the dissolution flux of CO2 is affected by anisotropic

408

porous media 68–70 for homogeneous or relatively simple (e.g., layered) heterogeneous media.

409

Figure 3 shows how γ < 1 affects the time evolution of σz in fully heterogeneous binary

410

media, with the 28% volume proportion of high-permeability facies.

411

We find that higher cell-scale anisotropy (i.e., lower γ) leads to earlier onset of convection

412

(i.e., tc ) in both homogeneous and heterogeneous media. However, this effect is more pro-

413

nounced for heterogeneous media. The degree of anisotropy considerably affects the onset

414

time and results in a shorter first diffusive regime. Following the instability onset, we find

415

that γ < 1 does not significantly affect the σz scaling in the early-time convective regime for

416

both homogeneous and heterogeneous media. Anisotropy does appear to lower the σz tem-

417

poral scaling for the last, fully convective, regime (again more pronounced for heterogeneous

418

media). The reason is that, although kv is the same for different γ, kh is higher for smaller

419

γ. This leads to more lateral mixing and diffusion, which lowers the CO2 concentration and

420

density within the sinking fingers. As a result, the (negative) buoyancy driving force and

421

the associated vertical spreading rate are weakened in the last (fourth) regime.

422

Correlation length and variability within facies

423

In the previous sections, we considered an isotropic integral scale (λ) of 10 m (X hereafter).

424

In this section we consider binary structures with 0.5 X and 0.1 X for a volume proportion

425

of 28%. Note that both k and φ are spatially heterogeneous. The results are shown in

426

Figure S4a, and suggest that changes in the correlation length have negligible effect on σz ,

427

and the scaling is similar for all cases.

428

Figure S4b shows σz (for a volume proportion of 28%) with different values of the per-

429

meability variance (σj2 ) within each facies. In previous analyses we did not incorporate any

430

variability in k since we assumed that the spatial organization of facies types is the major

431

source of heterogeneity, and that the spatial variabilities in k and φ within each facies are

432

of secondary importance. 72 Adding σj2 changes our binary systems to bimodal media. Fig-

19

ACS Paragon Plus Environment

Environmental Science & Technology

ACS Paragon Plus Environment

Page 20 of 33

Page 21 of 33

Environmental Science & Technology

441

respectively, 19 and 28% volume proportions of OFC cells are provided in Figure S5. We

442

also provide the results for homogeneous equivalent domains, i.e., with the same mean

443

permeability and porosity, in Figure S6. We note that a mean (scalar) permeability does not

444

capture effective anisotropy that may exist in the complex multiscale heterogeneity (defining

445

such anisotropy requires sophisticated upscaling techniques that are beyond the scope of this

446

work).

447

The difference in fingering behavior between a homogeneous (Figure S6) and a multi-

448

scale heterogeneous domain (Figure S5) appears to be more pronounced than for the binary

449

heterogeneity in terms of onset time and typical wavelength (with more and smaller fingers

450

in the homogeneous domain; see also Figure 4). This would not necessarily be surprising

451

given the wide range and non-trivial distribution of facies permeabilities. However, more

452

interestingly, the global measure of dispersion-width, σz , is surprisingly insensitive to the

453

degree and type of multiscale heterogeneity (Figure 4).

454

The inset in Figure 4 shows the same results but in a linear, rather than logarithmic,

455

scale. It is clear that cases with fully connected clusters (24 and 28%) exhibit identical time

456

evolution for σz . The architecture below the percolation threshold (the 19% case) has a

457

slightly lower growth rate, which reflects how the spatial distribution of disconnected flow

458

paths slows down the sinking of gravitational fingers to the bottom of domain.

459

Clearly, there are also some differences in the σz evolution between the heterogeneous

460

and the ‘equivalent’ homogeneous domains. The flow dynamics in the heterogeneous domain

461

represent essentially a superposition of the flow behavior in different clusters with varying

462

permeabilities. As such, the four flow regimes that are clearly visible in a homogeneous

463

medium are somewhat smeared out in the heterogeneous ones (though still identifiable).

464

A key take-away point is that despite the considerable complexity of the multiscale het-

465

erogeneous media, the scaling behavior for all heterogeneous and homogeneous domains is

466

remarkably similar. Most importantly, once the smallest early-time fingers have percolated

467

through potentially highly tortuous flow paths, temporal evolution of the σz always ends up

21

ACS Paragon Plus Environment

Environmental Science & Technology

ACS Paragon Plus Environment

Page 22 of 33

Page 23 of 33

Environmental Science & Technology

478

years, then it appears that a simple scaling relation is able to provide reasonable estimates,

479

regardless of the complexity of the formation heterogeneity.

480

Acknowledgement

481

The first author was supported by the U. S. Department of Energy (DOE) Office of Fossil

482

Energy funding to Oak Ridge National Laboratory (ORNL) under project FEAA-045. ORNL

483

is managed by UT-Battelle for the U.S. DOE under Contract DE-AC05-00OR22725. The

484

third and fifth authors was supported by funding as part of the Center for Geological Storage

485

of CO2 , an Energy Frontier Research Center funded by DOE, Office of Science, Basic Energy

486

Sciences under Award # DE-SC0C12504.

487

Supporting Information Available

488

The Supporting Information includes the following sections.

489

• Text S1: Governing equations

490

• Table S1: Mean homogeneous permeability and porosity in binary cases

491

• Figure S1: Heterogeneous binary models

492

• Figure S2: Permeability fields of multiscale heterogeneity model

493

• Figure S3: Visualizations of CO2 mole fraction at different simulation times for het-

494

495

496

497

498

erogeneous binary models • Figure S4: Sensitivity results of dispersion-width to correlation length and permeability variance within facies • Figure S5: Visualizations of CO2 mole fraction at different simulation times for multiscale heterogeneous models 23

ACS Paragon Plus Environment

Environmental Science & Technology

499

500

• Figure S6: Visualizations of CO2 mole fraction at different simulation times for homogeneous equivalents of multiscale models

501

• Text S2: References

502

This material is available free of charge via the Internet at http://pubs.acs.org/.

503

References

504

(1) Eccles, J. K.; Pratson, L.; Newell, R. G.; Jackson, R. B. Physical and economic potential

505

of geological CO2 storage in saline aquifers. Environ. Sci. Technol. 2009, 43, 1962–1969.

506

(2) Mathias, S. A.; Gluyas, J. G.; Goldthorpe, W. H.; Mackay, E. J. Impact of maximum

507

allowable cost on CO2 storage capacity in saline formations. Environ. Sci. Technol.

508

2015, 49, 13510–13518.

509

510

(3) Bourg, I. C.; Beckingham, L. E.; DePaolo, D. J. The nanoscale basis of CO2 trapping for geologic storage. Environ. Sci. Technol. 2015, 49, 10265–10284.

511

(4) Yang, C.; Treviño, R. H.; Zhang, T.; Romanak, K. D.; Wallace, K.; Lu, J.; Mickler, P. J.;

512

Hovorka, S. D. Regional assessment of CO2 –solubility trapping potential: A case study

513

of the coastal and offshore Texas Miocene Interval. Environ. Sci. Technol. 2014, 48,

514

8275–8282.

515

516

(5) Orr Jr, F. M. CO2 capture and storage: are we ready? Energy Environ. Sci. 2009, 2, 449–458.

517

(6) Soltanian, M. R.; Amooie, M. A.; Cole, D. R.; Graham, D. E.; Hosseini, S. A.; Hov-

518

orka, S.; Pfiffner, S. M.; Phelps, T. J.; Moortgat, J. Simulating the Cranfield geological

519

carbon sequestration project with high-resolution static models and an accurate equa-

520

tion of state. Int. J. Greenh. Gas Control. 2016, 54, 282–296.

24

ACS Paragon Plus Environment

Page 24 of 33

Page 25 of 33

521

522

Environmental Science & Technology

(7) Juanes, R.; Spiteri, E.; Orr, F.; Blunt, M. Impact of relative permeability hysteresis on geological CO2 storage. Water Resour. Res. 2006, 42, W12418.

523

(8) Gershenzon, N. I.; Ritzi, R. W.; Dominic, D. F.; Soltanian, M.; Mehnert, E.; Ok-

524

wen, R. T. Influence of small-scale fluvial architecture on CO2 trapping processes in

525

deep brine reservoirs. Water Resour. Res. 2015, 51, 8240–8256.

526

(9) Farajzadeh, R.; Salimi, H.; Zitha, P. L.; Bruining, H. Numerical simulation of density-

527

driven natural convection in porous media with application for CO2 injection projects.

528

Int. J. Heat Mass Transfer. 2007, 50, 5054–5064.

529

(10) Rapaka, S.; Chen, S.; Pawar, R. J.; Stauffer, P. H.; Zhang, D. Non-modal growth of

530

perturbations in density-driven convection in porous media. J. Fluid Mech. 2008, 609,

531

285–303.

532

(11) Farajzadeh, R.; Ranganathan, P.; Zitha, P. L. J.; Bruining, J. The effect of heterogeneity

533

on the character of density-driven natural convection of CO2 overlying a brine layer.

534

Adv. Water Resour. 2011, 34, 327–339.

535

536

(12) Hidalgo, J. J.; Dentz, M.; Cabeza, Y.; Carrera, J. Dissolution patterns and mixing dynamics in unstable reactive flow. Geophys. Res. Lett. 2015, 42, 6357–6364.

537

(13) Soltanian, M. R.; Amooie, M. A.; Dai, Z.; Cole, D.; Moortgat, J. Critical dynamics of

538

gravito-convective mixing in geological carbon sequestration. Sci. Rep. 2016, 6, 35921.

539

(14) Agartan, E.; Trevisan, L.; Cihan, A.; Birkholzer, J.; Zhou, Q.; Illangasekare, T. H. Ex-

540

perimental study on effects of geologic heterogeneity in enhancing dissolution trapping

541

of supercritical CO2 . Water Resour. Res. 2015, 51, 1635–1648.

542

543

(15) Daniel, D.; Riaz, A.; Tchelepi, H. A. Onset of natural convection in layered aquifers. J. Fluid Mech. 2015, 767, 763–781.

25

ACS Paragon Plus Environment

Environmental Science & Technology

544

545

(16) Riaz, A.; Hesse, M.; Tchelepi, H.; Orr, F. Onset of convection in a gravitationally unstable diffusive boundary layer in porous media. J. Fluid Mech. 2006, 548, 87–111.

546

(17) Hassanzadeh, H.; Pooladi-Darvish, M.; Keith, D. W. Scaling behavior of convective

547

mixing, with application to geological storage of CO2 . AlChE J. 2007, 53, 1121–1131.

548

(18) Pruess, K. Numerical modeling studies of the dissolution-diffusion-convection process

549

during CO2 storage in saline aquifers. LBNL 2008, 1243E

550

(19) Pau, G. S.; Bell, J. B.; Pruess, K.; Almgren, A. S.; Lijewski, M. J.; Zhang, K. High-

551

resolution simulation and characterization of density-driven flow in CO2 storage in

552

saline aquifers. Adv. Water Resour. 2010, 33, 443–455.

553

(20) Neufeld, J. A.; Hesse, M. A.; Riaz, A.; Hallworth, M. A.; Tchelepi, H. A.; Huppert, H. E.

554

Convective dissolution of carbon dioxide in saline aquifers. Geophys. Res. Lett. 2010,

555

37, L22404.

556

(21) Shahraeeni, E.; Moortgat, J.; Firoozabadi, A. High-resolution finite element methods

557

for 3D simulation of compositionally triggered instabilities in porous media. Comput.

558

Geosci. 2015, 19, 899–920.

559

(22) Soltanian, M. R.; Ritzi, R. W. A new method for analysis of variance of the hydraulic

560

and reactive attributes of aquifers as linked to hierarchical and multiscaled sedimentary

561

architecture. Water Resour. Res. 2014, 50, 9766–9776.

562

(23) Lu, J.; Cook, P.; Hosseini, S.; Yang, C.; Romanak, K.; Zhang, T.; Freifeld, B.;

563

Smyth, R.; Zeng, H.; Hovorka, S. Complex fluid flow revealed by monitoring CO2

564

injection in a fluvial formation. J. Geophys. Res. 2012, 117, B03208.

565

(24) Hosseini, S. A.; Lashgari, H.; Choi, J. W.; Nicot, J.-P.; Lu, J.; Hovorka, S. D. Static and

566

dynamic reservoir modeling for geological CO2 sequestration at Cranfield, Mississippi,

567

USA. Int. J. Greenh. Gas Control. 2013, 18, 449–462. 26

ACS Paragon Plus Environment

Page 26 of 33

Page 27 of 33

Environmental Science & Technology

568

(25) Ritzi, R. W.; Freiburg, J. T.; Webb, N. D. Understanding the (co)variance in petrophys-

569

ical properties of CO2 reservoirs comprising sedimentary architecture. Int. J. Greenh.

570

Gas Control. 2016, 51, 423–434.

571

(26) Gershenzon, N. I.; Soltanian, M. R.; Ritzi, R. W.; Dominic, D. F. Influence of small

572

scale heterogeneity on CO2 trapping processes in deep saline aquifers. Energy Procedia

573

2014, 59, 166–173.

574

(27) Krevor, S.; Pini, R.; Zuo, L.; Benson, S. M. Relative permeability and trapping of CO2

575

and water in sandstone rocks at reservoir conditions. Water Resour. Res, 2012, 48,

576

W02532.

577

(28) Dance, T.; Paterson, L. Observations of carbon dioxide saturation distribution and

578

residual trapping using core analysis and repeat pulsed-neutron logging at the CO2

579

CRC Otway site. Int. J. Greenh. Gas Control. 2016, 47, 210–220.

580

(29) Dai, Z.; Viswanathan, H.; Middleton, R.; Pan, F.; Ampomah,W.; Yang, C.; Jia, W.;

581

Xiao, T.; Lee, S.-Y.; McPherson, B. CO2 accounting and risk analysis for CO2 seques-

582

tration at enhanced oil recovery sites. Environ. Sci. Technol. 2016, 50, 7546–7554.

583

(30) Trevisan, L.; Pini, R.; Cihan, A.; Birkholzer, J. T.; Zhou, Q.; Illangasekare, T. H.

584

Experimental analysis of spatial correlation effects on capillary trapping of supercritical

585

CO2 at the intermediate laboratory scale in heterogeneous porous media. Water Resour.

586

Res. 2015, 51, 8791–8805.

587

(31) Trevisan, L.; Pini, R.; Cihan, A.; Birkholzer, J. T.; Zhou, Q.; González-Nicolás, A.;

588

Illangasekare, T. H. Imaging and quantification of spreading and trapping of carbon

589

dioxide in saline aquifers using meter-scale laboratory experiments. Water Resour. Res.

590

2016, 53, 485–502

591

(32) Trevisan, L.; Krishnamurthy, P.; Meckel, T. Impact of 3D capillary heterogeneity and

27

ACS Paragon Plus Environment

Environmental Science & Technology

592

bedform architecture at the sub-meter scale on CO2 saturation for buoyant flow in

593

clastic aquifers. Int. J. Greenh. Gas Control. 2017, 56, 237–249.

594

595

(33) Stauffer, D.; Aharony, A. Introduction to Percolation Theory; Taylor and Francis, Bristol, Pa., 1992.

596

(34) Sahimi, M. Applications of Percolation Theory; Taylor and Francis, Bristol, Pa., 1994.

597

(35) Guin, A.; Ritzi, R. W. Studying the effect of correlation and finite-domain size on

598

spatial continuity of permeable sediments. Geophys. Res. Lett. 2008, 35, L10402.

599

(36) Proce, C. J.; Ritzi, R. W.; Dominic, D. F.; Dai, Z. Modeling multiscale heterogeneity

600

and aquifer interconnectivity. Groundwater 2004, 42, 658–670.

601

(37) Ramanathan, R.; Guin, A.; Ritzi, R. W.; Dominic, D. F.; Freedman, V. L.;

602

Scheibe, T. D.; Lunt, I. A. Simulating the heterogeneity in braided channel belt de-

603

posits: 1. A geometric-based methodology and code. Water Resour. Res. 2010, 46,

604

W04516.

605

(38) Guin, A.; Ramanathan, R.; Ritzi, R. W.; Dominic, D. F.; Lunt, I. A.; Scheibe, T. D.;

606

Freedman, V. L. Simulating the heterogeneity in braided channel belt deposits: 2.

607

Examples of results and comparison to natural deposits. Water Resour. Res. 2010, 46,

608

W04515.

609

610

611

612

(39) Knudby, C.; Carrera, J. On the relationship between indicators of geostatistical, flow and transport connectivity. Adv. Water Resour. 2005, 28, 405–421. (40) Silliman, S. The influence of grid discretization on the percolation probability within discrete random fields. J. of Hydrol. 1990, 113, 177–191.

613

(41) Gershenzon, N.; Soltanian, M.; Ritzi, R. W.; Dominic, D. F. Understanding the impact

614

of open-framework conglomerates on water–oil displacements: the Victor interval of the

615

Ivishak Reservoir, Prudhoe Bay Field, Alaska. Pet. Geosci. 2015, 21, 43–54. 28

ACS Paragon Plus Environment

Page 28 of 33

Page 29 of 33

Environmental Science & Technology

616

(42) Massabó, M.; Bellin, A.; Valocchi, A. Spatial moments analysis of kinetically sorbing

617

solutes in aquifer with bimodal permeability distribution. Water Resour. Res. 2008,

618

44, W09424.

619

620

621

622

(43) Carle, S. F. T-PROGS: Transition probability geostatistical software. Univ. of Calif., Davis, Calif. 1999. (44) Bernabé, Y.; Mok, U.; Evans, B. in Thermo-Hydro-Mechanical Coupling in Fractured Rock, Birkhäuser Basel, 2003, 937–960.

623

(45) Deng, H.; Stauffer, P. H.; Dai, Z.; Jiao, Z.; Surdam, R. C. Simulation of industrial-scale

624

CO2 storage: Multi-scale heterogeneity and its impacts on storage capacity, injectivity

625

and leakage. Int. J. Greenh. Gas Control. 2012, 10, 397–418.

626

627

628

629

(46) Huang, L.; Ritzi, R. W.; Ramanathan, R. Conservative models: Parametric entropy vs. temporal entropy in outcomes. Groundwater 2012, 50, 199–206. (47) Lunt, I.; Bridge, J.; Tye, R. A quantitative, three-dimensional depositional model of gravelly braided rivers. Sedimentology 2004, 51, 377–414.

630

(48) Gershenzon, N. I.; Soltanian, M. R.; Ritzi, R. W.; Dominic, D. F.; Keefer, D.; Shaf-

631

fer, E.; Storsved, B. How does the connectivity of open-framework conglomerates within

632

multi-scale hierarchical fluvial architecture affect oil-sweep efficiency in waterflooding?

633

Geosp. 2015, 11, 2049–2066.

634

(49) Moortgat, J.; Sun, S.; Firoozabadi, A. Compositional modeling of three-phase flow

635

with gravity using higher-order finite element methods. Water Resour. Res. 2011, 47,

636

W05511.

637

638

(50) Moortgat, J.; Firoozabadi, A. Higher-order compositional modeling with Fickian diffusion in unstructured and anisotropic media. Adv. Water Resour. 2010, 33, 951–968.

29

ACS Paragon Plus Environment

Environmental Science & Technology

639

(51) Moortgat, J.; Li, Z.; Firoozabadi, A. Three-phase compositional modeling of CO2 in-

640

jection by higher-order finite element methods with CPA equation of state for aqueous

641

phase. Water Resour. Res. 2012, 48, W12511.

642

(52) Moortgat, J.; Firoozabadi, A. Mixed-hybrid and vertex-discontinuous-Galerkin finite el-

643

ement modeling of multiphase compositional flow on 3D unstructured grids. J. Comput.

644

Phys. 2016, 315, 476–500.

645

(53) Moortgat, J. B.; Firoozabadi, A.; Li, Z.; Espósito, R. CO2 injection in vertical and

646

horizontal cores: Measurements and numerical simulation. SPE J. 2013, 18, 331–344.

647

(54) Moortgat, J.; Amooie, M. A.; Soltanian, M. R. Implicit finite volume and discontinuous

648

Galerkin methods for multicomponent flow in unstructured 3D fractured porous media.

649

Adv. Water Resour. 2016, 96, 389–404.

650

651

(55) Amooie, M. A.; Soltanian, M. R.; Moortgat, J. Hydro-thermodynamic mixing of fluids across phases in porous media. Geophys. Res. Lett. 2017, 44,3624–3634.

652

(56) Amooie, M. A.; Soltanian, M. R.;Xiong, F.; Dai, Z.; Moortgat, J. Mixing and spread-

653

ing of multiphase fluids in heterogeneous bimodal porous media. Geomechanics and

654

Geophysics for Geo-Energy and Geo-Resources 2017, 1–20.

655

656

657

658

(57) Moortgat, J. Viscous and gravitational fingering in multiphase compositional and compressible flow. Adv. Water Resour. 2016, 89, 53–66 (58) Firoozabadi, A. Thermodynamics and Applications of Hydrocarbon Energy Production; McGraw-Hill Education, New York, 2016.

659

(59) Nasrabadi, H.; Moortgat, J.; Firoozabadi, A. New three-phase multicomponent com-

660

positional model for asphaltene precipitation during CO2 injection using CPA-EOS.

661

Energy Fuels 2016, 30, 3306–3319.

30

ACS Paragon Plus Environment

Page 30 of 33

Page 31 of 33

662

663

664

665

666

667

Environmental Science & Technology

(60) Li, Z.; Firoozabadi, A. Cubic-plus-association equation of state for water-containing mixtures: Is cross association necessary? AIChE J. 2009, 55, 1803–1813. (61) Li, Z.; Firoozabadi, A. General strategy for stability testing and phase-split calculation in two and three phases. SPE J. 2012, 17, 1–096. (62) Hidalgo, J. J.; Carrera, J.; Medina, A. Role of salt sources in density-dependent flow. Water Resour. Res. 2009, 45, W05503.

668

(63) Emami-Meybodi, H.; Hassanzadeh, H.; Green, C. P.; Ennis-King, J. Convective disso-

669

lution of CO2 in saline aquifers: Progress in modeling and experiments. Int. J. Greenh.

670

Gas Control. 2015, 40, 238–266.

671

672

673

674

675

676

(64) Martinez, M. J.; Hesse, M. A. Two-phase convective CO2 dissolution in saline aquifers. Water Resour. Res. 2016, 52, 585–599. (65) Aris, R. On the dispersion of a solute in a fluid flowing through a tube. Proc. R. Soc. Lond. A 1956, 235, 67–77. (66) Le Borgne, T.; Dentz, M.; Carrera, J. Lagrangian statistical model for transport in highly heterogeneous velocity fields. Phys. Rev. Lett. 2008, 101, 090601.

677

(67) Buscheck, T. A.; White, J. A.; Carroll, S. A.; Bielicki, J. M.; Aines, R. D. Managing

678

geologic CO2 storage with pre-injection brine production: a strategy evaluated with a

679

model of CO2 injection at Snøhvit. Energy Environ. Sci. 2016, 9, 1504–1512.

680

(68) Cheng, P.; Bestehorn, M.; Firoozabadi, A.Effect of permeability anisotropy on

681

buoyancy-driven flow for CO2 sequestration in saline aquifers.Water Resour. Res. 2012,

682

48, .

683

684

(69) Green, C.; Ennis-King, J. Steady dissolution rate due to convective mixing in anisotropic porous mediaAdv. Water Resour. 2014, 73, 65–73.

31

ACS Paragon Plus Environment

Environmental Science & Technology

685

686

687

688

689

690

(70) De Paoli, M.; Zonta, F.; Soldati, A. Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown. Physics of Fluids 2017, 29, . (71) Divins, D. Total Sediment Thickness of the World’s Oceans & Marginal Seas; NOAA National Geophysical Data Center, Boulder, CO, 2003 (72) Desbarats, A. J. Numerical estimation of effective permeability in sand-shale formations. Water Resour. Res. 1987, 23, 273–286.

32

ACS Paragon Plus Environment

Page 32 of 33

Page 33 of 33

Environmental Science & Technology

ACS Paragon Plus Environment