Combined inter- and intra-fracture free convection in fracture networks embedded in a low-permeability matrix

Accepted Manuscript

Combined inter- and intra-fracture free convection in fracture networks embedded in a low-permeability matrix Katharina Vujevi´c, Thomas Graf PII: DOI: Reference:

S0309-1708(15)00161-X 10.1016/j.advwatres.2015.07.014 ADWR 2426

To appear in:

Advances in Water Resources

Received date: Revised date: Accepted date:

28 October 2014 22 June 2015 23 July 2015

Please cite this article as: Katharina Vujevi´c, Thomas Graf, Combined inter- and intra-fracture free convection in fracture networks embedded in a low-permeability matrix, Advances in Water Resources (2015), doi: 10.1016/j.advwatres.2015.07.014

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • Inter- and intrafracture convection are simulated separately and in com-

CR IP T

bination • A critical Rayleigh number for interfracture convection is proposed • Interfracture convection dominates in continuous fracture circuits

• Combined inter- and intrafracture convection develops in complex frac-

AN US

ture networks

• Complex fracture networks require 3D simulation to capture free con-

AC

CE

PT

ED

M

vection

1

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Graphical Abstract

2

ACCEPTED MANUSCRIPT

CR IP T

Combined inter- and intra-fracture free convection in fracture networks embedded in a low-permeability matrix Katharina Vujevi´ca,∗, Thomas Grafa

Institute of Fluid Mechanics and Environmental Physics in Civil Engineering, Leibniz Universit¨ at Hannover. Appelstrasse 9A, 30167 Hannover, Germany

AN US

a

Abstract

M

Density-driven, haline free convection in fractured porous rock is studied in a 3-dimensional numerical model to evaluate the onset and interdependence of

ED

intra- and interfracture convection. When the rock matrix allows for solute diffusion, the most likely mode of convection in a regular fracture circuit is interfracture convection along the circuit such that regular fracture circuits can

PT

be modeled in a reduced dimensionality (2D) model. The critical Rayleigh number for the onset of interfracture convection in a regular fracture circuit

CE

is determined to be half of the critical Rayleigh number for convection in a single vertical fracture. In more complex fracture networks, vertical “dead-

AC

end fractures” and additional fracture intersections complicate convective patterns. Combined inter- and intrafracture convection is possible, such that ∗

Corresponding author Email addresses: [email protected] (Katharina Vujevi´c), [email protected] (Thomas Graf)

Preprint submitted to Advances in Water Resources

July 28, 2015

ACCEPTED MANUSCRIPT

convection is hard to predict.

1

CR IP T

Keywords: Groundwater; fractures; Rayleigh number; density-driven flow; 3-dimensional, HydroGeoSphere; numerical modelling 1. Introduction

In many deep groundwater aquifers, density-driven free convection is an

3

important process for solute transport [11, 36]. Free convection has to be con-

4

sidered when solute or heat influence groundwater density [31, 54]. When

5

dense brine is overlying less dense fresh water, this potentially unstable fluid

6

stratification may be dissipated by diffusion, or may lead to free convection in

7

case of dominating buoyancy forces [11, 24, 31, 54]. A number of authors em-

8

phasize the widespread importance of density-driven flow, as density-driven

9

flow transports solutes faster and farther into the aquifer then diffusion alone

10

[33, 37, 47]. Therefore, neglecting free convection will lead to a considerable

11

underestimation of solute and/or heat transport.

M

AN US

2

In low-permeability geological formations, such as shales or granites, frac-

13

tures may enable free convective flow even if the matrix is virtually imperme-

14

able to flow. Therefore, fracture networks are conjointly considered to play

15

an important role in solute and heat transport in low-permeability strata

16

and in determining the onset and patterns of free convection in such sys-

17

tems [19, 20, 26, 33, 34, 37, 43, 44]. Understanding the onset conditions of

18

free convection in fractured rock is increasingly important because fractured

AC

CE

PT

ED

12

19

rocks are widespread and gain importance as potential hosts for hazardous

20

waste disposal, heat and CO2 storage, natural gas extraction by fracking,

21

and are used for thermal energy exploitation. Furthermore, density-driven

22

flow in fractures and fault systems is known to play an important role in the 4

ACCEPTED MANUSCRIPT

23

formation of ore deposits [1, 8, 16, 49, 53, 55]. Convection in fractured low-permeability strata is possible along a circuit

25

of connected fractures (interfracture convection, “mode 1” as shown in Fig.

26

1a) or within a single fracture (intrafracture convection, “mode 2” as shown

27

in Fig. 1b) [37, 43]. Both types of convection have been examined separately

28

in a couple of theoretical studies with different considerations of fracture-

29

matrix interactions as discussed below.

b) isolated vertical fractures Intrafracture convection (mode 2B) (mode 2A)

AN US

a) 2-dimensional regular fracture circuit Interfracture convection (mode 1)

L =6 m

m

H=10 m

h=6 m

H=10 m

2B=6 m

h=6 m

L=0.2

CR IP T

24

z

y

M

x

W=10 m

W=10 m

ED

Figure 1: Conceptual models to study (a) interfracture convection and (b) intrafracture convection. Mode 2A was not considered in this study and is just shown for completeness.

Early studies on (thermal) convection within a single fracture approxi-

31

mated the fracture by a 3-dimensional fluid-filled vertical slab or by a sat-

32

urated porous box heated from below. Fracture walls were assumed to be

33

impermeable to flow and non-conducting, such that any fracture-matrix in-

34

teraction was neglected [7, 56]. In these cases, the onset of convection can be

35

predicted using the classical Rayleigh number Ra which compares buoyancy

AC

CE

PT

30

36

to diffusivity. If the critical Rayleigh number of Rac = 4π 2 is exceeded, con-

37

vection with streamlines parallel to the fracture plane (mode 2B) is expected

38

to establish [3, 9, 15, 30, 56].

39

Some studies accounted for heat-conducting fracture walls and applied a 5

ACCEPTED MANUSCRIPT

constant vertical temperature gradient to the lateral fracture walls [23, 25, 50,

41

52]. This boundary condition corresponds to assuming perfectly conducting

42

fracture wall and has a stabilizing effect on thermal convection [27]. Thus,

43

higher vertical temperature differences are required to initiate free convection

44

and the critical Rayleigh number greatly exceeds 4π 2 [8, 23]. However, heat

45

transport within the matrix was not considered [23, 25, 50, 52]. In reality,

46

this would correspond to a case where the fracture gains or looses heat along

47

its vertical walls, but the adjacent matrix does not change temperature.

CR IP T

40

More advanced models for convection within a fracture embedded in a

49

conducting matrix were not restricted to simulating processes within the

50

fracture, but also accounted for heat transport in the adjacent matrix [1, 21,

51

26, 27, 28, 44, 53]. The stabilizing effect of the matrix is less distinct in

52

this case, such that resulting critical Rayleigh numbers exceed 4π 2 , but are

53

below those determined with vertical fracture walls at constant temperature

54

gradient [27].

M

AN US

48

Interfracture convection within a fracture network was usually studied in

56

2-dimensional domains with fractures represented as 1-dimensional line ele-

57

ments such that convection within a fracture was neglected. Even without

58

considering intrafracture convection (mode 2), it was shown that fracture net-

59

works considerably affect the onset, strength and patterns of free convection

60

in fractured porous media [20, 34, 38, 43, 48]. Depending on the geometry

61

of individual fracture networks, convection was enhanced or inhibited, such

AC

CE

PT

ED

55

62

that making predictions is a complex task.

63

Very few studies link inter- and intrafracture convection. A theoretical

64

comparison of the likelihood of each mode of free (haline) convection is pro-

65

vide by Simmons et al. [37], who considered a fracture circuit embedded 6

ACCEPTED MANUSCRIPT

in an impermeable rock matrix. A comparison of critical concentration dif-

67

ferences (derived from critical Rayleigh numbers) led to the conclusion that

68

intrafracture convection parallel to the fracture plane (mode 2B) is the most

69

likely mode of convection and occurs for even very small density differences.

70

According to Simmons et al. [37], interfracture convection around an im-

71

permeable rock matrix block (mode 1) is likely as well, while intrafracture

72

convection normal to the fracture plane (mode 2A) requires large density dif-

73

ferences in combination with large fracture apertures and is hence unlikely.

CR IP T

66

Ramazanov [32] analytically determined the onset of interfracture ther-

75

mal convection (called transverse convection) along a fracture circuit and

76

intrafracture convection (called longitudinal convection) within the vertical

77

fractures of a fracture circuit embedded in an impermeable, conductive rock

78

matrix. Ramazanov’s [32] results suggest that the onset of convection in a

79

fracture circuit depends on fracture aperture and aspect ratio of the fracture

80

circuit. Interfracture convection sets in at lower Rayleigh numbers than in-

81

trafracture convection for small fracture apertures and fracture circuits where

82

the aspect ratio (height to width) is between 0.25 and 8.67.

83

Contribution of this study

PT

ED

M

AN US

74

The onset of possible free convection in fractures has been theoretically

85

studied and numerically simulated separately, either within a single fracture

86

(e.g. [26, 44]) or along a fracture network (e.g. [43, 48]). One of the reasons

87

for the lack of studies on combined intra- and interfracture convection is that

88

a combined consideration requires a fully 3-dimensional approach, which is

89

demanding especially for problems with complex geometries. Studies that

90

theoretically considered both inter- and intrafracture convection simplified

91

the fracture-matrix interaction by either averaging over fracture and

AC

CE

84

7

ACCEPTED MANUSCRIPT

matrix properties (Simmons et al. [37]) or assuming diffusion in the

93

fracture and the matrix to be identical (Ramazanov [32]). Thus, the

94

resulting probabilities for each mode of convection may not be applicable to a

95

discrete fracture circuit embedded in a porous matrix that allows for diffusive

96

transport and that has different diffusivity parameters than the fracture.

CR IP T

92

The purpose of this study is to examine intra- and interfracture haline

98

convection separately (Fig. 1) and in combination for simple fracture network

99

geometries shown in Fig. 2 and to determine under which conditions which

100

mode dominates. The effect of the adjacent matrix will be accounted for by

101

allowing for diffusion within the rock matrix, within the fractures, and across

102

the fracture-matrix interface. Thus, we want to contribute to a profound un-

103

derstanding of the interdependence of intra- and interfracture convection,

104

the onset of combined convection as well as the convective patterns which

105

can be expected in fractured low-permeability media. We further aim at

106

encompassing in which type of fracture networks a fully 3-dimensional repre-

107

sentation is required or when a reduced dimensionality representation for the

108

simulation of free convection yields reasonable results. We focus our study

109

on haline free convection on the basis of findings from Graf and Therrien

110

[21], Sharp et al. [33] and Alt-Epping and Zhao [1], who stated

111

that in low-permeability shale layers and geological fault zones,

112

salinity gradients are the dominating driving force for free convec-

113

tion (in contrast to temperature gradients), because geothermal

AC

CE

PT

ED

M

AN US

97

114

convection requires much higher fracture and matrix permeabili-

115

ties to establish than haline convection due to the large stabilizing

116

thermal diffusivity. This is the first numerical study of free convection in

117

fractured-porous rock that considers a combination of intra- and interfracture 8

ACCEPTED MANUSCRIPT

convection and allows for matrix diffusion. a) 3-dimensional regular fracture circuit Intra-, interfracture or combined convection (mode 1 or mode 2B or mode 12B)

H=10 m

h=6 m

H=10 m

2B=6 m

h=4 m

L=4 m

L=6 m

CR IP T

b) 3-dimensional complex fracture circuit Intra-, interfracture or combined convection (mode 1 or mode 2B or mode 12B)

circuit height

118

W=10 m

AN US

W=10 m

AC

CE

PT

ED

M

Figure 2: Conceptual models to study combined inter- and intrafracture convection. a) regular fracture circuit consisting of 4 fractures. Arrows show convection patterns of combined inter- and intrafracture convection (mode 12B) b) more complex fracture circuit consisting of 6 fractures.

9

ACCEPTED MANUSCRIPT

119

2. Methodology

120

2.1. Mathematical model In density-driven flow and transport problems, flow and trans-

122

port are coupled by fluid density which varies with concentration

123

and affects flow by adding a buoyancy component to the flow equa-

124

tion. In fully saturated porous media, the coupled governing equa-

125

tions are the Darcy equation and mass conservation equations for

126

fluid and solute mass. The non-linear system of equations is com-

127

pleted by the constitutive equations that define porous medium

128

and fracture freshwater hydraulic conductivity and relationships

129

between fluid density and salinity. The governing and constitutive

130

equations can be found in the Appendix.

AN US

CR IP T

121

Free convection in fractured porous rock is investigated numerically using

132

the variable-density groundwater flow and transport model HydroGeoSphere

133

[39] which applies the control volume finite element method to the flow equa-

134

tion and the Galerkin finite element method to the transport equation. The

135

primary unknown for flow is equivalent fresh water head h0 [L]

136

defined as [14]:

ED

PT

CE

h0 =

p +z ρ0 g

(1)

where p [ML−1 T−2 ] is fluid pressure, ρ0 [ML−3 ] is reference fluid density,

AC

137

M

131

138

g [LT−2 ] is gravitational acceleration and z [L] is elevation. The primary

139

unknown for transport is relative concentration c [-].

140

Fractures are represented as discrete 2-dimensional planes that

141

share common nodes with 3-dimensional elements of the rock ma10

ACCEPTED MANUSCRIPT

trix. The common node approach ensures continuity of the un-

143

knowns at the fracture-matrix interface. By superimposing the

144

contributions from fracture and matrix at common nodes, fluid and

145

solute mass exchange do not have to be explicitly calculated in the

146

locally mass conservative discretized flow and transport equations.

147

A linear relationship between fluid density and relative concentra-

148

tion is assumed as described in the Appendix. Further details on

149

the HydroGeoSphere model can be found in the Appendix, in Graf

150

and Therrien [20], and in Therrien et al. [39].

151

2.2. Dimensionless numbers

AN US

CR IP T

142

The Rayleigh number is used to predict the onset of convection [35, 42].

153

The Sherwood number is used to characterize the strength of free convection

154

[4, 43].

155

Rayleigh Number

ED

M

152

The Rayleigh number Ra [-] relates buoyancy that promotes free con-

157

vection to diffusion that dissipates free convection [11, 35, 42, 43, 45]. The

158

classical Rayleigh number Ra for salinity-driven free convection in homoge-

159

neous porous media reads:

CE

H 2 kgρ0 β ∆c H Ra = Dm µ

(2)

where H [L] is the height of the domain, k [L2 ] is the permeability of

AC

160

PT

156

161

the porous medium, g [L T−2 ] is gravitational acceleration, µ [M L−1 T−1 ] is

162

fluid dynamic viscosity, and Dm [L2 T−1 ] is the effective diffusion coefficient

163

in the porous matrix. If Ra is smaller than the critical Rayleigh number

164

Rac , any small initial perturbation (trigger) decays with time, such that the 11

ACCEPTED MANUSCRIPT

system becomes purely diffusive and stable [15, 31, 54]. If Ra is larger than

166

Rac , the system is dominated by buoyancy-induced free convection, and it is

167

called unstable [15, 31, 54]. The value of the critical Rayleigh number Rac

168

depends on the applied flow and transport boundary conditions [30]. The

169

critical Rayleigh number is Rac = 4π 2 for a closed box filled with saturated

170

homogeneous porous media, with fixed concentrations cmax , cmin at the top

171

and bottom of the domain, respectively.

CR IP T

165

When a single 2-dimensional open vertical fracture is considered instead

173

of a box of porous medium, the permeability k has to be replaced by kfr , the

174

diffusion coefficient in the fracture corresponds to D0 in Eq. (22), and the

175

height of the fracture h [L] is introduced, such that Rafr reads [44]:

AN US

172

176

(3)

M

h2 kfr gρ0 β ∆c H Rafr = D0 µ

Using linear stability analysis, Malkovsky and Pek [26] and Wang et al.

[44] determined a critical Rayleigh number Rac-fr

178

for convection in a single vertical fracture embedded in a fully

179

3-dimensional fracture-matrix system. The matrix was assumed to

180

be impermeable, but heat-conducting, and upper and lower bound-

181

aries of the setup were impermeable and of constant temperature.

182

Malkovsky and Pek [26] and Wang et al. [44] pointed out that

183

a conducting matrix has a stabilizing effect, due to heat transfer

184

along the fracture-matrix interface. The magnitude of heat ex-

185

change depends on the ratio of fracture aperture (2b) to fracture

186

height h, and on the thermal conductivity ratio  between matrix

187

and fracture [44]. Consequentially, the critical Rayleigh number

188

Rac-fr can be formulate as a function of those two ratios, where

AC

CE

PT

ED

177

12

ACCEPTED MANUSCRIPT

Malkovsky and Pek [26] determined Rac-fr = 16.38·h/(2b) and Wang

190

et al. [44] derived Rac-fr = 16.32 · h/(2b). The formulation of the critical

191

Rayleigh number determined by Wang et al. [44] is adapted here for haline

192

convection by replacing the thermal conductivity ratio  between matrix and

193

fracture by the solutal diffusivity ratio rD = Dm /D0 : Rac-fr = 16.32rD

h (2b)

CR IP T

189

(4)

We combine Eq. (3) with Eq. (4) and replace kfr by (2b)2 /12 [5, 20] to

195

obtain a critical fracture aperture (2b)c in Eq. (5), which is the parameter

196

that is usually unknown. With parameters given in Table 1 that are typical

197

for limestone, claystone or shale [6, 29], (2b)c becomes:

3

rD µ ∆c β H hρg

· 12 · 16.32 = 1.74 · 10−5 m

M

(2b)c =

r

AN US

194

(5)

A critical Rayleigh number Rac-circ for convection along a fracture circuit

199

has been formulated by Ramazanov [32]. As that formulation is quite com-

200

plex and does not account for diffusivity differences between fracture and

201

matrix, it is not shown here. A new formulation is proposed in this paper in

202

section 3.1.

203

Sherwood Number

PT

The dimensionless Sherwood number Sh is used as an indicator for the

strength of convection [4, 43]:

AC

205

CE

204

ED

198

206

Sh =

j Dm ρ0 ∆c H

(6)

where j [M T−1 L−2 ] is the total mass flux through the top layer of the 13

ACCEPTED MANUSCRIPT

Table 1: Simulation parameters that are typical for a claystone or shale

9.81 m s−2 0.1 m−1 0.07 1000 kg m−3 1.1×10−3 kg m−1 s−1 1.0×10−18 m2 0.1 0.1 10−9 m2 s−1 10−11 m2 s−1 0.01 0.1 m 0.005 m

CR IP T

g ∆c/H β ρ0 µ km θ τ D0 Dm rD αL αT

AN US

Gravitational acceleration Concentration gradient Expansion coefficient Reference fluid density (fresh water) Fluid dynamic viscosity Rock matrix permeability Rock matrix porosity Rock matrix tortuosity Free solution diffusion coefficient Diffusion coefficient in the matrix (D0 θτ ) Diffusivity ratio matrix to fracture Longitudinal dispersivity Transverse dispersivity

simulation domain, and the denominator describes the expected mass flux

208

through this layer for a purely diffusive regime. Therefore, Sh is 1 in the

209

steady state of stable and purely diffusive scenarios, and Sh is greater than

210

1 once free convection begins [17].

211

2.3. Conceptual models of simulation scenarios

ED

M

207

Figs. 1 and 2 depict the basic conceptual models we examined, where Fig.

213

2 shows the scenarios that are novel. We assume a 10 m×6 m×10 m (width×

214

length× height) matrix block, in which fractures of height h = 6 m (vertical

215

fractures) and width w = 6 m (horizontal fractures) are embedded. Thus,

216

the fractures do not traverse the entire height nor width of the block, but

217

they traverse the entire length L of the block in y-direction. The matrix has

AC

CE

PT

212

218

a very low hydraulic permeability (10−18 m2 ), but allows for diffusive mass

219

transport. As a reference, we separately consider each mode of convection:

220

interfracture convection mode 1 (Fig. 1a) and intrafracture convection mode

221

2B (Fig. 1b) as described below. Fig. 1a represents a 2-dimensional cross

14

ACCEPTED MANUSCRIPT

section of the fracture network, which is a common approach used to simulate

223

free convection in fractured rock: The continuous fracture circuit represents

224

the most fundamental building block (as of [37]) of an orthogonal fracture

225

network. Fig. 1b refers to free convection in an isolated vertical fracture,

226

e.g. in fractured rock where almost uniform fracture orientations avoid

227

inter-connection of fractures. As a novelty, we then consider a 3-dimensional

228

fracture circuit (Fig. 2a), where both modes of convection are theoretically

229

possible. Finally, we show how our results can be applied to other fracture

230

network geometries and discuss the predictability of convection in a more

231

complex 3-dimensional fracture network (Fig. 2b). Thus, we approach the

232

reality of fractured rocks with a variety of geometries of connected fracture

233

clusters.

AN US

CR IP T

222

For each geometry (Fig. 1a, Fig. 1b, Fig. 2a, Fig. 2b), we modify

235

the fracture aperture (2b) to numerically determine the critical (2b) for the

236

onset of convection. The concentration gradient along the vertical axis and

237

all other parameters will be kept constant as listed in Table 1 unless other-

238

wise stated. The onset of convection is determined calculating Sh for each

239

scenario, were Sh = 1 indicates the absence of convection, and Sh > 1 indi-

240

cates the existence of free convection. We determine the convective patterns

241

visually by evaluating the shape of the streamlines and isochlors.

ED

PT

It should be made clear that we distinguish between the geometry of a

fracture network as shown in Figs. 1 and 2 and the modes of convection

AC

243

CE

242

M

234

244

defined by the occurring flow patterns. The 2-dimensional fracture circuit

245

(Fig. 1a) only allows for interfracture convection (named mode 1). Isolated

246

vertical fractures (Fig. 1b) only allow for intrafracture convection (named

247

mode 2B). Only the 3-dimensional fracture circuits (Figs. 2a, b) allow for 15

ACCEPTED MANUSCRIPT

248

either interfracture convection (mode 1) or intrafracture convection (mode

249

2B) or for combined convection (named mode 12B). It should be acknowledged that we do not study intrafracture convection

251

perpendicular to the fracture plane (mode 2A as of Simmons et al. [37] or

252

2-D slender circle flow as of Zhao et al. [50]). According to Simmons

253

et al. [37] and several earlier studies [3, 9, 15, 56], convection parallel to the

254

fracture plane is the preferred mode within a finite vertical fracture, such

255

that neglecting intrafracture convection perpendicular to the fracture plane

256

is acceptable without loss of generality of this study. In some studies where

257

conduction/diffusion across the fracture walls was allowed [1, 26, 27, 44, 50],

258

intrafracture convection cells parallel to the fracture plane are described as

259

(weakly) 3-dimensional: Assuming a linear vertical temperature dis-

260

tribution along the impermeable vertical fracture walls, Kassoy and

261

Cotte [23] and Zhao et al. [50, 54] theoretically showed that con-

262

vection cells with rotation axes parallel to the fault width (fracture

263

aperture) dominate thermal convective patterns. The convective

264

pattern was described as 3-dimensional (called “standard convec-

265

tive mode” or “finger-like convective mode”) because the fracture

266

walls have a dampening effect on the flow field, such that flow ve-

267

locities are highest in the central vertical plane of the fault. In

268

contrast to the studies of Kassoy and Cotte [23] and Zhao et al.

269

[50, 54], concentration at vertical fracture walls in our study is not

AC

CE

PT

ED

M

AN US

CR IP T

250

270

fixed, but transport in the rock matrix surrounding the vertical

271

fracture is fully considered, similar to the theoretical analyses by

272

Wang et al. [44] and Malkovsky and Pek [26, 27]. Wang et al.

273

[44] and Malkovsky and Pek [26, 27] showed that convection ap16

ACCEPTED MANUSCRIPT

proaches a 2-dimensional pattern for low ratios of matrix diffusiv-

275

ity to fracture diffusivity, and for low ratios of fracture aperture to

276

fracture height, which are both given in our study. Thus, averaging

277

over the fracture width seems a valid approach here, even if that

278

2-dimensional approach does not resolve the flow field within the

279

fracture in the third dimension. The main advantage of the present

280

2-dimensional fracture approximation is that computational efforts

281

are significantly reduced. As the first step of dealing with compli-

282

cated geoscience problems, such an approximation is commonly

283

acceptable for numerical models [55]. Thus, in the present study,

284

the term “intrafracture convection” corresponds to convective flow

285

with streamlines parallel to the fracture plane (mode 2B).

286

Boundary and initial conditions

M

AN US

CR IP T

274

All domain boundaries are impermeable to flow. The initial state is in hy-

288

drostatic equilibrium. We assign constant concentration boundary conditions

289

to the bottom and top of the domain, while lateral boundaries are assigned

290

a zero-mass flux boundary condition. Concentration at the top is cmax = 1,

291

and bottom concentration is cmin = 0. Initially, solute is distributed linearly

292

between upper and lower boundary. A physical trigger was applied, where

293

the central node of the vertical fracture was assigned a 10% higher concentra-

294

tion relative to the ambient fluid. In linear stability analyses, a system

295

is defined as stable, when a small initial perturbation (trigger) de-

296

cays with time [15, 31]. Using numerical simulations to determine

297

the onset and mode of convection limits the number of possible

298

initial perturbations that may be tested. We varied the magnitude

299

of the applied trigger as well as the location of the trigger with-

AC

CE

PT

ED

287

17

ACCEPTED MANUSCRIPT

out observing any change in the critical Rayleigh number and the

301

patterns of free convection. The fact that results for the onset of

302

convection in a single vertical fracture are in good agreement with

303

the results of linear stability analyses obtained by Wang et al. [44]

304

and Malkovsky and Pek [26] further substantiates the reliability of

305

our numerical simulations. The advantage of using numerical sim-

306

ulations to determine onset and patterns of free convection is that

307

this approach allows to test more complex scenarios than linear

308

stability analysis.

309

Spatial and temporal discretization

AN US

CR IP T

300

The domain is discretized using regular blocks with an edge length of 0.2

311

m. Finer uniform grids (minimum edge length 0.05 m) as well as refined grids

312

close to fractures as proposed by Weatherill [46] were tested, too (results not

313

shown). The resulting convective pattern using all 3 different grids is visually

314

identical and the steady-state Sh deviate by less than 0.1%. Thus, the grid

315

consisting of 0.2 m blocks is sufficient.

ED

M

310

An adaptive time stepping scheme is used, where the time step size is

317

calculated based on the maximum concentration change during a time step

318

[12, 13, 39]. A maximum concentration change of 0.1 at any node is allowed.

319

Reducing the tolerance criteria to a maximum concentration change of 0.01

320

did not affect the results. Simulations were run until a steady state was

CE

PT

316

reached.

AC 321

18

ACCEPTED MANUSCRIPT

322

3. Results and discussion

323

3.1. Interfracture convection (mode 1) in a 2-dimensional regular fracture circuit

CR IP T

324

Interfracture convection is examined in a regular fracture circuit as shown

326

in Fig. 1a. We define a “regular fracture circuit” to consist of 4 fractures

327

that form an isolated continuous rectangle without any dead-end fractures.

328

A 2-dimensional representation is sufficient here, and the domain thickness

329

in y-direction can be reduced to one element (0.2 m). Thereby, intrafracture

330

convection (mode 2B) is deliberately precluded. Critical fracture aperture

331

for the onset of interfracture free convection (mode 1) in a 2-dimensional,

332

square regular fracture circuit as well as the strength of convection (Sh)

333

for various fracture apertures are evaluated numerically. The results are

334

later compared to those obtained from a 3-dimensional representation in

335

section 3.3. Table 2 shows the required fracture aperture (2b) for the onset

336

of interfracture convection (mode 1). Convection was detected at the critical

337

fracture aperture (2b)c = 1.4 · 10−5 m.

AC

CE

PT

ED

M

AN US

325

19

ACCEPTED MANUSCRIPT

Table 2: Onset and patterns of inter- and intrafracture convection

3d fracture circuit (Fig. 2a) −         

CR IP T

3d isolated fracture (Fig. 1b) − − − − −  (1)  (2)  (2)  (2)  (4)

AN US

(2b) 2d fracture circuit [10−5 m] (Fig. 1a) 1.0...1.3 − 1.4  1.5  1.6  1.7  1.8  1.9  2.0  2.5  3.0...5.0 

ED

M

- : no convection : interfracture convection (mode 1) along the fracture circuit as shown in e.g. Figs. 3 and 6,  (1): intrafracture convection mode 2B and number of convection cells per vertical fracture in brackets, as e.g. Fig. 5.

No convection occurs for (2b) < 1.4 · 10−5 m (Fig. 3a). For (2b) ≥

338

1.4 · 10−5 m however, circular flow along the fracture circuit occurs, shown

340

in Figs. 3b and 3c. With increasing fracture aperture, fracture velocity and

341

strength of convection increase. Therefore, solute is transported farther along

342

the fracture circuit, thereby increasing the concentration difference between

343

fracture and matrix at a transient stage of the simulation (transient results

AC

CE

PT

339

344

not shown). Thus, diffusion along the fracture-matrix interface increases. As

345

a result, the steady-state concentration along the fracture circuit is almost

346

uniform for high fracture apertures of (2b) = 5.0 · 10−5 m as shown in Fig.

347

3c. 20

ACCEPTED MANUSCRIPT

b)

c)

CR IP T

a)

AN US

Figure 3: Interfracture convection (mode 1): Solute distribution in the matrix is indicated by transparent colors (red high concentration, blue zero concentration) and grey isochlors. a) - no convection for (2b) < 1.4 · 10−5 m b)  Clockwise interfracture convection for (2b) = 1.4 · 10−5 m c)  Clockwise interfracture convection for (2b) = 5.0 · 10−5 m 348

Critical Rayleigh number for convection in a regular fracture circuit (mode

349

1)

The critical Rayleigh number for the onset of interfracture convection

351

(mode 1) in a 2-dimensional fracture circuit of variable size and aspect ra-

352

tio was determined numerically. A variety of regular fracture circuits with

353

variable circuit height h (distance between top and bottom fracture of the

354

circuit) and circuit width w (lateral distance between vertical fractures of

355

the circuit) were tested such that circuit area (h × w) and circuit aspect ratio

356

(w/h) varied. Fig. 4 shows the strength of convection (Sh) depending on the

357

number Rafr × (2b)/(hDr), where Rafr is calculated with Eq. (3) and where

PT

ED

M

350

fracture aperture is varied from (2b) = 5 · 10−6 m to 6 · 10−5 m. Varying

359

the fracture aperture of a 2-dimensional regular fracture circuit, all scenarios

360

below the critical number of 8.16 are stable, while convection was found in

361

all scenarios with a critical number larger than 8.16.

AC

CE

358

21

1.3

1.05

1.25

1.04

1.2 1.15 1.1 1.05

1.03 1.02 1.01 1

1 10

20

30

40

50

Rafr x (2b)/(hθτ)

0.99

60

5

CR IP T

Sherwood number

Sherwood number

ACCEPTED MANUSCRIPT

8.16

10

Rafr x (2b)/(hθτ)

15

362

AN US

Figure 4: Onset of interfracture convection mode 1: Strength of convection depending on Rayleigh number constant (Rafr × (2b)/(hrD )). Data points show circuits of different aspect ratio (width/height) and area (width×height). Black points represent circuits with large aspect ratio, white points represent circuits with small aspect ratio. Area is indicated by datapoint size.

Interestingly, the numerically determined critical number Rafr ×(2b)/(hDr) for the onset of interfracture convection (mode 1) is equal to half of the critical

364

number of 16.32 for the onset of intrafracture convection (mode 2B) as of Eq.

365

(4) determined by Wang et al. [44]. The critical Rayleigh number Rac-circ for

366

the onset of interfracture free convection (mode 1) in a 2-dimensional regular

367

fracture circuit can therefore be formulated as:

PT

ED

M

363

370

h (2b)

(7)

The corresponding formulation of the critical fracture aperture found here

is therefore:

AC

369

CE

368

Rac-circ = 8.16rD

(2b)c-circ =

s 3

rD D0 µ · 12 · 8.16 = 1.38 · 10−5 m β ∆c hρg H

(8)

That number is in excellent agreement with the numerically determined 22

ACCEPTED MANUSCRIPT

value of 1.4 · 10−5 m given in Tab. 2. Eqs. (7) and (8) are valid for a

372

regular fracture circuit embedded in a low permeable matrix, which inhibits

373

fluid flow but allows for diffusive solute transport, where upper and lower

374

boundaries of the matrix block are of fixed concentration.

CR IP T

371

Compared to previous studies which examined the onset of con-

376

vection in fracture circuits, our results are similar to those obtained

377

by Ramazanov [32] who gave a ratio between the critical Rayleigh

378

numbers for intra- and interfracture convection (mode 1) in a

379

square fracture circuit of approximately 2. While Ramazanov [32]

380

assumed that the thermal conductivity of fracture-filling material

381

and surrounding rock matrix are of the same order of magnitude,

382

his formulation of the critical Rayleigh number does not include

383

the thermal equivalent of the diffusivity ratio rD . For the special

384

case of a square fracture circuit with thin fractures, Ramazanov’s

385

[32] formulation can be approximated by Rac-circ = 8h/(2b).

M

Simmons et al.

[37] applied the classical formulation of the

ED

386

AN US

375

Rayleigh number for convection in a homogeneous porous medium

388

using the average permeability of the fractured medium and com-

389

pared it to the critical Rayleigh number of 4π 2 to determine the

390

onset of interfracture convection (mode 1).

391

that the Rayleigh number formulation corresponds to Eq. (3), the

392

critical Rayleigh number proposed by Simmons et al. [37] would

Reformulated such

AC

CE

PT

387

393

be almost 5 times higher than what was determined here, and

394

would thus underestimate the occurrence of interfracture convec-

395

tion. However, using the average permeability is appropriate to

396

determine the onset of free convection under the assumption of 23

ACCEPTED MANUSCRIPT

large fracture density and uniformly distributed fractures. For less

398

dense fracture networks, it was numerically shown by Vujevi´ c et

399

al. [43] and Yang et al. [48] that convection may occur in a frac-

400

tured medium containing interconnected fractures such as contin-

401

uous fracture circuits even if the Rayleigh number (Ra) averaged

402

over the entire medium is below the critical value [43, 48].

403

3.2. Intrafracture convection (mode 2B) in isolated vertical fractures

CR IP T

397

Intrafracture convection (mode 2B) is studied in two isolated parallel

405

vertical fractures embedded in a 3-dimensional matrix (Fig. 1b). The vertical

406

fractures are 6 m high, 6 m long in the y-direction and 6 m apart from each

407

other such that they are hydraulically disconnected by the low permeable

408

matrix. Fracture aperture is identical in both fractures and was varied in

409

order to determine the critical fracture aperture.

M

AN US

404

Onset of convection was detected at (2b) = 1.8·10−5 m (Table 2), which is

411

in very good agreement with the critical fracture aperture of (2b) = 1.74·10−5

412

m predicted by Eq. (5) based on the results by Wang et al. [44]. The

413

critical Rayleigh number Rac-fr as of Eq. (4) and the corresponding critical

414

fracture aperture for intrafracture convection (mode 2B) in Eq. (5), were

415

originally formulated for impermeable upper and lower fracture boundaries

416

of constant temperature (equivalent to constant concentration for haline con-

417

vection). This is to say, the vertical fracture is in contact with both, upper

418

and lower boundary. A new result obtained here is that Eqs. (4) and (5)

419

are also applicable if the upper and lower end of the fracture is separated

420

from the constant concentration boundary condition by an impermeable rock

421

matrix block that allows for solute diffusion. As a further result, it can be

422

expected that Eqs. (4) and (5) are also applicable to thermal convection

AC

CE

PT

ED

410

24

ACCEPTED MANUSCRIPT

423

in vertical fractures separated from the boundary conditions, because the

424

equations were originally formulated for a heat convection problem. The critical fracture aperture determined here for a single verti-

426

cal fracture can be compared to results obtained theoretically with

427

different considerations of the fracture-matrix solute exchange us-

428

ing (critical) Rayleigh numbers determined in previous studies.

429

Applying the classical Rayleigh criterion (Rac = 4π 2 ) using only

430

the fracture properties neglects the stabilizing effect of the ma-

431

trix and would considerably underestimate the required fracture

432

aperture for the onset of convection as demonstrated by Graf and

433

Therrien [21]. Thus, applying the closed box model by Beck [3]

434

leads to a decrease of the critical fracture aperture by almost two

435

orders of magnitude.

436

no-transport boundary conditions at fracture walls in combination

437

with the lower diffusion coefficient of the matrix (Dm ) in the formu-

438

lation of the Rayleigh number. Thus, the critical fracture aperture

439

further decreases, such that intrafracture convection seems even

440

more likely. When, on the other hand, solute loss along the vertical

441

fracture walls is considered by applying a Dirichlet fixed concen-

442

tration at the vertical fracture walls which decreases with depth

443

as e.g. by Kassoy and Cotte [23], the critical fracture aperture for

444

intrafracture convection increases by the factor of 4.

AN US

CR IP T

425

[37] used these no-flow,

AC

CE

PT

ED

M

Simmons et al.

25

ACCEPTED MANUSCRIPT

b)

c)

d)

CR IP T

a)

AN US

Figure 5: Intrafracture convection mode 2B: solute distribution and flow path within the fracture. a) − no convection (2b) < 1.8 · 10−5 m b)  (1) single convection cell per vertical fracture (2b) = 1.8 · 10−5 m c)  (2) two convection cells per vertical fracture (2b) = 2.0 · 10−5 m d)  (4) four convection cells per vertical fracture (2b) = 5.0 · 10−5 m

The simulated convective patterns of intrafracture convection (mode 2B)

446

are displayed in Fig. 5 and listed in Table 2. The convective pattern just

447

above the critical Rayleigh number is a single convection cell in

448

each vertical fracture (Fig. 5b). This is in agreement with findings

449

by Wang et al. [44], who predicted that a single cell is the most

450

likely mode of convection in a square fracture when the length ratio

451

between fracture aperture and adjacent matrix block width as well

452

as the conductivity ratio between fracture and matrix are small,

453

both of which is the case in the scenario displayed in Fig. 5b. Sim-

454

ilarly, Tournier et al. [41] numerically showed that a single cell is

455

the most likely steady-state convective pattern in a 2-dimensional

456

vertical fracture embedded in a heat-conducting matrix. This is

457

because conduction into and within the matrix dampens higher or-

458

der convective modes (many vertically elongated side-by-side cells)

459

such that higher Rayleigh numbers are required to increase the

460

number of convection cells. This is also demonstrated in Figs. 5c

AC

CE

PT

ED

M

445

26

ACCEPTED MANUSCRIPT

461

and 5d where the convective patterns for larger fracture apertures

462

and thus increasing Rayleigh numbers are shown, producing up to

463

four vertically elongated side-by-side convection cells (Fig. 5d). For brevity, the effect of fracture inclination is neglected in this

465

study. Theoretical investigations by Zhao et al. [51, 54] of con-

466

vection in an inclined fracture with constant temperature gradient

467

along the inclined fracture walls suggest that the critical Rayleigh

468

number of a vertical fracture is remarkably different from that of

469

an inclined fracture [54]. On the other hand, Tournier et al. [41],

470

who considered a conducting matrix adjacent to the fracture, as-

471

sumed that the effect of fracture inclination is smooth for moderate

472

inclinations. However, the effect of the fracture inclination on the

473

intrafracture convection mode should be considered in future nu-

474

merical models.

475

3.3. Convection in a 3-dimensional regular fracture circuit

ED

M

AN US

CR IP T

464

We combine geometries shown in Figs. 1a and 1b in a 3-dimensional

477

regular fracture circuit (Fig. 2a). This is the first numerical 3-dimensional

478

combined convection scenario that allows for interfracture convection mode

479

1, intrafracture convection mode 2B, and for a combination of both modes

480

(named “mode 12B”).

CE

481

PT

476

The onset of convection in this 3-dimensional regular fracture circuit was

determined numerically. Table 2 shows that, similar to interfracture convec-

483

tion (mode 1) in a 2-dimensional regular fracture circuit, the critical fracture

484

aperture for the onset of convection in a 3-dimensional regular fracture circuit

485

is 1.4 · 10−5 m. As the critical fracture aperture for intrafracture convection

AC 482

27

ACCEPTED MANUSCRIPT

486

(mode 2B) is much higher (1.8 · 10−5 m) than that for interfracture convec-

tion (1.4 · 10−5 m), it can be concluded that interfracture convection (mode

488

1) will dominate over intrafracture convection (mode 2B). It can further be

489

concluded that the formulation of the new critical Rayleigh number as of Eq.

490

(7) found for convection in a 2-dimensional fracture circuit is also valid for a

491

regular 3-dimensional fracture circuit. b)

c)

d)

AN US

a)

CR IP T

487

ED

M

Figure 6: Convection in combined fracture scenario: solute distribution and flow path in the fracture circuit. a) - no convection for (2b) < 1.4 · 10−5 m. b) to d)  interfracture convection only for (2b) = 1.4 · 10−5 m, 2.0 · 10−5 m, 5.0 · 10−5 m, respectively.

Fig. 6 shows convective patterns in 3-dimensional regular fracture cir-

493

cuits. Streamlines within each fracture are straight with no component in

494

the y-direction. This indicates that the regarded convection problems are

495

purely 2-dimensional and convective patterns correspond to those of purely

496

interfracture convection (mode 1).

CE

497

PT

492

Convective patterns remained to be 2-dimensional interfracture convec-

tion (mode 1) if the aspect ratio of the fracture circuit is modified by mul-

499

tiplying and dividing (one at a time) the extent of the fracture circuit in

500

each dimension by up to 1.5. Different scenarios and geometries of the re-

501

sulting fracture circuit dimensions with the highest and lowest aspect ratios

502

in each dimension are listed in Table 3. The range of aspect ratios given

AC

498

28

ACCEPTED MANUSCRIPT

503

in Table 3 covers the critical wave lengths for the lowest critical

506

number given in Eq. (4), which corresponds to a single intrafrac-

507

ture (mode 2B) convection cell in a vertical fracture with a fracture

508

length to fracture height ratio (L/h) of 0.7. However, the critical

509

Rayleigh numbers for larger aspect ratios (including aspect ratio

510

1 for a square fracture as used before) hardly deviate and do not

511

change the critical fracture apertures determined in section 3.2.

512

Ramazanov [32] determined the lowest critical Rayleigh number

513

for interfracture convection (mode 1) around a fracture circuit of

514

aspect ratio (2B/h) of 1.2. Thus, allowing both modes of convection

515

(mode 1 and mode 2B) in the potentially most favourable fracture

516

circuit geometries for convection ensures that no mode is favoured

517

by the geometry.

M

AN US

CR IP T

505

Rayleigh numbers in previous studies. Wang et al. [44] give a mini√ mum critical wave length of 2π for the minimum critical Rayleigh

504

Table 3 confirms that the strength of convection (Sh) increases with in-

519

creasing height, increasing area, decreasing width and decreasing aspect ratio

520

of the fracture circuit as previously shown by Vujevi´c et al. [43]. A new find-

521

ing is that the strength of convection is invariant to changes of the length of

522

a fracture circuit, such that Sh is identical for the 2-dimensional represen-

523

tation and for the 3-dimensional fracture circuit. Therefore, a 2-dimensional

524

representation captures all relevant processes for a regular fracture circuit,

AC

CE

PT

ED

518

525

because intrafracture convection appears to be absent in all cases.

29

ACCEPTED MANUSCRIPT

Table 3: Convection in a fracture circuit with varying geometry

strength of convection Sh 1.28 1.28 1.28 1.28 1.16 1.00 1.93 1.39 1.35

CR IP T

fracture circuit convective aspect ratio area pattern 2 (2B/h) m 1 36  1 36  1 36  1 36  0.67 24  0.67 24 0.75 48  1.33 48  1.2 48 

AN US

Fracture circuit fracture width×height×length aspect ratio in [m] (L/h) 6 × 6 × 0.2 (2D) 0.03 6 × 6 × 6 (3D) 1 6 × 6 × 4 (3D) 0.67 6 × 6 × 8 (3D) 1.5 4 × 6 × 6 (3D) 1 6 × 4 × 6 (3D) 1.5 6 × 8 × 6 (3D) 0.75 8 × 6 × 6 (3D) 1 7.2 × 6 × 4.2(3D) 0.7

ED

M

Strength of convection Sh depending on fracture circuit geometry for a 3-dimensional regular fracture circuit. Fracture aperture is 1.5 · 10−5 m. The 2-dimensional fracture circuit like in Fig. 1a is listed as reference.

Fig. 7 compares the strength of interfracture convection (mode 1) in a

527

2-dimensional regular fracture circuit, intrafracture convection (mode 2B) in

528

isolated vertical fractures, and convection in a 3-dimensional regular fracture

529

circuit using the Sherwood number Sh (Eq. (6)). Fig. 7 shows that the

530

Sh curve of the 3-dimensional regular fracture circuit is identical to that of

531

the 2-dimensional circuit. This indicates that all features of interfracture

532

convection in a regular fracture circuit can be captured with a 2-dimensional

AC

CE

PT

526

533

representation. Intrafracture convection (mode 2B) in two isolated vertical

534

fractures is weaker (smaller Sh) and therefore of subordinate importance.

30

1.3

1.05

1.25

1.04

1.2 1.15 1.1 1.05

1.03 1.02 1.01 1

1 10

20

30

40

Rafr x (2b)/(h rD)

50

60

0.99

8.16

5

10

Rafr x (2b)/(h rD)

CR IP T

Sherwood number

Sherwood number

ACCEPTED MANUSCRIPT

15

535

536

3.4. Effect of matrix permeability

AN US

Figure 7: Strength of convection (Sherwood number Sh) for inter- and intrafracture convection separately and for the combined convection scenario. Colored symbols indicate convective patterns in each case as in Tab. 2.

In all cases described above, the rock matrix permeability is km = 1 ·

10−18 m2 such that the matrix is essentially impermeable. We increase the

538

permeability of the matrix around the 3-dimensional regular fracture circuit

539

(Fig. 2a) and its 2-dimensional representation (Fig. 1a) to km = 1 · 10−17

M

537

m2 , km = 5 · 10−17 m2 and km = 1 · 10−16 m2 , respectively, to demonstrate

541

the effect of matrix permeability.

ED

540

Fig. 8 shows the strength of convection (Sh) depending on km and (2b).

543

Higher matrix permeability enhances convection because convective flow is

544

no longer restricted to fractures and to the fracture circuit, but can also

545

develop in the matrix. Fig. 8 shows that critical fracture aperture remains

546

at 1.4 · 10−5 m for km = 1 · 10−17 m2 and decreases to 0.9 · 10−5 m for

CE

PT

542

km = 5 · 10−17 m2 . For km = 1 · 10−16 m2 , the Rayleigh number for the

548

matrix is above the critical Rayleigh number, such that a critical fracture

549

aperture can not be determined as the occurrence of free convection does

550

not depend on the presence of fractures. In the last case of high matrix

551

permeability, the classical Rayleigh criterion may be applied to

552

determine the onset of convection in fractured rock using the bulk

AC

547

31

ACCEPTED MANUSCRIPT

permeability as proposed by Simmons et al. [37]. This is only

554

possible in a system where fracture permeability is small relative

555

to matrix permeability such that flow is predominately controlled

556

by matrix permeability as stated by Graf and Therrien [21].

CR IP T

553

2.5

2

AN US

Sherwood number (Sh)

3

-16

km=1·10 m² -17

km=5·10 m²

1.5

-17

km=1·10 m² -18

km=1·10 m²

onset of convection 1

1.5

2

2.5

3

3.5

4

fracture aperture (2b)

ED

0.5

M

1

4.5

5 5.5 -5 · 10 m

558

3.5. More complex scenarios In addition to studying a regular fracture circuit, we also give an example

CE

557

PT

Figure 8: Strength of convection (Sh) in the combined convection scenario depending on km and (2b).

of a more complex geometry. We thereto reduce the extent of the fractures

560

to 4 m and add two more fractures, such that the fracture network consists

AC

559

561

of 6 instead of 4 fractures (Fig. 2b). The 6 fractures are all connected and

562

form a fracture circuit with some “dead-end fractures” at the corners of the

563

fracture network. All fractures traverse the entire domain in y-direction. The

564

location of the center of each fracture is given in Table 4. 32

ACCEPTED MANUSCRIPT

Table 4: Fracture location in the complex fracture circuit consisting of 6 fractures

CR IP T

Fracture x-coordinate z-coordinate orientation of fracture center [m] vertical 2.2 6.6 vertical 4.6 4.6 vertical 5.8 7.2 horizontal 3.0 5.0 horizontal 6.0 6.0 horizontal* 4.0 8.0/9.0

AN US

* the vertical position of this fracture is variable and is 8.0 m in the 3 m high fracture circuit and 9.0 m in the 4 m high fracture circuit.

Based on the observations from interfracture convection (mode 1) and in-

566

trafracture convection (mode 2B), we hypothesize that the critical Rayleigh

567

number depends on either the fracture height or the height of the fracture

568

circuit. Studying a case of constant fracture height and varying circuit height

569

allows us to evaluate whether (i) the fracture circuit is the dominating struc-

570

ture in the fracture network such that predicting the onset of convection is

571

possible using a 2-dimensional representation, or (ii) the vertical fractures

572

in the network promote interfracture convection such that a 3-dimensional

573

representation is required. To answer this question, the location of the up-

574

permost horizontal fracture is variable (red arrow in Fig. 2b), such that the

575

fracture circuit is either 3 m high (Fig. 9a), or 4 m high (Fig. 9b).

AC

CE

PT

ED

M

565

576

If the circuit height is only 3 m and thus smaller than the height of each

577

vertical fracture, intrafracture convection (mode 2B) develops (Fig. 9a). The

578

resulting convection pattern is combined inter- and intrafracture convection

579

(mode 12B) where intrafracture convection cells in the vertical fractures are 33

ACCEPTED MANUSCRIPT

hydraulically connected along the horizontal fractures as shown in Fig. 9a.

581

However, if the fracture circuit height equals 4 m like the individual fracture

582

height, interfracture convection (mode 1) dominates (Fig. 9b). Intrafrac-

583

ture convection (mode 2B) is only observed in the dead-end fractures when

584

the fracture aperture exceeds the critical value for intrafracture convection

585

(2b)c = 2.0 · 10−5 m as of Eq. (5). b)

M

AN US

a)

CR IP T

580

PT

ED

Figure 9: More complex fracture network with 6 fractures: solute distribution and flow paths in the fracture circuit. Fracture length is 4 m, fracture aperture (2b) = 2.0 · 10−5 m. a),  (1) Fracture circuit height is 3 m, and lower than the height of individual vertical fractures, combined intra- and interfracture convection (mode 12B) develops. b) Fracture circuit height is 4 m, identical to the height of individual vertical fractures, interfracture convection mode 1 dominates.

For both heights (3 and 4 m) of the fracture circuit, a 2-dimensional cross-

587

section of the 3-dimensional network is simulated in order to compare the

588

results with respect to convective pattern and strength of convection (Sh).

589

The results are summarized in Table 5. If fracture and circuit height are not

AC

CE

586

590

identical (4 m vs. 3 m), predictions based on a 2-dimensional representation

591

are critical because the 2-dimensional model neglects intrafracture convec-

592

tion (mode 2B) in dead-end fractures and combined convection (mode 12B)

593

such that the occurrence and the strength of convection are underestimated. 34

ACCEPTED MANUSCRIPT

This is because in the 3 m high fracture circuit with combined convection

595

(mode 12B), the dead-end fractures vertically elongate the fracture circuit

596

and promote convection by increasing the vertical density contrast. If, how-

597

ever, fracture and circuit height are identical (both 4 m), the 2-dimensional

598

cross-section is an appropriate simplification showing the same convective

599

pattern and the same Sherwood number for (2b) ≤ 2 · 10−5 m because there

600

is no convection in the dead-end fractures.

CR IP T

594

AN US

Table 5: Observed convective patterns in a fracture circuit consisting of 6 fractures

3-dimensional network deviance in Sh fracture circuit height between 2d and 3d 3m 4m 3m 4m  0% ,  (1)  4% 0% ,  (1)  5% 0% ,  (1) * 6% 1%

ED

M

Fracture 2-dimensional network aperture fracture circuit height ·10−5 [m] 3 m 4m 1.5 a 1.6 1.7  b 1.8  1.9   c 2.0   a

onset of mode 1 convection in a 4 m high circuit predicted by Eq. (7). onset of mode 1 convection in a 3 m high circuit predicted by Eq. (7). c onset of mode 2B convection in a 4 m high fracture predicted by Eq. (4). : interfracture convection along the fracture circuit (mode 1). *: only interfracture convection in the circuit, but intrafracture convection in the dead-end fractures. ,  (1): combination of interfracture and intrafracture convection (mode 12B).

AC

CE

PT

b

601

Results in Table 5 show that predicting the onset of convection in com-

602

plex fracture networks using theoretical considerations (Rayleigh criterion)

603

are difficult. Beside the effect of long vertical fractures promoting convection,

604

additional fracture intersections along the fracture circuit disturb interfrac35

ACCEPTED MANUSCRIPT

ture flow and thus inhibit convection, such that complex fracture net-

606

works impede the prediction of free convective patterns and vigor

607

of convection as already stated by Shikaze et al. [34].

608

4. Summary and conclusions

CR IP T

605

Continuing prior work on convective flow modes in fractured rock by

610

Wang et al. [44], Malkovsky and Pek [26], Ramazanov [32], or Simmons et

611

al. [37], this study examines inter- and intrafracture haline free convection

612

and their combination. The onset of convection, convective patterns and the

613

dominating mode of convection are determined. The key new findings of this

614

study are:

AN US

609

1. The critical Rayleigh number for the onset of convection in a 2- and

616

3-dimensional regular fracture circuit is determined numerically to be

617

half of the critical Rayleigh number for the onset of convection in a

618

single vertical fracture determined by Wang et al. [44].

ED

M

615

2. The most likely mode of convection in a regular fracture circuit is inter-

620

fracture convection along the circuit (mode 1). When both inter- and

621

intrafracture convection are theoretically possible, the onset of con-

622

vection is marked by the onset of interfracture convection (mode 1). Convective patterns correspond to streamlines along the circuit.

3. Interfracture convection (mode 1) is stronger in terms of mass transport

AC

624

CE

623

PT

619

625

than intrafracture convection (mode 2B). Sherwood numbers for in-

626

trafracture convection in two vertical fractures are considerably smaller

627

than for interfracture convection in a regular fracture circuit of the same

628

height. 36

ACCEPTED MANUSCRIPT

4. In more complex fracture networks, convection is hard to be theoret-

630

ically predicted using Rayleigh criteria. Vertical “dead-end fractures”

631

exceeding the fracture circuit promote convection, whereas additional

632

fracture intersections along the fracture circuit constrain convection.

CR IP T

629

5. A 2-dimensional representation of a 3-dimensional fracture network

634

only yields good results in predicting the onset and strength of con-

635

vection if the network is composed of regular fracture circuits, where

636

the circuit height equals the individual fracture height. Otherwise, 2-

637

dimensional representations tend to underestimate the occurrence and

638

strength of convection, as intrafracture convection (mode 2B) and com-

639

bined convection (mode 12B) are neglected.

AN US

633

Considering the easier onset and the higher mass transport rates of in-

641

terfracture convection (mode 1) compared to intrafracture convection (mode

642

2B), interfracture convection along a fracture circuit (mode 1) is the most

643

likely mode of convection in a fracture network consisting of regular frac-

644

ture circuits. However, more complex fracture networks, which not exclu-

645

sively consist of regular fracture circuits, challenge this rule and require 3-

646

dimensional representations. A fully 3-dimensional modelling approach of

647

free convection in fractured rock remains challenging in terms of computa-

648

tional power. It took up to 10 times more CPU time to run a 3-dimensional

649

complex scenario relative to its 2-dimensional representation. A series of

AC

CE

PT

ED

M

640

650

2-dimensional cross-sections perpendicular to the fracture planes of a frac-

651

ture network (covers mode 1) combined with a theoretical evaluation of the

652

existence of intrafracture free convection in isolated vertical fractures using

653

Eq. (4) might be a practical solution to at least approximate the probability 37

ACCEPTED MANUSCRIPT

654

of inter- and intrafracture convection in fractured rock. Nevertheless, this

655

approach will not be able to examine combined convection (mode 12). This study focusses on simple deterministic fracture networks,

657

consisting of few fractures and therefore represents an academic

658

configuration rather than a real-world field setting. However, the

659

new findings presented may help to better assess the probability

660

of free convection in various fracture networks. Furthermore, the

661

conclusions drawn in this study lay the foundation for further re-

662

search on more complex fracture networks.

AN US

CR IP T

656

On the basin scale, interfracture convection enables convective

664

flow and transport over larger distances than intrafracture convec-

665

tion, especially if fracture length is small compared to the consid-

666

ered scale. Cui et al. [8] stated that isolated vertical faults do not

667

affect basin-scale free convection, and Yang et al. [49] showed that

668

open vertical faults may determine convective patterns in a sedi-

669

mentary basin only if they are hydraulically connected via highly

670

permeable horizontal sediment layers. Thus, further investigation

671

of free convection in complex fracture networks is required.

PT

ED

M

663

In field studies, the exact geometry of fracture networks is often unknown.

673

Parameters like average and maximum vertical extent and permeability of

674

fractures should be estimated to be able to calculate a (critical) fracture

675

Rayleigh number (Eqs. (3) and (4)), such that the existence of intrafracture

AC

CE

672

676

convection (mode 2) can be assessed. Assessment of interfracture convection

677

(mode 1) is more difficult, because the existence and geometry of fracture

678

circuits is harder to detect than the existence of fractures. An elementary

679

relationship between the strength of interfracture convection and average 38

ACCEPTED MANUSCRIPT

fracture length and fracture density is given in Vujevi´c et al. [43]. How-

681

ever, more statistically oriented research is required to better quantify these

682

relationships.

AC

CE

PT

ED

M

AN US

CR IP T

680

39

ACCEPTED MANUSCRIPT

683

Acknowledgement This research was supported by the Deutsche Forschungsgemeinschaft

685

(DFG) under grant number GR 3463/2-1. We thank the editorial board of

686

AWR (D. Andrew Barry) for handling our manuscript. We also thank four

687

anonymous reviewers whose constructive comments have helped to improve

688

the manuscript.

AC

CE

PT

ED

M

AN US

CR IP T

684

40

ACCEPTED MANUSCRIPT

APPENDIX Nomenclature

AC

CR IP T

fracture aperture critical fracture aperture in a single vertical fracture critical fracture aperture in the fracture circuit fracture spacing relative solute concentration hydrodynamic dispersion tensor in the matrix hydrodynamic dispersion tensor in the fracture effective diffusion coefficient in the matrix free solution diffusion coefficient unit vector in vertical direction gravitational acceleration height of fracture domain height equivalent fresh water head identity matrix permeability of (isotropic) porous medium fracture permeability matrix permeability freshwater hydraulic conductivity tensor of porous medium freshwater hydraulic conductivity tensor of porous medium length of the domain (y-direction) fluid pressure Darcy flux vector in the porous matrix Darcy flux vector in the fracture diffusivity ratio between matrix and fracture Rayleigh number (homogeneous porous medium) Rayleigh number in the fracture critical Rayleigh number critical Rayleigh number in a regular fracture circuit critical Rayleigh number in a single vertical fracture

PT

ED

[L] [L] [L] [L] [-] [L2 T−1 ] [L2 T−1 ] [L2 T−1 ] [L2 T−1 ] [-] [LT−2 ] [L] [L] [L] [-] [L2 ] [L2 ] [L2 ] [LT−1 ] [LT−1 ] [L] [ML−1 T−2 ] [LT−1 ] [LT−1 ] [-] [-] [-] [-] [-] [-]

CE

691

(2b) (2b)c (2b)c,circ (2B) c D Dfr Dm D0 ez g h H h0 I k kfr km K0 K0,fr L p q qfr rD Ra Rafr Rac Rac,circ Rac,fr

AN US

690

M

689

41

ACCEPTED MANUSCRIPT

CR IP T

AN US

M

Governing equations

AC

694

CE

693

specific storage of porous medium specific storage of the fracture Sherwood number time total solute mass flux through top of the domain width of the domain in x-direction horizontal coordinates vertical coordinate longitudinal dispersion coefficient longitudinal dispersion coefficient in the fracture transverse dispersion coefficient transverse dispersion coefficient in the fracture solutal expansion coefficient fluid compressibility matrix compressibility matrix porosity dynamic viscosity of the fluid fluid density reference fluid density tortuosity incline of fracture (ϕ for a vertical fracture is 0) exchange term for fluid in the porous matrix exchange term for fluid in the fracture exchange term for solute in the porous matrix exchange term for solute in the fracture

ED

[L−1 ] [L−1 ] [-] [T] [ML−2 T−1 ] [L] [L] [L] [L] [L] [L] [L] [-] [LT2 M−1 ] [LT2 M−1 ] [-] [ML−1 T−1 ] [ML−3 ] [ML−3 ] [-] [-] [T−1 ] [LT−1 ] [T−1 ] [LT−1 ]

PT

692

SS SS,fr Sh t j W x, y z αL αL,fr αT αT,fr β γ fl γm θ µ ρ ρ0 τ ϕ Ωfluid Ωfluid,fr Ωsolute Ωsolute,fr

695

The Darcy equation describes the Darcy flux in a porous medium as a

696

function of both, the equivalent freshwater head h0 and the relative fluid

42

ACCEPTED MANUSCRIPT

density (ρ − ρ0 )/ρ0 [-] [2]: q = −K0



ρ(c) − ρ0 ez ∇h0 + ρ0



(9)

CR IP T

697

698

Assuming Darcy flow in the fracture, the Darcy equation for density-driven

699

flow in the fracture reads [21]:

700 701

AN US

  ρ(c) − ρ0 qfr = −K0,fr ∇h0 + ez cosϕ ρ0

The conservation equations for fluid mass in the 3-dimensional porous medium (11) and in the 2-dimensional fracture (12) are [10, 21]: 

 ρ(c) − ρ0 ∂h0 ∇ · K0 (∇h0 + ez ) − SS = Ωfluid ρ0 ∂t

(11)

M

   ρ(c) − ρ0 ∂h0 (2b) ∇ · K0,fr (∇h0 + ez cosγ) − SS,fr = (2b)Ωfluid,fr (12) ρ0 ∂t The specific storage in matrix (SS ) and fracture (SS,fr ) are defined as [14, 21]:

ED

703



PT

702

(10)

(13)

SS,fr = ρ0 gγ fl

(14)

CE

704

SS = ρ0 g(γ m + θγ fl )

The conservation equations for solute mass in the porous medium (15) and

706

in the fracture (16) are [10, 21]:

AC

705

707

∇ · (θD∇c − qc) − 

∂(θc) = Ωsolute ∂t

 ∂c (2b) ∇ · (Dfr ∇c − qfr c) − = (2b)Ωsolute,fr ∂t 43

(15)

(16)

ACCEPTED MANUSCRIPT

708

The hydrodynamic dispersion tensor in the matrix (D) and in the fracture

709

(Dfr ) are defined as [2, 20]:

710

and Dfr = (αL,fr − αT,fr )

(17)

qfr qT fr + αT,fr |qfr |I + D0 I |qfr |

(18)

where the effective diffusion coefficient in the porous matrix is [39]:

AN US

711

qqT + αT |q|I + Dm I |q|

CR IP T

θD = (αL − αT )

Dm = θτ D0

713 714

Constitutive equations

Porous medium freshwater hydraulic conductivity tensor in the matrix K0 and in the fracture K0,fr are given by [2, 5, 20]:

M

712

ED

K0 =

715

km Iρ0 g µ

(20)

ρ0 g µ

(21)

PT

K0,fr = kfr I

where fracture permeability kfr for an open fracture is defined as [2, 5] kfr =

(2b)2 12

(22)

AC

CE

716

(19)

717 718

Fluid density ρ is determined by relative concentration c and is calculated

using a linear density function [11, 37]: ρ(c) = ρ0 (1 + βc)

44

(23)

ACCEPTED MANUSCRIPT

The dependency of viscosity on salinity is neglected as it has been shown that

720

the effect is negligible for density-dependent flow and transport compared to

721

the density effect of salinity [18]. Thus, fluid viscosity is assumed to be

722

constant in accordance with most studied applying Rayleigh theory [26, 32,

723

34, 37, 44, 50].

AC

CE

PT

ED

M

AN US

CR IP T

719

45

ACCEPTED MANUSCRIPT

724

References [1] P. Alt-Epping, C. Zhao, Reactive mass transport modelling of a three-

726

dimensional vertical fault zone with a finger-like convective flow regime,

727

J. Geochemical Explor.(2010),106(1-3):8-23.

728

http://dx.doi.org/10.1016/j.gexplo.2009.12.007.

730

731

[2] J. Bear, Dynamics of fluids in porous media, American Elsevier Publishing Co., New York, 1972.

AN US

729

CR IP T

725

[3] J.L. Beck, Convection in a box of porous material saturated with fluid,

732

Phys. Fluids(1972),15:1377-1383.

733

http://dx.doi.org/10.1063/1.1694096

[4] A. Bejan, K.R. Khairy, Heat and mass transfer by natural convection

735

in a porous medium, Int. J. Heat Mass Tran.(1985),28(5):909-918.

736

http://dx.doi.org/10.1016/0017-9310(85)90272-8.

ED

737

M

734

[5] B. Berkowitz, Characterizing flow and transport in fractured geological media: A review, Adv. Water Resour.(2002),25:861-884.

739

http://dx.doi.org/10.1016/S0309-1708(02)00042-8.

741

rocks: Correlation to porosity and hydraulic conductivity, J. Contam. Hydrol.(2001),53:85100.

AC

742

[6] T.B. Boving, P. Grathwohl, Tracer diffusion coefficients in sedimentary

CE

740

PT

738

743

hhtp://dx.doi.org/10.1016/S0169-7722(01)00138-3.

744

[7] J.P. Caltagirone, Convection in a porous medium, in: J. Zierep, H.

745

Oertel jr. (Eds.), Convective transport and instability phenomena, G.

746

Braun, Karlsruhe(1982), pp.199-232. 46

ACCEPTED MANUSCRIPT

[8] T. Cui, J. Yang, I.M. Samson, Numerical modeling of hydrothermal

748

fluid flow in the Paleoproterozoic Thelon Basin, Nunavut, Canada, J.

749

Geochem. Explor.(2010),106,69-76.

750

http://dx.doi.org/10.1016/j.gexplo.2009.12.008.

751 752

CR IP T

747

[9] S.H. Davis, Convection in a box: linear theory, J. Fluid Mech.(1967),30:465478. http://dx.doi.org/10.1017/S0022112067001545

[10] H.-J.G. Diersch, O. Kolditz, Coupled groundwater flow and trans-

754

port: 2. Thermohaline and 3D convection systems, Adv. Water Re-

755

sour.(1998),21:401425. http://dx.doi.org/10.1016/S0309-1708(97)00003-

756

1.

757

AN US

753

[11] H.-J.G. Diersch, O. Kolditz, Variable-density flow and transport in porous media: approaches and challenges, Adv. Water Resour.(2002),25:899-

759

944.

760

http://dx.doi.org/10.1016/S0309-1708(02)00063-5.

ED

M

758

[12] P.A. Forsyth, P.H. Sammon, Practical considerations for adaptive im-

762

plicit methods in reservoir simulation, J. Comput. Phys.(1986),62:265-

763

281.

764

http://dx.doi.org/10.1016/0021-9991(86)90127-0.

CE

[13] P.A. Forsyth, M.C. Kropinski, Monotonicity considerations for saturated-

AC

765

PT

761

766

unsaturated subsurface flow, SIAM J. Sci. Comp.(1997),18:1328-1354.

767

http://dx.doi.org/10.1137/S1064827594265824.

768 769

[14] E.O. Frind, Simulation of long-term transient density-dependent transport in groundwater, Adv. Water Resour.(1982),5:7388. http://dx.doi.org/10.1016/030947

ACCEPTED MANUSCRIPT

770

1708(82)90049-5. [15] G.Z. Gershuni, E.M. Zhukhovitskii, Convective stability of incompress-

772

ible fluids, translated from Russian by D. Louvish, Keter Publishing

773

House Jerusalem Ltd., Jerusalem, 1976.

CR IP T

771

[16] G. Garven, S.W. Bull, R.R. Large, Hydrothermal fluid flow models

775

of stratiform ore genesis in the McArthur Basin, Northern Territory,

776

Australia, Geofluids(2001),1:289-311.

777

http://dx.doi.org/10.1046/j.1468-8123.2001.00021.x.

AN US

774

[17] T. Graf, Modeling coupled thermohaline flow and reactive transport in

779

discretely-fractured porous media, PhD thesis, Universit´e Laval(2005),

780

209pp. Available at .

782

[18] T. Graf, M.C. Boufadel, Effect of viscosity, capillarity and grid spacing on thermal variable-density flow, J. Hydrol.(2011),400:4157.

ED

781

M

778

[19] T. Graf, R. Therrien, Variable-density groundwater flow and solute

784

transport in porous media containing nonuniform discrete fractures,

785

Adv. Water Resour.(2005),28:1351-1367.

786

http://dx.doi.org/10.1016/j.advwatres.2005.04.011. [20] T. Graf, R. Therrien, Variable-density groundwater flow and solute transport in irregular 2D fracture networks, Adv. Water Resour.(2007),

AC

788

CE

787

PT

783

789

790 791

30:455-468. http://dx.doi.org/10.1016/j.advwatres.2006.05.003.

[21] T. Graf, R. Therrien, Stable-unstable flow of geothermal fluids in fractured rock, Geofluids(2009),9:138-152. 48

ACCEPTED MANUSCRIPT

792

http://dx.doi.org/10.1111/j.1468-8123.2008.00233.x. [22] C.W. Horton, F.T. Rogers, Convection currents in a porous medium,

794

J. Appl. Phys.(1945),16:367370. http://dx.doi.org/10.1063/1.1707601.

CR IP T

793

795

796

[23] D.R. Kassoy, B. Cotte, The effect of sidewall heat loss on convection in a saturated porous vertical slab, J. Fluid Mech.(1985),152:361-378.

798

http://dx.doi.org/10.1017/S0022112085000738.

799

AN US

797

[24] E.R. Lapwood, Convection of a fluid in a porous media, Proc. Cam-

800

bridge Philos. Soc.(1948),44(4):508-521.

801

http://dx.doi.org/10.1017/S030500410002452X.

[25] R.P. Lowell, C.-T. Shyu, On the onset of convection in a water-saturated

803

porous box: effect of conducting walls, Lett. heat mass trans.(1978),5:

804

371-378. http://dx.doi.org/10.1016/0094-4548(78)90046-2.

ED

805

M

802

[26] V.I. Malkovsky, A.A. Pek, Conditions for the onset of thermal convection of a homogeneous fluid in a vertical fault, Petrology+(1997),5(4):381-

807

387.

809

fluid in open vertical faults, Izvestiya, Physics of the Solid Earth(2004),40(8):672679.

AC

810

[27] V.I. Malkovsky, A.A. Pek, Onset of thermal convection of single-phase

CE

808

PT

806

811

[28] H.D. Murphy, Convective instabilities in vertical fractures and faults,

812

J. Geophys. Res.(1979),84:6121-6130.

813

http://dx.doi.org/10.1029/JB084iB11p06121. 49

ACCEPTED MANUSCRIPT

[29] C.E. Neuzil, How permeable are clays and shales?, Water Resour.

815

Res.(1994),30:145150.

816

http://dx.doi.org/10.1029/93WR02930.

817

[30] D.A. Nield, Onset of thermohaline convection in a porous medium,

818

Water Resour. Res.(1968),4(3):553-560.

819

http://dx.doi.org/10.1029/WR004i003p00553.

821

822

[31] D.A. Nield, A. Bejan, Convection in porous media, Springer, New York,

AN US

820

CR IP T

814

2013.

[32] M. Ramazanov, Convective stability of a fluid in two hydrodynamically connected vertical faults, Izvestiya, Physics of the Solid Earth(2003),39(12):1014-

824

1020.

M

823

[33] J.M. Sharp, Jr., T.R. Fenstemaker, C.T. Simmons, T.E. McKenna, J.K.

826

Dickinson, Potential thermohaline convection in a shale-rich sedimen-

827

tary basin: Example from the Gulf of Mexico Basin in South Texas.

828

AAPG Bull.(2001),85:2089-2110.

PT

ED

825

[34] S.G. Shikaze, E.A. Sudicky, F.W. Schwartz, Density-dependent solute

830

transport in discretely-fractured geological media: is prediction possi-

831

ble? J. Contam. Hydrol.(1998),34:273-291. http://dx.doi.org/10.1016/S0169-7722(98)00080-1.

AC

832

CE

829

833

[35] C.T. Simmons, K.A. Narayan, R.A. Wodding, On a test case for density-

834

dependent groundwater flow and solute transport models: The salt lake

835

problem, Water Resour. Res.(1999),35:3607-3620.

836

http://dx.doi.org/10.1029/1999WR900254. 50

ACCEPTED MANUSCRIPT

837

[36] C.T. Simmons, Variable density groundwater flow: From current challenges to future possibilities, Hydrogeol. J.(2005),13:116-119.

839

http://dx.doi.org/10.1007/s10040-004-0408-3

840

CR IP T

838

[37] C.T. Simmons, J.M. Sharp, D.A. Nield, Modes of free convection in

841

fractured low-permeability media, Water Resour. Res.(2008),44:W03431.

842

http://dx.doi.org/10.1029/2007WR006551.

[38] M.A. Simms, G. Garven, Thermal convection in faulted extensional

844

sedimentary basins: theoretical results from finite-element modeling,

845

Geofluids(2004),4,109-130.

846

http://dx.doi.org/10.1111/j.1468-8115.2004.00069.x.

847

AN US

843

[39] R. Therrien, R.G. McLaren, E.A. Sudicky, S.M. Panday, HydroGeoSphere - A Three-dimensional numerical model describing fully-integrated

849

subsurface and surface flow and solute transport, University of Wa-

850

terloo, Canada: Groundwater Simulations Group, 2010. Available at

851

http://www.ggl.ulaval.ca/fileadmin/ggl/documents/rtherrien/hydrogeosphere.pdf.

ED

M

848

PT

852

[40] R. Thiesen, K. Scheel, H. Diesselhorst, Untersuchungen u ¨ber die ther-

854

mische Ausdehnung von festen und tropfbar fl¨ ussigen K¨orpern - Bestim-

855

mung der Ausdehnung des Wassers f¨ ur die zwischen 0 und 40 liegenden Temperaturen, Wiss. Abhandlungen Der Phys. Reichsanst.(1900),3:170.

AC

856

CE

853

857

858

[41] C. Tournier, P. Genthon, M. Rabinowicz, The onset of natural con-

859

vection in vertical fault planes:consequences for the thermal regime in

51

ACCEPTED MANUSCRIPT

crystalline basements and for heat recovery experiments, Geophys. J.

861

Int.(2000),140:500-508.

862

http://dx.doi.org/10.1046/j.1365-246X.2000.00041.x.

CR IP T

860

[42] C.I. Voss, C.T. Simmons, N.I. Robinson, Three-dimensional bench-

864

mark for variable-density flow and transport simulation: matching

865

semi-analytical stability modes for steady unstable convection in an

866

inclined porous box, Hydrogeol. J.(2009),18:5-23.

867

http://dx.doi.org/10.1007/s10040-009-0556-6.

AN US

863

868

[43] K. Vujevi´c, T. Graf, C.T. Simmons, A.D. Werner, Impact of frac-

869

ture network geometry on free convective flow patterns, Adv. Water

870

Res.(2014),71:65-80. http://dx.doi.org/10.1016/j.advwatres.2014.06.001.

M

871

[44] M. Wang, D.R. Kassoy, P.D. Weidman, Onset of convection in a vertical

873

slab of saturated porous media between two conducting impermeable

874

blocks, Int. J. Heat Mass Tran.(1987),30(7):1331-1341.

875

http://dx.doi.org/10.1016/0017-9310(87)90165-7.

877 878

PT

dependent groundwater models: two-dimensional steady state unstable convection in infinite, finite and inclined porous layers, Adv. Water Resour.(2004),27:547-562. http://dx.doi.org/10.1016/j.advwatres.2004.01.003.

AC

879

[45] D. Weatherill, C.T. Simmons, C.I. Voss, N.T. Robinson, Testing density-

CE

876

ED

872

880

881

[46] D. Weatherill, T. Graf, C.T. Simmons, P.G. Cook, R. Therrien, D.A.

882

Reynolds, Discretizing the fracture-matrix interface to simulate solute

52

ACCEPTED MANUSCRIPT

883

transport. Ground Water(2008),46(4):606-615.

884

http://dx.doi.org/10.1111/j.1745-6584.2007.00430.x. [47] Y. Xie, C.T. Simmons, A.D. Werner, Speed of free convective fingering

CR IP T

885 886

in porous media, Water Resour. Res.(2011),47:W11501,16. http://dx.doi.org/

887

10.1029/2011WR010555.

[48] J. Yang, K. Latychev, R.N. Edwards, Numerical computation of hy-

889

drothermal fluid circulation in fractured earth structures, Geophys. J.

890

Int.(1998),135:627-649.

891

http://dx.doi.org/10.1046/j.1365-246X.1998.00669.x

AN US

888

[49] J. Yang, R.R. Large, S.W. Bull, Factors controlling free thermal con-

893

vection in faults in sedimentary basins: implications for the formation

894

of zinc-lead mineral deposits, Geofluids(2001),4:237-247.

895

http://dx.doi.org/10.1111/j.1468-8123.2004.00084.x

ED

M

892

[50] C. Zhao, B.E. Hobbs, H. M¨ uhlhaus, A. Ord, G. Lin, Convective insta-

897

bility of 3-D fluid-saturated geological fault zones heated from below,

898

Geophys. J. Int.(2003),155:213-220. http://dx.doi.org/10.1046/j.1365-

899

246X.2003.02032.x

[51] C. Zhao, B.E. Hobbs, A. Ord, S. Peng, H.B. M¨ uhlhaus, L. Liu,

AC

901

CE

900

PT

896

902

Theoretical investigation of convective instability in inclined

903

and fluid-saturated three-dimensional fault zones, Tectono-

904

physics (2004),387:47-64. http://dx.doi.org/10.1016/j.tecto.2004.06.007.

905

53

ACCEPTED MANUSCRIPT

906

[52] C. Zhao, B.E. Hobbs, A. Ord, S. Peng, H. M¨ uhlhaus, L.Liu, Doublediffusive-driven convective instability of three-dimensional fluid-saturated

908

geological fault zones heated from below, Math. Geol.(2005),37(4):373-

909

391. http://dx.doi.org/10.1007/s11004-005-5954-2

CR IP T

907

[53] C. Zhao, B.E. Hobbs, A. Ord, M.K¨ uhn, H. M¨ uhlhaus, S. Peng, Nu-

911

merical simulation of double-diffusive driven convective flow and rock

912

alteration in three-dimensional fluid-saturated geological fault zones,

913

Comp. Methods Appl. Mech. Engrg.(2006),195:2816-2840.

914

http://dx.doi.org/10.1016/j.cma.2005.07.008

916

917

[54] C. Zhao, B.E. Hobbs, A. Ord, Convective and advective heat transfer in geological systems, Springer, Berlin, 2008. [55] C. Zhao, B.E. Hobbs, A. Ord, Fundamentals of computational

M

915

AN US

910

geoscience: Numerical methods and algorithms, Springer, Berlin,

919

2009.

920

ED

918

[56] A. Zebib, D.R. Kassoy, Onset of natural convection in a box watersaturated porous media with large temperature variation, Phys. Fluids(1977),20:4-

922

9. http://dx.doi.org/10.1063/1.861695.

AC

CE

PT

921

54